From 8eb74ffc67a8cf25424c6ae05e997e6adcea708d Mon Sep 17 00:00:00 2001
From: Francesco Massimo
 <francescomassimo@client-172-18-119-209.eduroam.universite-paris-saclay.fr>
Date: Tue, 21 Nov 2023 11:50:23 +0100
Subject: [PATCH] improve script

---
 InputNamelist.py                              | 232 +++++++++---------
 .../Laser_waist_theory_vs_Smilei.py           |   6 +-
 2 files changed, 116 insertions(+), 122 deletions(-)

diff --git a/InputNamelist.py b/InputNamelist.py
index 20c4b00..90fe017 100755
--- a/InputNamelist.py
+++ b/InputNamelist.py
@@ -78,28 +78,28 @@
     random_seed = smilei_mpi_rank
 )
 
-# ######################### Define the laser pulse
-# 
-# #### laser parameters
-# laser_fwhm        = 25.5*math.sqrt(2)*fs                              # laser FWHM duration in field, i.e. FWHM duration in intensity*sqrt(2)
-# laser_waist       = 12*um                                             # laser waist, conversion from um
-# center_laser      = Lx-1.7*laser_fwhm                                 # laser position at the start of the simulation
-# x_focus           = (center_laser+0.1*laser_fwhm)                     # laser focal plane position
-# a0                = 1.8                                               # laser peak field, normalized by E0 defined above
-# 
-# #### Define a Gaussian bunch with Gaussian temporal envelope
-# LaserEnvelopeGaussianAM(
-#   a0              = a0, 
-#   omega           = (2.*math.pi/lambda0*c)/reference_frequency,       # laser frequency, normalized
-#   focus           = [x_focus,0.],                                     # laser focus, [x,r] position
-#   waist           = laser_waist,                                      # laser waist
-#   time_envelope   = tgaussian(center=center_laser, fwhm=laser_fwhm),  # time profile of the laser pulse
-#   envelope_solver = 'explicit_reduced_dispersion',
-#   Envelope_boundary_conditions     = [ ["reflective"],["PML"] ],
-#   Env_pml_sigma_parameters         = [[0.9 ,2     ],[80.0,2]     ,[80.0,2     ]],
-#   Env_pml_kappa_parameters         = [[1.00,1.00,2],[1.00,1.00,2],[1.00,1.00,2]],
-#   Env_pml_alpha_parameters         = [[0.90,0.90,1],[0.65,0.65,1],[0.65,0.65,1]]
-# )
+######################### Define the laser pulse
+
+#### laser parameters
+#laser_fwhm        = 25.5*math.sqrt(2)*fs                              # laser FWHM duration in field, i.e. FWHM duration in intensity*sqrt(2)
+#laser_waist       = 12*um                                             # laser waist, conversion from um
+#center_laser      = Lx-1.7*laser_fwhm                                 # laser position at the start of the simulation
+#x_focus           = (center_laser+0.1*laser_fwhm)                     # laser focal plane position
+#a0                = 1.8                                               # laser peak field, normalized by E0 defined above
+
+#### Define a Gaussian bunch with Gaussian temporal envelope
+#LaserEnvelopeGaussianAM(
+#  a0              = a0, 
+#  omega           = (2.*math.pi/lambda0*c)/reference_frequency,       # laser frequency, normalized
+#  focus           = [x_focus,0.],                                     # laser focus, [x,r] position
+#  waist           = laser_waist,                                      # laser waist
+#  time_envelope   = tgaussian(center=center_laser, fwhm=laser_fwhm),  # time profile of the laser pulse
+#  envelope_solver = 'explicit_reduced_dispersion',
+#  Envelope_boundary_conditions     = [ ["reflective"],["PML"] ],
+#  Env_pml_sigma_parameters         = [[0.9 ,2     ],[80.0,2]     ,[80.0,2     ]],
+#  Env_pml_kappa_parameters         = [[1.00,1.00,2],[1.00,1.00,2],[1.00,1.00,2]],
+#  Env_pml_alpha_parameters         = [[0.90,0.90,1],[0.65,0.65,1],[0.65,0.65,1]]
+#)
 
 
 ######################### Define a moving window
