# ----------------------------------------------------------------------------------------
#                     SIMULATION PARAMETERS FOR THE PIC-CODE SMILEI
# ----------------------------------------------------------------------------------------

import math
import numpy as np
from numpy import pi, sqrt

# calculation of electron skin-depth
def esd_cgs(density):
    ne = density*1.1e21
    return 5.31e5 * ne**(-0.5) * 1e4

# Mean velocity
gamma= 2. # relativistic lorentz factor
beta = math.sqrt(1.-1/gamma**2)


ppc = 1

l0 = 2.0*math.pi       # laser wavelength
t0 = l0                # optical cycle

deni = 25.0
chgi = 10.0
dene = deni * chgi 

temp_e = 2.4e-5       
temp_i = 2.4e-5  

# To reproduce the issue, just divide the following by 200      
temp_h = 0.2         #  issue appear: 100 keV for hot electron beam
# temp_h = 0.2/200   #  issue gone:    50  eV for hot electron beam 

reso = 1.0
cell_length = [reso*esd_cgs(dene)*l0, reso*esd_cgs(dene)*l0]
# Number of patches
N_patches = 128
number_of_patches =[N_patches, N_patches]
# Cells per patches (patch shape)
mult = 1.
N_cells = 20
cells_per_patch = [mult*N_cells, mult*N_cells]
# 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]

Tsim = ((1.5 - 0.125)*grid_length[0])/beta
timestep = 0.95/np.sqrt(1./ cell_length[0]**2 + 1./ cell_length[1]**2)
globalEvery = int(Tsim / timestep/ 30)

c = 3.0e8
lambdar = 1.0e-6
wr = 2*pi*c/lambdar                # Normalization frequency

Main(
    geometry = "2Dcartesian",
    interpolation_order = 4 ,
    cell_length  = cell_length,
    grid_length  = grid_length,
    number_of_patches = number_of_patches,
    timestep = timestep,
    simulation_time = Tsim,
    EM_boundary_conditions = [
        ['silver-muller'],
        ["reflective"]
    ],
    random_seed = smilei_mpi_rank,
    reference_angular_frequency_SI = wr,
    print_every                    = 100,

    solve_poisson = False,
    solve_relativistic_poisson = True,
)

Species(
    name = 'ion',
    position_initialization = 'random',
    momentum_initialization = 'mj',
    atomic_number = 79,
    particles_per_cell = ppc,
    mass = 1836.0*197.0,
    charge = chgi,
    number_density = trapezoidal(deni,xvacuum=0.,xplateau=grid_length[0]),
    temperature = [temp_i]*3,
    boundary_conditions = [
        ["reflective"],
        ["reflective"]
    ],
    # time_frozen = 1000000         # this will crash the code, don't know why
)


Species(    
    name = 'eon',
    position_initialization = 'random',
    momentum_initialization = 'mj',
    particles_per_cell = ppc,
    mass = 1.0,
    charge = -1.0,
    number_density = trapezoidal(dene,xvacuum=0.,xplateau=grid_length[0]),
    temperature = [temp_e]*3, 
    boundary_conditions = [
        ["reflective"],
        ["reflective"]
    ], 
)

Species(    
    name = 'colle',
    position_initialization = 'random',
    momentum_initialization = 'cold',
    particles_per_cell = 0,
    mass = 1.0,
    charge = -1.0,
    number_density = trapezoidal(0.0,xvacuum=0.,xplateau=grid_length[0]),
    boundary_conditions = [
        ["reflective"],
        ["reflective"]
    ], 
)

Species(    
    name = 'heb',
    position_initialization = 'random',
    momentum_initialization = 'mj',
    particles_per_cell = 0,
    mass = 1.0,
    charge = -1.0,
    # number_density = trapezoidal(dene,xvacuum=tgtf,xplateau=tgtw),
    number_density = trapezoidal(0.0,xvacuum=0.,xplateau=grid_length[0]),
    temperature = [temp_h]*3, 
    boundary_conditions = [
        ["remove"],
        ["reflective"]
    ], 
)

ParticleInjector(
    name      = "heb_inj",
    species   = "heb",
    box_side  = "xmin",
    position_initialization = "random",
    momentum_initialization = "mj",
    mean_velocity = [beta,0.,0.],
    temperature = [temp_h]*3,
    number_density = trapezoidal(1.0,yvacuum=0.20*grid_length[1],yslope1=0.05*grid_length[1],
                                     yplateau=0.5*grid_length[1],yslope2=0.05*grid_length[1]),
    particles_per_cell = ppc,
)

Collisions(
    species1 = ["heb"],
    species2 = ["ion"],
    # coulomb_log = 5.,
    debug_every = globalEvery,
    ionizing = "colle",
#      nuclear_reaction = [],
)

DiagScalar(every=globalEvery)