How to use the fbpic.main.Simulation function in fbpic

def simulate_periodic_plasma_wave( particle_shape, show=False ):
    "Simulate a periodic plasma wave and check its fields"

    # Initialization of the simulation object
    sim = Simulation( Nz, zmax, Nr, rmax, Nm, dt,
                  p_zmin, p_zmax, p_rmin, p_rmax, p_nz, p_nr,
                  p_nt, n_e, n_order=n_order, use_cuda=use_cuda,
                  particle_shape=particle_shape )

    # Save the initial density in spectral space, and consider it
    # to be the density of the (uninitialized) ions
    sim.deposit('rho_prev', exchange=True)
    sim.fld.spect2interp('rho_prev')
    rho_ions = [ ]
    for m in range(len(sim.fld.interp)):
        rho_ions.append( -sim.fld.interp[m].rho.copy() )

    # Impart velocities to the electrons
    # (The electrons are initially homogeneous, but have an
    # intial non-zero velocity that develops into a plasma wave)
    impart_momenta( sim.ptcl[0], epsilons, k0, w0, wp )

Returns
    -------
    A dictionary containing :
    - 'E' : 1d array containing the values of the amplitude
    - 'w' : 1d array containing the values of waist
    - 'fld' : the Fields object at the end of the simulation.
    - 'z_centroid' : 1d array containing the values of the centroid position
    """

    # Choose the timestep, so that the simulation advances by
    # one cell at every timestep
    dt = Lz*1./Nz * 1./c
    
    # Initialize the simulation object
    sim = Simulation( Nz, Lz, Nr, Lr, Nm, dt, p_zmin=0, p_zmax=0,
                    p_rmin=0, p_rmax=0, p_nz=2, p_nr=2, p_nt=2, n_e=0.,
                    use_cuda=True, boundaries='open' )
    # Remove the particles
    sim.ptcl = []
    # Create moving window object
    sim.set_moving_window( v=c )

    # Initialize the laser fields
    z0 = Lz/2
    init_fields( sim.fld, w0, ctau, k0, z0, E0, m )

    # Create the arrays to get the waist and amplitude
    w = np.zeros(Nt)
    E = np.zeros(Nt)
    z_center = np.zeros(Nt)

v_comoving : float
        Velocity at which the currents are assumed to move

    use_galilean: bool
        Whether to use a galilean frame that moves at the speed v_comoving

    Returns
    -------
    A dictionary containing :
    - 'E' : 1d array containing the values of the amplitude
    - 'w' : 1d array containing the values of waist
    - 'fld' : the Fields object at the end of the simulation.
    """

    # Initialize the simulation object
    sim = Simulation( Nz, zmax, Nr, Lr, Nm, dt,
                    n_order=n_order, zmin=zmin, use_cuda=use_cuda,
                    boundaries=boundaries, v_comoving=v_comoving,
                    exchange_period = 1, use_galilean=use_galilean )
    # Remove the particles
    sim.ptcl = []
    # Set the moving window object
    if v_window !=0:
        sim.set_moving_window( v=v_window )

    # Initialize the laser fields
    z0 = (zmax+zmin)/2
    init_fields( sim, w0, ctau, k0, z0, zf, E0, m )

    # Create the arrays to get the waist and amplitude
    w = np.zeros(N_diag)
    E = np.zeros(N_diag)

N_step = 300     # Number of iterations to perform
N_show = 300     # Number of timestep between every plot

# The particles
gamma0 = 25.
p_zmin = 15.e-6  # Position of the beginning of the bunch (meters)
p_zmax = 25.e-6  # Position of the end of the bunch (meters)
p_rmin = 0.      # Minimal radial position of the bunch (meters)
p_rmax = 5.e-6   # Maximal radial position of the bunch (meters)
n_e = 4.e18*1.e6 # Density (electrons.meters^-3)
p_nz = 2         # Number of particles per cell along z
p_nr = 2         # Number of particles per cell along r
p_nt = 4         # Number of particles per cell along theta

# Initialize the simulation object
sim = Simulation( Nz, zmax, Nr, rmax, Nm, dt,
    p_zmin, p_zmax, p_rmin, p_rmax, p_nz, p_nr, p_nt, n_e, n_order=n_order,
    v_comoving=-0.999999*c, use_galilean=True )

# Suppress the particles that were intialized by default and add the bunch
sim.ptcl = [ ]
add_elec_bunch( sim, gamma0, n_e, p_zmin, p_zmax,
            p_rmin, p_rmax, direction='backward' )


# Show the initial fields
plt.figure(0)
sim.fld.interp[0].show('Ez')
plt.figure(1)
sim.fld.interp[0].show('Er')
plt.show()
print( 'Done' )

