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

import math

patchesx = 2**(patchx_int-1)
patchesy = 2**(patchy_int-1)
#l0 = 2.0*math.pi	# laser wavelength
#t0 = l0		# optical cicle
ncell = 160
dx = 0.2
Lsim = [ncell*dx,ncell*dx]		#length of the simulation
#Tsim = 30.*t0		# duration of the simulation
#resx = 50.		# nb of cells in on laser wavelength
#rest = 75.		# time of timestep in one optical cycle 

# Lr = c / wr
# for Lr = 1e-6; wr = 3e14, Tr = 3.333e-15 (three femtoseconds per 'unit time')
# nr = (0.0003149*wr*wr) = 2.8344e25
# so to get ne = 1e24 we need a faction of the ref density.
# ne = 1e24 / 2.8344e25 = 0.03528

ne = 0.03528
ppc = 4
ppc_core = 128

def circleJ(max, radius=0, centerx=0, centery=0):
    global Main
    if len(Main)==0:
        raise Exception("circleJ profile has been defined before `Main()`")
    def fxy(x,y):
        if ((x-centerx)**2+(y-centery)**2)>radius**2: return 0.0
        elif ((x-centerx)**2+(y-centery)**2)<radius**2: return max
        else: return 0.0
    f=fxy
    return f

def antiCircleJ(max, radius=0, centerx=0, centery=0):
    global Main
    if len(Main)==0:
        raise Exception("circleJ profile has been defined before `Main()`")
    def fxy(x,y):
        if ((x-centerx)**2+(y-centery)**2)>radius**2: return max
        elif ((x-centerx)**2+(y-centery)**2)<radius**2: return 0.0
        else: return 0.0
    f=fxy
    return f

Main(
    geometry = "2d3v",
    
    interpolation_order = 2 ,
    
    cell_length = [dx,dx],
    sim_length  = Lsim,

    number_of_patches = [ patchesx , patchesy ],
    
    timestep = 0.95*dx/math.sqrt(2),
    #timestep = (0.95*dx)/2.0,
    #timestep = 8.35,
    #number_of_timesteps = 4400000, # ~4ns 
    number_of_timesteps = 2000, # ~4ns 
    #number_of_timesteps = 10000, # ~0.1 / 3 ns 
     
    bc_em_type_x = ['silver-muller'],
    bc_em_type_y = ['silver-muller'],
    
    random_seed = 0,

#    referenceAngularFrequency_SI = 3.0e14
)

core = circleJ(ne, radius=1.5, centerx=(ncell*dx)/2, centery=(ncell*dx)/2)
clad = antiCircleJ(ne, radius=1.5, centerx=(ncell*dx)/2, centery=(ncell*dx)/2) 

Species(
	species_type = 'elecCore',
	initPosition_type = 'random',
	initMomentum_type = 'mj',
#	ionization_model = 'none',
	n_part_per_cell = ppc_core,
	c_part_max = 1.0,
	mass = 1.0,
	charge = 1.0,
	nb_density = core,
	mean_velocity = [0.,0.,0.],
	temperature = [0.03718],
	time_frozen = 0.0,
        bc_part_type_xmin  = 'refl',
        bc_part_type_xmax  = 'refl',
        bc_part_type_ymin  = 'refl',
        bc_part_type_ymax  = 'refl'
)

Species(
	species_type = 'ionCore',
	initPosition_type = 'random',
	initMomentum_type = 'mj',
#	ionization_model = 'none',
	n_part_per_cell = ppc_core,
	c_part_max = 1.0,
	mass = 1836.0,
	charge = -1.0,
	nb_density = core,
	mean_velocity = [0.,0.,0.],
	temperature = [0.03718],
	time_frozen = 0.0,
	bc_part_type_xmin  = 'refl',
	bc_part_type_xmax  = 'refl',
	bc_part_type_ymin  = 'refl',
	bc_part_type_ymax  = 'refl'
)

#Species(
#	species_type = 'hydrogen',
#        atomic_number = 1.0,
##        ionization_model = 'tunnel',
##        ionization_electrons = 'elecIon',
#	initPosition_type = 'random',
#	initMomentum_type = 'mj',
#	n_part_per_cell = ppc,
#	c_part_max = 1.0,
#	mass = 1836.0,
#	charge = 0.0,
#	nb_density = clad,
#	mean_velocity = [0.,0.,0.],
#	temperature = [0.0000004385],
#	time_frozen = 0.0,
#	bc_part_type_xmin  = 'refl',
#	bc_part_type_xmax  = 'refl',
#	bc_part_type_ymin  = 'refl',
#	bc_part_type_ymax  = 'refl'
#)

Species(
	species_type = 'elecIon',
	initPosition_type = 'random',
	initMomentum_type = 'mj',
#	ionization_model = 'none',
	n_part_per_cell = ppc,
	c_part_max = 1.0,
	mass = 1.0,
	charge = 1.0,
	nb_density = 0.0,
	mean_velocity = [0.,0.,0.],
	time_frozen = 0.0,
	bc_part_type_xmin  = 'refl',
	bc_part_type_xmax  = 'refl',
	bc_part_type_ymin  = 'refl',
	bc_part_type_ymax  = 'refl'
)

globalEvery = 2000.0

#Collisions(
#    species1 = ["elecCore",  "elecIon"],
#    species2 = ["hydrogen"],
#    coulomb_log = 0.0,
#    debug_every = globalEvery,
#    ionizing = True,
#)

DiagScalar(every=globalEvery)

DiagFields(
    every = globalEvery,
#    fields = ['Ex','Ey','Ez','Bx','By','Bz','Rho_elecCore','Rho_ionCore','Rho_hydrogen','Rho_elecIon']
#    fields = ['Ex','Ey','Rho_elecCore','Rho_ionCore','Rho_hydrogen','Rho_elecIon']
    fields = ['Ex','Ey','Rho_elecCore','Rho_ionCore']
)

#DiagParticles(
#
#    output = "density",
#    every = globalEvery,
#    time_average = 1,
#    species = ["elecCore"],
#    axes = [ ["x",    0.,    ncell*dx,    ncell],
#             ["y",    0.,    ncell*dx,    ncell] ]
#)