How to use the fbpic.fields.Fields function in fbpic

Check that the relation div(E) - rho/epsilon_0 is satisfied, with a
    relative precision close to the machine precision (directly in spectral space)

    Parameters
    ----------
    sim: Simulation object
    rho_ions: list of 2d complex arrays (one per mode)
        The density of the ions (which are not explicitly present in the `sim`
        object, since they are motionless)
    """
    # Create a global field object across all subdomains, and copy the fields
    global_Nz, _ = sim.comm.get_Nz_and_iz(
            local=False, with_damp=False, with_guard=False )
    global_zmin, global_zmax = sim.comm.get_zmin_zmax(
            local=False, with_damp=False, with_guard=False )
    global_fld = Fields( global_Nz, global_zmax,
            sim.fld.Nr, sim.fld.rmax, sim.fld.Nm, sim.fld.dt,
            zmin=global_zmin, n_order=sim.fld.n_order, use_cuda=False)
    # Gather the fields of the interpolation grid
    for m in range(sim.fld.Nm):
        # Gather E
        for field in ['Er', 'Et', 'Ez' ]:
            local_array = getattr( sim.fld.interp[m], field )
            gathered_array = sim.comm.gather_grid_array( local_array )
            setattr( global_fld.interp[m], field, gathered_array )
        # Gather rho
        global_fld.interp[m].rho = \
            sim.comm.gather_grid_array( sim.fld.interp[m].rho + rho_ions[m] )

    # Loop over modes and check charge conservation in spectral space
    if sim.comm.rank == 0:
        global_fld.interp2spect('E')

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 fields object
    fld = Fields( Nz, Lz, Nr, Lr, Nm, dt, v_galilean = 0.999*c )
    z0 = Lz/2
    init_fields( 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)
        
    # Get the fields in spectral space
    fld.interp2spect('E')
    fld.interp2spect('B')
    
    #Create moving window object
    mov_win = MovingWindow( fld.interp[0], v= 0. )

    # Loop over the iterations

# Project the charge and currents onto the local subdomain
    sim.deposit( 'rho', exchange=True, species_list=[ptcl],
                    update_spectral=False )
    sim.deposit( 'J', exchange=True, species_list=[ptcl],
                    update_spectral=False )

    # Create a global field object across all subdomains, and copy the sources
    # (Space-charge calculation is a global operation)
    # Note: in the single-proc case, this is also useful in order not to
    # erase the pre-existing E and B field in sim.fld
    global_Nz, _ = sim.comm.get_Nz_and_iz(
                    local=False, with_damp=True, with_guard=False )
    global_zmin, global_zmax = sim.comm.get_zmin_zmax(
                    local=False, with_damp=True, with_guard=False )
    global_fld = Fields( global_Nz, global_zmax,
            sim.fld.Nr, sim.fld.rmax, sim.fld.Nm, sim.fld.dt,
            n_order=sim.fld.n_order, smoother=sim.fld.smoother,
            zmin=global_zmin, use_cuda=False)
    # Gather the sources on the interpolation grid of global_fld
    for m in range(sim.fld.Nm):
        for field in ['Jr', 'Jt', 'Jz', 'rho']:
            local_array = getattr( sim.fld.interp[m], field )
            gathered_array = sim.comm.gather_grid_array(
                                            local_array, with_damp=True )
            setattr( global_fld.interp[m], field, gathered_array )

    # Calculate the space-charge fields on the global grid
    # (For a multi-proc simulation: only performed by the first proc)
    if sim.comm.rank == 0:
        # - Convert the sources to spectral space
        global_fld.interp2spect('rho_prev')

self.boost = None
        # Register time step
        self.dt = dt

        # Initialize the boundary communicator
        cdt_over_dr = c*dt / (rmax/Nr)
        self.comm = BoundaryCommunicator( Nz, zmin, zmax, Nr, rmax, Nm, dt,
            self.v_comoving, self.use_galilean, boundaries, n_order,
            n_guard, n_damp, cdt_over_dr, None, exchange_period,
            use_all_mpi_ranks )
        # Modify domain region
        zmin, zmax, Nz = self.comm.divide_into_domain()
        Nr = self.comm.get_Nr( with_damp=True )
        rmax = self.comm.get_rmax( with_damp=True )
        # Initialize the field structure
        self.fld = Fields( Nz, zmax, Nr, rmax, Nm, dt,
                    n_order=n_order, zmin=zmin,
                    v_comoving=v_comoving,
                    use_pml=self.use_pml,
                    use_galilean=use_galilean,
                    current_correction=current_correction,
                    use_cuda=self.use_cuda,
                    smoother=smoother,
                    # Only create threading buffers when running on CPU
                    create_threading_buffers=(self.use_cuda is False) )

        # Initialize the electrons and the ions
        self.grid_shape = self.fld.interp[0].Ez.shape
        self.particle_shape = particle_shape
        self.ptcl = []
        if n_e is not None:
            # - Initialize the electrons

# (This is done in preparation for gathering among procs)
    saved_Er = []
    saved_Et = []
    for m in range(sim.fld.Nm):
        saved_Er.append( sim.fld.interp[m].Er.copy() )
        sim.fld.interp[m].Er[:,:] = laser_Er[:,:,m]
        saved_Et.append( sim.fld.interp[m].Et.copy() )
        sim.fld.interp[m].Et[:,:] = laser_Et[:,:,m]

    # Create a global field object across all subdomains, and copy the fields
    # (Calculating the self-consistent Ez and B is a global operation)
    global_Nz, _ = sim.comm.get_Nz_and_iz(
                    local=False, with_damp=True, with_guard=False )
    global_zmin, global_zmax = sim.comm.get_zmin_zmax(
                    local=False, with_damp=True, with_guard=False )
    global_fld = Fields( global_Nz, global_zmax,
            sim.fld.Nr, sim.fld.rmax, sim.fld.Nm, sim.fld.dt,
            zmin=global_zmin, n_order=sim.fld.n_order, use_cuda=False)
    # Gather the fields of the interpolation grid
    for m in range(sim.fld.Nm):
        for field in ['Er', 'Et']:
            local_array = getattr( sim.fld.interp[m], field )
            gathered_array = sim.comm.gather_grid_array(
                                local_array, with_damp=True)
            setattr( global_fld.interp[m], field, gathered_array )

    # Now that the (gathered) laser fields are stored in global_fld,
    # copy the saved field back into the local grid
    for m in range(sim.fld.Nm):
        sim.fld.interp[m].Er[:,:] = saved_Er[m]
        sim.fld.interp[m].Et[:,:] = saved_Et[m]