# ----------------------------------------------------------------------------------------
# SIMULATION PARAMETERS FOR THE PIC-CODE SMILEI
#
# Particle injection from the Xmin and Xmax boundaries
#
# ----------------------------------------------------------------------------------------

import math
import numpy as np

# Mean velocity
mean_velocity = 0.999
# Electron temperature
Te = 0.01
# Ion temperature
Ti = 0.001
# Ion charge
Zi = 1
# Density
n0 = 1
# Debye length
Debye_length = 1. / np.sqrt( n0 / Te + Zi * n0 / Ti )
# Cell length
cell_length = [Debye_length*0.5, Debye_length*0.5]
# Number of patches
number_of_patches =[4, 16]
# Cells per patches (patch shape)
cells_per_patch = [32., 8.]
# Grid length
grid_length = [0.,0.]
for i in range(2):
    grid_length[i] = number_of_patches[i] * cell_length[i] * cells_per_patch[i]
# Number of particles per cell
particles_per_cell = 25
# Position init
position_initialization = 'random'
# Time step
timestep = 0.95/np.sqrt(1./ cell_length[0]**2 + 1./ cell_length[1]**2 )
# Total simulation time
simulation_time = ((1.5 - 0.125)*grid_length[0])/mean_velocity          # duration of the simulation
# Period of output for the diags
diag_every = int(simulation_time / timestep)
# Boundary conditions for particles
particle_boundary_conditions = [["remove", "remove"],["remove", "remove"]]
# Boundary conditions for fields
field_boundary_conditions = [['silver-muller'],['silver-muller']]

Main(
    geometry = "2Dcartesian",
    interpolation_order = 2 ,
    cell_length = cell_length,
    grid_length  = grid_length,
    number_of_patches = number_of_patches,
    timestep = timestep,
    simulation_time = simulation_time,
    EM_boundary_conditions = field_boundary_conditions,
    solve_relativistic_poisson = True,
  
    relativistic_poisson_max_iteration = 500000,

)

LoadBalancing(
	every = 100
)

# Initial plasma shape
fp = trapezoidal(1., xvacuum=grid_length[0]*7.0/8.0                 ,xplateau=grid_length[0]/8.)
def rho(x,y):
    if x<grid_length[0]*7.0/8.0: return 0.
    else: return 1.0*200


Species(
	name = 'pon1',
	position_initialization = 'regular',
	momentum_initialization = 'cold',
	ionization_model = 'none',
	particles_per_cell = particles_per_cell,
	c_part_max = 1.0,
	mass = 1836.0,
	charge = 1.0,
	number_density = rho,
	mean_velocity = [-0.99,0.,0.],
	temperature = [0.0],
	time_frozen = 0.0,
	boundary_conditions = particle_boundary_conditions,
    pusher = "higueracary",
    relativistic_field_initialization=True,
)

ParticleInjector(
    species = 'pon1',
    box_side = 'xmax',
    regular_number=[5,5]
)

DiagFields(
    every = 1,
    fields = ['Rho_pon1']
)

#DiagProbe(
#    every = 5,
#    origin = [0., Main.grid_length[1]/2.],
#    corners = [[Main.grid_length[0],Main.grid_length[1]/2.]],
#    number = [256],
#    fields = ['Ex','Ey','Ez','Rho_pon1','Jx_pon1','Jy_pon1','Jz_pon1']
#)