@@ -109,97 +109,95 @@
     velocity_x = c_normalized,
 )
 
-# ######################### Define the plasma
-# 
-# ##### plasma parameters
-# plasma_plateau_density_1_ov_cm3    = 1.3e18
-# n0 = plasma_plateau_density_1_ov_cm3*1e6/ncrit  # plasma plateau density in units of critical density defined above
-# Radius_plasma = 30.*um                          # Radius of plasma
-# Lramp         = 15*um                           # Plasma density upramp length
-# Lplateau      = 1*mm                            # Length of density plateau
-# Ldownramp     = 15*um                           # Length of density downramp
-# x_begin_upramp  = Lx                            # x coordinate of the start of the density upramp
-# x_begin_plateau = x_begin_upramp+Lramp          # x coordinate of the end of the density upramp / start of density plateau
-# x_end_plateau   = x_begin_plateau+Lplateau      # x coordinate of the end of the density plateau start of the density downramp
-# x_end_downramp  = x_end_plateau+Ldownramp       # x coordinate of the end of the density downramp
-# 
-# #### plasma density profile
-# longitudinal_profile = polygonal(xpoints=[x_begin_upramp,x_begin_plateau,x_end_plateau,x_end_downramp],xvalues=[0.,n0,n0,0.])
-# def plasma_density(x,r):
-# 	profile_r = 0.
-# 	if ((r)**2<Radius_plasma**2):
-# 		profile_r = 1.
-# 	return profile_r*longitudinal_profile(x,r)
-# 
-# ###### define the plasma electrons
-# Species(
-#   name = "plasmaelectrons",
-#   position_initialization = "regular",
-#   momentum_initialization = "cold",
-#   particles_per_cell = 4,
-#   regular_number = [2,2,1],
-#   c_part_max = 1.0,
-#   mass = 1.0,
-#   charge = -1.0,
-#   number_density = plasma_density,
-#   mean_velocity = [0.0, 0.0, 0.0],
-#   temperature = [0.,0.,0.],
-#   pusher = "ponderomotive_borisBTIS3",
-#   time_frozen = 0.0,
-#   boundary_conditions = [
-#      ["remove", "remove"],
-#      ["remove", "remove"],
-#   ],
-# )
+########################## Define the plasma
+#
+###### plasma parameters
+#plasma_plateau_density_1_ov_cm3    = 1.3e18
+#n0 = plasma_plateau_density_1_ov_cm3*1e6/ncrit  # plasma plateau density in units of critical density defined above
+#Radius_plasma = 30.*um                          # Radius of plasma
+#Lramp         = 15*um                           # Plasma density upramp length
+#Lplateau      = 1*mm                            # Length of density plateau
+#Ldownramp     = 15*um                           # Length of density downramp
+#x_begin_upramp  = Lx                            # x coordinate of the start of the density upramp
+#x_begin_plateau = x_begin_upramp+Lramp          # x coordinate of the end of the density upramp / start of density plateau
+#x_end_plateau   = x_begin_plateau+Lplateau      # x coordinate of the end of the density plateau start of the density downramp
+#x_end_downramp  = x_end_plateau+Ldownramp       # x coordinate of the end of the density downramp
+
+##### plasma density profile
+#longitudinal_profile = polygonal(xpoints=[x_begin_upramp,x_begin_plateau,x_end_plateau,x_end_downramp],xvalues=[0.,n0,n0,0.])
+#def plasma_density(x,r):
+#	profile_r = 0.
+#	if ((r)**2<Radius_plasma**2):
+#		profile_r = 1.
+#	return profile_r*longitudinal_profile(x,r)
+
+####### define the plasma electrons
+#Species(
+#  name = "plasmaelectrons",
+#  position_initialization = "regular",
+#  momentum_initialization = "cold",
+#  particles_per_cell = 4,
+#  regular_number = [2,2,1],
+#  mass = 1.0,
+#  charge = -1.0,
+#  number_density = plasma_density,
+#  mean_velocity = [0.0, 0.0, 0.0],
+#  temperature = [0.,0.,0.],
+#  pusher = "ponderomotive_borisBTIS3",
+#  time_frozen = 0.0,
+#  boundary_conditions = [
+#     ["remove", "remove"],
+#     ["remove", "remove"],
+#  ],
+#)
 
 
-# ######################## Define the electron bunch
-# 
-# ##### electron bunch parameters
-# Q_bunch                    = -15*pC                          # Total charge of the electron bunch
-# sigma_x                    = 1.5*um                          # initial longitudinal rms size
-# sigma_r                    = 2*um                            # initial transverse/radial rms size (cylindrical symmetry)
-# bunch_energy_spread        = 0.01                            # initial rms energy spread / average energy (not in percent)
-# bunch_normalized_emittance = 3.*mm_mrad                      # initial rms emittance, same emittance for both transverse planes
-# delay_behind_laser         = 18.5*um                         # distance between center_laser and center_bunch
-# center_bunch               = center_laser-delay_behind_laser # initial position of the electron bunch in the window   
-# gamma_bunch                = 200.                            # initial relativistic Lorentz factor of the bunch
-# 
-# npart                      = 50000                           # number of computational macro-particles to model the electron bunch 
-# normalized_species_charge  = -1                              # For electrons
-# Q_part                     = Q_bunch/npart                   # charge for every macroparticle in the electron bunch
-# weight                     = Q_part/((c/omega0)**3*ncrit*normalized_species_charge)
-# 
-# #### initialize the bunch using numpy arrays
-# #### the bunch will have npart particles, so an array of npart elements is used to define the x coordinate of each particle and so on ...
-# array_position = np.zeros((4,npart))                         # positions x,y,z, weight
-# array_momentum = np.zeros((3,npart))                         # momenta x,y,z
-# 
-# #### The electron bunch is supposed at waist. To make it convergent/divergent, transport matrices can be used
-# array_position[0,:] = np.random.normal(loc=center_bunch, scale=sigma_x, size=npart)                        # generate random number from gaussian distribution for x position
-# array_position[1,:] = np.random.normal(loc=0., scale=sigma_r, size=npart)                                  # generate random number from gaussian distribution for y position
-# array_position[2,:] = np.random.normal(loc=0., scale=sigma_r, size=npart)                                  # generate random number from gaussian distribution for z position
-# array_momentum[0,:] = np.random.normal(loc=gamma_bunch, scale=bunch_energy_spread*gamma_bunch, size=npart) # generate random number from gaussian distribution for px position
-# array_momentum[1,:] = np.random.normal(loc=0., scale=bunch_normalized_emittance/sigma_r, size=npart)       # generate random number from gaussian distribution for py position
-# array_momentum[2,:] = np.random.normal(loc=0., scale=bunch_normalized_emittance/sigma_r, size=npart)       # generate random number from gaussian distribution for pz position
-# 
-# array_position[3,:] = np.multiply(np.ones(npart),weight)
-
-# #### define the electron bunch
-# Species( 
-#   name = "electronbunch",
-#   position_initialization = array_position,
-#   momentum_initialization = array_momentum,
-#   c_part_max = 1.0,
-#   mass = 1.0,
-#   charge = -1.0,
-#   relativistic_field_initialization = True,
-#   pusher = "ponderomotive_borisBTIS3", 
-#   boundary_conditions = [
-#   	["remove", "remove"],
-#   	["remove", "remove"], 
-#   ],
-# )
+######################## Define the electron bunch
+
+###### electron bunch parameters
+#Q_bunch                    = -15*pC                          # Total charge of the electron bunch
+#sigma_x                    = 1.5*um                          # initial longitudinal rms size
+#sigma_r                    = 2*um                            # initial transverse/radial rms size (cylindrical symmetry)
+#bunch_energy_spread        = 0.01                            # initial rms energy spread / average energy (not in percent)
+#bunch_normalized_emittance = 3.*mm_mrad                      # initial rms emittance, same emittance for both transverse planes
+#delay_behind_laser         = 18.5*um                         # distance between center_laser and center_bunch
+#center_bunch               = center_laser-delay_behind_laser # initial position of the electron bunch in the window   
+#gamma_bunch                = 200.                            # initial relativistic Lorentz factor of the bunch
+
+#npart                      = 50000                           # number of computational macro-particles to model the electron bunch 
+#normalized_species_charge  = -1                              # For electrons
+#Q_part                     = Q_bunch/npart                   # charge for every macroparticle in the electron bunch
+#weight                     = Q_part/((c/omega0)**3*ncrit*normalized_species_charge)
+
+##### initialize the bunch using numpy arrays
+##### the bunch will have npart particles, so an array of npart elements is used to define the x coordinate of each particle and so on ...
+#array_position = np.zeros((4,npart))                         # positions x,y,z, weight
+#array_momentum = np.zeros((3,npart))                         # momenta x,y,z
+
+##### The electron bunch is supposed at waist. To make it convergent/divergent, transport matrices can be used
+#array_position[0,:] = np.random.normal(loc=center_bunch, scale=sigma_x, size=npart)                        # generate random number from gaussian distribution for x position
+#array_position[1,:] = np.random.normal(loc=0., scale=sigma_r, size=npart)                                  # generate random number from gaussian distribution for y position
+#array_position[2,:] = np.random.normal(loc=0., scale=sigma_r, size=npart)                                  # generate random number from gaussian distribution for z position
+#array_momentum[0,:] = np.random.normal(loc=gamma_bunch, scale=bunch_energy_spread*gamma_bunch, size=npart) # generate random number from gaussian distribution for px position
+#array_momentum[1,:] = np.random.normal(loc=0., scale=bunch_normalized_emittance/sigma_r, size=npart)       # generate random number from gaussian distribution for py position
+#array_momentum[2,:] = np.random.normal(loc=0., scale=bunch_normalized_emittance/sigma_r, size=npart)       # generate random number from gaussian distribution for pz position
+
+#array_position[3,:] = np.multiply(np.ones(npart),weight)
+
+##### define the electron bunch
+#Species( 
+#  name = "electronbunch",
+#  position_initialization = array_position,
+#  momentum_initialization = array_momentum,
+#  mass = 1.0,
+#  charge = -1.0,
+#  relativistic_field_initialization = True,
+#  pusher = "ponderomotive_borisBTIS3", 
+#  boundary_conditions = [
+#  	["remove", "remove"],
+#  	["remove", "remove"], 
+#  ],
+#)
 
 
 ######################### Current filter
