Finalize thin dipole example.

examples/thin_dipole/analysis_thin_dipole.py

@@ -0,0 +1,104 @@
+#!/usr/bin/env python3
+#
+# Copyright 2022-2023 ImpactX contributors
+# Authors: Axel Huebl, Chad Mitchell
+# License: BSD-3-Clause-LBNL
+#
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import openpmd_api as io
+from scipy.stats import moment
+
+
+def get_moments(beam):
+    """Calculate standard deviations of beam position & momenta
+    and emittance values
+
+    Returns
+    -------
+    sigx, sigy, sigt, emittance_x, emittance_y, emittance_t
+    """
+    sigx = moment(beam["position_x"], moment=2) ** 0.5  # variance -> std dev.
+    sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
+    sigy = moment(beam["position_y"], moment=2) ** 0.5
+    sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
+    sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
+
+    epstrms = beam.cov(ddof=0)
+    emittance_x = (
+        sigx**2 * sigpx**2 - epstrms["position_x"]["momentum_x"] ** 2
+    ) ** 0.5
+    emittance_y = (
+        sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
+    ) ** 0.5
+    emittance_t = (
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
+    ) ** 0.5
+
+    return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
+
+
+# initial/final beam
+series = io.Series("diags/openPMD/monitor.h5", io.Access.read_only)
+last_step = list(series.iterations)[-1]
+initial = series.iterations[1].particles["beam"].to_df()
+final = series.iterations[last_step].particles["beam"].to_df()
+
+# compare number of particles
+num_particles = 10000
+assert num_particles == len(initial)
+assert num_particles == len(final)
+
+print("Initial Beam:")
+sigx, sigy, sigt, emittance_x, emittance_y, emittance_t = get_moments(initial)
+print(f"  sigx={sigx:e} sigy={sigy:e} sigt={sigt:e}")
+print(
+    f"  emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
+)
+
+atol = 0.0  # ignored
+rtol = 1.8 * num_particles**-0.5  # from random sampling of a smooth distribution
+print(f"  rtol={rtol} (ignored: atol~={atol})")
+
+assert np.allclose(
+    [sigx, sigy, sigt, emittance_x, emittance_y, emittance_t],
+    [
+        1.0e-3,
+        1.0e-3,
+        3.0e-1,
+        2.0e-7,
+        2.0e-7,
+        6.0e-5,
+    ],
+    rtol=rtol,
+    atol=atol,
+)
+
+
+print("")
+print("Final Beam:")
+sigx, sigy, sigt, emittance_x, emittance_y, emittance_t = get_moments(final)
+print(f"  sigx={sigx:e} sigy={sigy:e} sigt={sigt:e}")
+print(
+    f"  emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
+)
+
+atol = 0.0  # ignored
+rtol = 1.8 * num_particles**-0.5  # from random sampling of a smooth distribution
+print(f"  rtol={rtol} (ignored: atol~={atol})")
+
+assert np.allclose(
+    [sigx, sigy, sigt, emittance_x, emittance_y, emittance_t],
+    [
+        1.0e-3,
+        1.0e-3,
+        3.0e-1,
+        2.0e-7,
+        2.0e-7,
+        6.0e-5,
+    ],
+    rtol=rtol,
+    atol=atol,
+)

examples/thin_dipole/input_thin_dipole.in

@@ -3,7 +3,7 @@
 ###############################################################################
 beam.npart = 10000
 beam.units = static
-beam.kin_energy = 4.0e-3
+beam.kin_energy = 5.0
 beam.charge = 1.0e-9
 beam.particle = proton
 beam.distribution = waterbag
@@ -12,7 +12,7 @@ beam.sigmaY = 1.0e-3
 beam.sigmaT = 0.3
 beam.sigmaPx = 2.0e-4
 beam.sigmaPy = 2.0e-4
-beam.sigmaPt = 2.0e-5
+beam.sigmaPt = 2.0e-4
 beam.muxpx = 0.0
 beam.muypy = 0.0
 beam.mutpt = 0.0
@@ -21,28 +21,28 @@ beam.mutpt = 0.0
 ###############################################################################
 # Beamline: lattice elements and segments
 ###############################################################################
-lattice.elements = monitor bend monitor
+lattice.elements = monitor bend inverse_bend monitor
 lattice.nslice = 1
 
 monitor.type = beam_monitor
 monitor.backend = h5
 
-#a 360 deg sbend using drift-kick-drift:
+#90 degree sbend using drift-kick-drift:
 bend.type = line
 bend.elements = dr kick dr
-bend.repeat = 100
+bend.repeat = 200
 
 dr.type = drift
-dr.ds = 0.0314159265358979
+dr.ds = 0.003926990816987
 
 kick.type = thin_dipole
-kick.theta = 3.6
+kick.theta = 0.45
 kick.rc = 1.0
 
-#a 360 sbend using exact nonlinear map:
-#bend.type = sbend_exact
-#bend.ds = 6.283185307179586
-#bend.phi = 360.0
+#inverse of a 90 degree sbend using exact nonlinear map:
+inverse_bend.type = sbend_exact
+inverse_bend.ds = -1.570796326794897
+inverse_bend.phi = -90.0
 
 ###############################################################################
 # Algorithms

src/particles/elements/ThinDipole.H

@@ -83,7 +83,7 @@ namespace impactx
 
             // compute the function expressing dp/p in terms of pt (labeled f in Ripken etc.)
             amrex::ParticleReal f = -1.0_prt + sqrt(1.0_prt - 2.0_prt*pt/beta_ref + pow(pt,2));
-            amrex::ParticleReal fprime = (1.0_prt + beta_ref*pt)/(beta_ref*(1.0_prt + f));
+            amrex::ParticleReal fprime = (1.0_prt - beta_ref*pt)/(beta_ref*(1.0_prt + f));
 
             // compute the effective (equivalent) arc length and curvature
             amrex::ParticleReal ds = m_theta*m_rc;
@@ -96,7 +96,7 @@ namespace impactx
             p.pos(RealAoS::y) = y;
             pyout = py;
 
-            p.pos(RealAoS::t) = t - kx*x*ds*fprime; //eq (3.2e)
+            p.pos(RealAoS::t) = t + kx*x*ds*fprime; //eq (3.2e)
             ptout = pt;
 
             // assign updated momenta