Docs: Thomson Parabola Spectrometer example

Docs/source/refs.bib

@@ -441,3 +441,18 @@ @article{Vranic2015
   issn = {0010-4655},
   doi = {https://doi.org/10.1016/j.cpc.2015.01.020},
 }
+
+@article{Rhee1987,
+    author = {Rhee, M. J. and Schneider, R. F. and Weidman, D. J.},
+    title = "{Simple time‐resolving Thomson spectrometer}",
+    journal = {Review of Scientific Instruments},
+    volume = {58},
+    number = {2},
+    pages = {240-244},
+    year = {1987},
+    month = {02},
+    issn = {0034-6748},
+    doi = {10.1063/1.1139314},
+    url = {https://doi.org/10.1063/1.1139314},
+    eprint = {https://pubs.aip.org/aip/rsi/article-pdf/58/2/240/19154912/240\_1\_online.pdf},
+}

Docs/source/usage/examples.rst

@@ -44,6 +44,7 @@ Particle Accelerator & Beam Physics
 
    examples/gaussian_beam/README.rst
    examples/beam-beam_collision/README.rst
+   examples/thomson_parabola_spectrometer/README.rst
 
 
 High Energy Astrophysical Plasma Physics

Examples/Physics_applications/thomson_parabola_spectrometer/README.rst

@@ -0,0 +1,54 @@
+.. _examples-thomson_parabola_spectrometer:
+
+Thomson Parabola Spectrometer
+=============================
+
+This example simulates a Thomson parabola spectrometer (TPS) :cite:t:`ex-Rhee1987`.
+
+A TPS is a type of detector that separates incoming ions according to their charge-to-mass ratio (:math:`q/m`) and initial velocity (hence energy :math:`E_0 = 1/2 m v_0^2` if we assume non-relativistic dynamics).
+TPSs are often used in laser-driven ion acceleration experiments, where different ion species are accelerated at once. To mimic this, we initialize a point-like source of 3 different ion species with different :math:`q/m` and :math:`E_0` (i.e. all ions have the same initial position, representative of a pinhole).
+
+The ions propagate along :math:`z` through 4 subsequent regions:
+
+ - a vacuum region, the distance between the pinhole and the TPS (0.1 m)
+ - a region of constant electric field along :math:`x`, (0.19 m, 1e5 V/m)
+ - a region of constant magnetic field along :math:`x`, (0.872 T, 0.12 m)
+ - a vacuum region, the distance between the TPS and the screen of the detector (0.2 m)
+
+The initial particle velocity :math:`v_0` is sampled from a uniform distribution in the range :math:`[v_{min}, v_{max}]` where :math:`v_{min} = \sqrt{E_{max}/m}`, :math:`v_{max} = \sqrt{2E_{max}/m}`, and :math:`E_{max}` is an input parameter for each species. We assume zero transverse momentum.
+
+The detector is modeled using a ``BoundaryScrapingDiagnostic`` at the upper :math:`z` boundary of the domain, which stores the attributes of the particles when they exit the simulation box from the corresponding edge. Note that the transverse box size is large enough such that all particles exit the domain from the upper :math:`z` side.
+
+Run
+---
+
+The PICMI input file is not available for this example yet.
+
+For `MPI-parallel <https://www.mpi-forum.org>`__ runs, prefix these lines with ``mpiexec -n 4 ...`` or ``srun -n 4 ...``, depending on the system.
+
+.. literalinclude:: inputs
+   :language: ini
+   :caption: You can copy this file from ``Examples/Physics_applications/thomson_parabola_spectrometer/inputs``.
+
+Visualize
+---------
+
+This figure below shows the ion trajectories starting from the pinhole (black star), entering the E and B field regions (purple box), up to the detector (gray plane).
+The colors represent the different species: protons in blue, C :sup:`+4` in red, and C :sup:`+6` in green.
+The particles are accelerated and deflected through the TPS.
+
+.. figure:: https://gist.github.com/assets/17280419/3e45e5aa-d1fc-46e3-aa24-d9e0d6a74d1a
+   :alt: Ion trajectories through a synthetic TPS.
+   :width: 100%
+
+In our simulation, the virtual detector stores all the particle data once entering it (i.e. exiting the simulation box).
+The figure below shows the ions colored according to their species (same as above) and shaded according to their initial energy.
+The :math:`x` coordinate represents the electric deflection, while :math:`y` the magnetic deflection.
+
+.. figure:: https://gist.github.com/assets/17280419/4dd1adb7-b4ab-481d-bc24-8a7ca51471d9
+   :alt: Synthetic TPS screen.
+   :width: 100%
+
+.. literalinclude:: plot_detector.py
+   :language: ini
+   :caption: You can copy this file from ``Examples/Physics_applications/thomson_parabola_spectrometer/plot_detector.py``.

Examples/Physics_applications/thomson_parabola_spectrometer/inputs.txt

@@ -0,0 +1,135 @@
+##############
+#### CONSTANTS
+##############
+my_constants.MeV = 1e6*q_e
+my_constants.d1 = 0.1 # m
+my_constants.d2 = 0.19 # m
+my_constants.d3 = 0.12 # m
+my_constants.d4 = 0.2 # m
+my_constants.E0 = 1e5 # V/m
+my_constants.B0 = 0.872 # T
+my_constants.xmin = -0.4
+my_constants.xmax =  0.4
+my_constants.ymin = -0.4
+my_constants.ymax =  0.4
+my_constants.zmin= -1e-3 # m
+my_constants.zmax = d1+d2+d3+d4
+my_constants.N_real_particles = 1e3
+my_constants.N_macro_particles = 1e3
+my_constants.Emax_hydrogen1_1 = 40*MeV
+my_constants.Emax_carbon12_6 = 10*MeV
+my_constants.vz = sqrt(2*1*MeV/(12*m_p))
+my_constants.max_steps = 400
+my_constants.max_time = (-zmin+d1+d2+d3+d4) / vz
+my_constants.dt = max_time / max_steps
+
+#############
+#### NUMERICS
+#############
+algo.particle_shape = 1
+algo.maxwell_solver = none
+algo.particle_pusher = boris
+amr.max_level = 0
+warpx.verbose = 1
+
+########
+#### BOX
+########
+amr.n_cell = 8 8 8
+geometry.dims = 3
+geometry.prob_hi = xmax ymax zmax
+geometry.prob_lo = xmin ymin zmin
+
+#########
+#### TIME
+#########
+stop_time = max_time
+warpx.const_dt = dt
+
+#############
+#### BOUNDARY
+#############
+boundary.particle_hi = absorbing absorbing absorbing
+boundary.particle_lo = absorbing absorbing absorbing
+
+##############
+#### PARTICLES
+##############
+particles.species_names = hydrogen1_1 carbon12_6
+
+hydrogen1_1.charge = q_e
+hydrogen1_1.initialize_self_fields = 0
+hydrogen1_1.injection_style = gaussian_beam
+hydrogen1_1.mass = m_p
+hydrogen1_1.momentum_distribution_type = uniform
+hydrogen1_1.npart = N_macro_particles
+hydrogen1_1.q_tot = N_real_particles*q_e
+hydrogen1_1.ux_min = 0
+hydrogen1_1.uy_min = 0
+hydrogen1_1.uz_min = sqrt(2*Emax_hydrogen1_1/m_p)/clight*0.5
+hydrogen1_1.ux_max = 0
+hydrogen1_1.uy_max = 0
+hydrogen1_1.uz_max = sqrt(2*Emax_hydrogen1_1/m_p)/clight
+hydrogen1_1.x_m = 0
+hydrogen1_1.x_rms = 0
+hydrogen1_1.y_m = 0
+hydrogen1_1.y_rms = 0
+hydrogen1_1.z_m = 0
+hydrogen1_1.z_rms = 0
+hydrogen1_1.do_not_gather = 1
+hydrogen1_1.do_not_deposit = 1
+
+carbon12_6.charge = 6*q_e
+carbon12_6.initialize_self_fields = 0
+carbon12_6.injection_style = gaussian_beam
+carbon12_6.mass = 12*m_p
+carbon12_6.momentum_distribution_type = uniform
+carbon12_6.npart = N_macro_particles
+carbon12_6.q_tot = N_real_particles*6*q_e
+carbon12_6.ux_min = 0
+carbon12_6.uy_min = 0
+carbon12_6.uz_min = sqrt(2*Emax_carbon12_6/(12*m_p))/clight*0.5
+carbon12_6.ux_max = 0
+carbon12_6.uy_max = 0
+carbon12_6.uz_max = sqrt(2*Emax_carbon12_6/(12*m_p))/clight
+carbon12_6.x_m = 0
+carbon12_6.x_rms = 0
+carbon12_6.y_m = 0
+carbon12_6.y_rms = 0
+carbon12_6.z_m = 0
+carbon12_6.z_rms = 0
+carbon12_6.do_not_gather = 1
+carbon12_6.do_not_deposit = 1
+
+###########
+#### FIELDS
+###########
+particles.E_ext_particle_init_style = parse_E_ext_particle_function
+particles.Ex_external_particle_function(x,y,z,t) = "E0*(z>d1)*(z<(d1+d2))"
+particles.Ey_external_particle_function(x,y,z,t) = 0
+particles.Ez_external_particle_function(x,y,z,t) = 0
+
+particles.B_ext_particle_init_style = parse_B_ext_particle_function
+particles.Bx_external_particle_function(x,y,z,t) = "B0*(z>d1+d2)*(z<(d1+d2+d3))"
+particles.By_external_particle_function(x,y,z,t) = 0
+particles.Bz_external_particle_function(x,y,z,t) = 0
+
+################
+#### DIAGNOSTICS
+################
+diagnostics.diags_names = diag1 diag2
+
+diag1.diag_type = Full
+diag1.fields_to_plot = none
+diag1.format = openpmd
+diag1.intervals = 1
+diag1.openpmd_backend = h5
+diag1.write_species = 1
+diag1.species = hydrogen1_1 carbon12_6
+
+diag2.diag_type = BoundaryScraping
+diag2.format = openpmd
+diag2.intervals = 1
+diag2.openpmd_backend = h5
+hydrogen1_1.save_particles_at_zhi = 1
+carbon12_6.save_particles_at_zhi = 1

Examples/Physics_applications/thomson_parabola_spectrometer/plot_detector.py

@@ -0,0 +1,77 @@
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.cm import ScalarMappable
+from openpmd_viewer import OpenPMDTimeSeries
+from scipy.constants import c, eV, m_p
+
+mpl.use('Agg')
+mpl.rcParams.update({'font.size': 18})
+
+MeV=1e6*eV
+
+# some parameters from the inputfile
+Emax_hydrogen1_1 = 40*MeV
+Emax_carbon12_6 = 20*MeV
+Emax_carbon12_5 = 20*MeV
+
+# open the BoundaryScrapingDiagnostic that represents the detector
+series = OpenPMDTimeSeries('./diags/screen/particles_at_zhi/')
+# open the Full diagnostic at time zero
+series0 = OpenPMDTimeSeries('./diags/diag0/')
+# we use the data at time 0 to retrieve the initial energy
+# of all the particles the boundary
+
+# timesteps and real times
+it = series.iterations
+time = series.t # s
+N_iterations = len(it)
+
+# list of species names
+species = series.avail_species
+N_species = len(species)
+
+fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(10,8), dpi=300)
+
+# some stuff for plotting
+vmin=0
+vmax=50
+cmap=['Reds', 'Greens', 'Blues']
+
+# loop through the species
+for s in range(N_species):
+    print(species[s])
+
+    # arrays of positions and energies
+    X, Y, E = [], [], []
+    for i in range(N_iterations):
+        # get particles at detector location
+        x,y,z,ids = series.get_particle( ['x','y','z','id'],
+                                         iteration=it[i],
+                                         species=species[s],
+                                         plot=False )
+        # get particles at initialization
+        uz0, ids0, m = series0.get_particle( ['uz','id','mass'],
+                                             iteration=series0.iterations[0],
+                                             species=species[s],
+                                             plot=False )
+
+        indeces = np.where(np.in1d(ids0,ids))[0]
+
+        E = np.append(E, 0.5*m[indeces]*(uz0[indeces]*c)**2/MeV)
+        X = np.append(X, x)
+        Y = np.append(Y, y)
+    print(np.min(E), np.max(E))
+    # sort particles according to energy for nicer plot
+    sorted_indeces = np.argsort(E)
+    im = ax.scatter(X[sorted_indeces], Y[sorted_indeces], c=E[sorted_indeces], vmin=vmin, vmax=vmax, cmap=cmap[s])
+
+# dummy plot just to have a neutral colorbar
+im = ax.scatter(np.nan, np.nan, c=np.nan, cmap='Greys_r', vmin=vmin, vmax=vmax)
+plt.colorbar(im, label='E [MeV]')
+ax.set_xlabel('x [m]')
+ax.set_ylabel('y [m]')
+
+plt.tight_layout()
+fig.savefig(f'detect.png', dpi=300)
+plt.close()