@@ -230,12 +228,12 @@
     fields = ['Ex','Ey','Rho','Env_A_abs','Env_E_abs','Bz','BzBTIS3']
 )
 
-# ##### Diagnostic for the electron bunch macro-particles
-# DiagTrackParticles(
-#   species = "electronbunch",
-#   every = 200,
-#   attributes = ["x", "y", "z", "px", "py", "pz", "w"]
-# )
+##### Diagnostic for the electron bunch macro-particles
+#DiagTrackParticles(
+#  species = "electronbunch",
+#  every = 200,
+#  attributes = ["x", "y", "z", "px", "py", "pz", "w"]
+#)
 
 
 ######################### Load balancing (for parallelization)                                                                                                                                                     
diff --git a/Postprocessing_Scripts/Laser_waist_theory_vs_Smilei.py b/Postprocessing_Scripts/Laser_waist_theory_vs_Smilei.py
index dac37a6..0997940 100755
--- a/Postprocessing_Scripts/Laser_waist_theory_vs_Smilei.py
+++ b/Postprocessing_Scripts/Laser_waist_theory_vs_Smilei.py
@@ -12,10 +12,6 @@
 lambda0           = S.namelist.lambda0*1e6                        # laser wavelength, um 
 conversion_factor = lambda0/(2.*math.pi)
 
-
-####### load simulation
-S    = happi.Open(".")
-
 ####### create grid
 dx   = S.namelist.Main.cell_length[0]                             # resolution on the x direction
 nx   = S.namelist.Main.number_of_cells[0]                         # number of grid points in the x direction
@@ -62,7 +58,7 @@
         half_waist_simulated.append(  half_waist )
                  
         # analytical waist, Rayleigh formula 
-        waist_analytical.append(waist0*math.sqrt(1.+((iter*dt-S.namelist.x_focus)/Zr)**2))
+        waist_analytical.append(waist0*math.sqrt(1.+((S.namelist.Lx+iter*dt-S.namelist.x_focus)/Zr)**2))
 
 waist_analytical= np.asarray(waist_analytical)                
 waist_simulated = 2.*np.asarray(half_waist_simulated)