v_comoving : float
        Velocity at which the currents are assumed to move

    use_galilean: bool
        Whether to use a galilean frame that moves at the speed v_comoving

    Returns
    -------
    A dictionary containing :
    - 'a' : 1d array containing the values of the amplitude
    - 'w' : 1d array containing the values of waist
    - 'fld' : the Fields object at the end of the simulation.
    """

    # Initialize the simulation object
    sim = Simulation( Nz, zmax, Nr, Lr, Nm, dt, p_zmin=0, p_zmax=0,
                    p_rmin=0, p_rmax=0, p_nz=2, p_nr=2, p_nt=2, n_e=0.,
                    n_order=n_order, zmin=zmin, use_cuda=use_cuda,
                    boundaries=boundaries, v_comoving=v_comoving,
                    exchange_period = 1, use_galilean=use_galilean,
                    use_envelope = True )
    # Remove the particles
    sim.ptcl = []
    # Set the moving window object
    if v_window !=0:
        sim.set_moving_window( v=v_window )

    # Initialize the laser fields
    z0 = (zmax+zmin)/2
    init_fields( sim, w0, ctau, k0, z0, zf, a0, m )

    # Create the arrays to get the waist and amplitude

interpolation='nearest')
    plt.xlabel('z')
    plt.ylabel('r')
    cb = plt.colorbar()
    cb.set_label('Imaginary part')

    plt.show()

Nm = 1
dt = (zmax-zmin)*1./c/Nz
dt = 0.06e-6/c
print(c*dt)
print(L_prop)
print(L_prop / c / dt)

sim = Simulation( Nz, zmax, Nr, rmax, Nm, dt,
    p_zmin=p_zmin, p_zmax=p_zmax, p_rmin=p_rmin, p_rmax=p_rmax, p_nz=p_nz,
    p_nr=p_nr, p_nt=Nm+1, n_e=n_e, n_order=n_order, zmin=zmin,
    boundaries='open', v_comoving=0.999*c, use_galilean=True, initialize_ions=True,
    use_cuda=use_cuda, use_envelope=True )
sim.set_moving_window(v=c)
tau = ctau/c
lambda0 = 2*np.pi/k0
# Create the relevant laser profile
z0 = 0
sim.diags = [ FieldDiagnostic( diag_period, sim.fld, comm=sim.comm, fieldtypes = ["a", "rho", "E"] ),
        ParticleDiagnostic(diag_period, {'electrons': sim.ptcl[0]}, sim.comm )]

profile = GaussianLaser( a0=a0, waist=w0, tau=tau,
            lambda0=lambda0, z0=z0, zf=zf )

add_laser_pulse( sim, profile, method = 'direct_envelope' )

for scheme in [ 'standard', 'galilean', 'pseudo-galilean']:

        # Choose the correct parameters for the scheme
        if scheme == 'standard':
            v_comoving = 0.
            use_galilean = False
        else:
            v_comoving = 0.9999*c
            if scheme == 'galilean':
                use_galilean = True
            else:
                use_galilean = False

        # Initialize the simulation object
        sim = Simulation( Nz, zmax, Nr, rmax, Nm, dt,
            p_zmin, p_zmax, p_rmin, p_rmax, p_nz, p_nr, p_nt, n_e,
            zmin=zmin, initialize_ions=True,
            v_comoving=v_comoving, use_galilean=use_galilean,
            boundaries='periodic', use_cuda=use_cuda )

        # Give a relativistic velocity to the particle, with some noise
        sim.ptcl[0].uz[:] = uz_m
        sim.ptcl[0].inv_gamma[:] = 1./np.sqrt( 1 + sim.ptcl[0].uz**2 )
        sim.ptcl[1].uz[:] = uz_m
        sim.ptcl[1].inv_gamma[:] = 1./np.sqrt( 1 + sim.ptcl[1].uz**2 )

        # Perform the simulation;
        # record the rms electric field every 50 timestep
        Er_rms = np.zeros(int(N_step/30)+1)
        t = np.zeros(int(N_step/30+1))
        Er_rms[0] = get_Er_rms(sim)

# The interaction length of the simulation (meters)
L_interact = 50.e-6 # increase to simulate longer distance!
# Interaction time (seconds) (to calculate number of PIC iterations)
T_interact = ( L_interact + (zmax-zmin) ) / v_window
# (i.e. the time it takes for the moving window to slide across the plasma)

# ---------------------------
# Carrying out the simulation
# ---------------------------

# NB: The code below is only executed when running the script,
# (`python lwfa_script.py`), but not when importing it (`import lwfa_script`).
if __name__ == '__main__':

    # Initialize the simulation object
    sim = Simulation( Nz, zmax, Nr, rmax, Nm, dt, zmin=zmin,
        n_order=n_order, use_cuda=use_cuda,
        boundaries={'z':'open', 'r':'reflective'})
        # 'r': 'open' can also be used, but is more computationally expensive

    # Create the plasma electrons
    elec = sim.add_new_species( q=-e, m=m_e, n=n_e,
        dens_func=dens_func, p_zmin=p_zmin, p_zmax=p_zmax, p_rmax=p_rmax,
        p_nz=p_nz, p_nr=p_nr, p_nt=p_nt )

    # Load initial fields
    # Add a laser to the fields of the simulation
    add_laser( sim, a0, w0, ctau, z0 )

    if use_restart is False:
        # Track electrons if required (species 0 correspond to the electrons)
        if track_electrons: