Add Spin Polarization to Lifetime Calculation.

pyaccel/lifetime.py

@@ -20,7 +20,7 @@
 
 
 class Lifetime:
-    """Class which calculates the lifetime for a given accelerator."""
+    """Class to calculate lifetime contributions for the accelerator."""
 
     # Constant factors
     _MBAR_2_PASCAL = 1.0e-3 / _u.pascal_2_bar
@@ -387,7 +387,8 @@ def touschek_data(self):
 
         If touschek_model = 'FlatBeam', the calculation follows the formulas
             presented in Ref. [1], where the vertical betatron beam size and
-            vertical dispersion are not taken into account
+            vertical dispersion are not taken into account. Dependence of spin
+            polarization is also calculated following Ref. [3].
         If touschek_model = 'Piwinski', the calculation follows the formulas
             presented in Ref. [2], Eqs. 32-42. This formalism describes the
             most general case with respect to the horizontal and vertical
@@ -399,6 +400,11 @@ def touschek_data(self):
                 CERN Acccelerator School: Accelerator Physics (pp. 114–130).
             [2] Piwinski, A. (1999). The Touschek Effect in Strong Focusing
                 Storage Rings. November. http://arxiv.org/abs/physics/9903034
+            [3] J. Zhang, et al. (2013). Precise beam energy measurement using
+                resonant spin depolarization in the SOLEIL storage ring.
+                Nuclear Instruments and Methods in Physics Research Section A:
+                Accelerators, Spectrometers, Detectors and Associated
+                Equipment. https://doi.org/10.1016/j.nima.2012.09.003.
 
         parameters used in calculation:
 
@@ -484,6 +490,7 @@ def touschek_data(self):
         const = (_cst.electron_radius**2 * _cst.light_speed) / (8*_np.pi)
 
         touschek_coeffs = dict()
+        A_spin_num, A_spin_den = None, None
         if self.touschek_model == self.TOUSCHEKMODEL.FlatBeam:
             fator = beta1*eta1l + alpha1*eta1
             a_var = 1 / (4*espread**2) + (eta1**2 + fator**2) / (4*sig1b2)
@@ -512,6 +519,36 @@ def touschek_data(self):
             ratep = const * nr_part/gamma**2 / d_accp**3 * d_pos / vol
             raten = const * nr_part/gamma**2 / d_accn**3 * d_neg / vol
             rate = (ratep+raten)/2
+
+            spin_const = _cst.electron_radius**2
+            spin_const *= _cst.light_speed
+            spin_const /= 2*_np.pi
+            gamma1 = (1 + alpha1**2)/beta1
+            sig1l = _np.sqrt(eta1l**2*espread**2 + gamma1*emit1)
+            a_spinp = spin_const / gamma**3 / vol / sig1l / d_accp**2
+            a_spinn = spin_const / gamma**3 / vol / sig1l / d_accn**2
+
+            ksip = (d_accp/gamma/sig1l)**2
+            ksin = (d_accn/gamma/sig1l)**2
+
+            c_ksip, f_ksip = [], []
+            c_ksin, f_ksin = [], []
+            for en, ep in zip(ksin, ksip):
+                ckp, fkp = self._calc_spin_funcs(ep)
+                ckn, fkn = self._calc_spin_funcs(en)
+                c_ksip.append(ckp)
+                f_ksip.append(fkp)
+                c_ksin.append(ckn)
+                f_ksin.append(fkn)
+
+            c_ksip, f_ksip = _np.array(c_ksip), _np.array(f_ksip)
+            c_ksin, f_ksin = _np.array(c_ksin), _np.array(f_ksin)
+
+            A_spinp_num, A_spinp_den = a_spinp*f_ksip, a_spinp*c_ksip
+            A_spinn_num, A_spinn_den = a_spinn*f_ksin, a_spinn*c_ksin
+            A_spin_num = (A_spinp_num + A_spinn_num)/2
+            A_spin_den = (A_spinp_den + A_spinn_den)/2
+
         elif self.touschek_model == self.TOUSCHEKMODEL.Piwinski:
             eta1til2 = (alpha1*eta1 + beta1*eta1l)**2
             eta2til2 = (alpha2*eta2 + beta2*eta2l)**2
@@ -553,13 +590,20 @@ def touschek_data(self):
             rate = (ratep + raten)/2
 
         rate = _np.array(rate).ravel()
+        circum = s_calc[-1] - s_calc[0]
         # Average inverse Touschek Lifetime
-        avg_rate = _np.trapz(rate, x=s_calc) / (s_calc[-1] - s_calc[0])
-        dit = dict(
+        avg_rate = _np.trapz(rate, x=s_calc)/circum
+        touschek_dict = dict(
             rate=rate, avg_rate=avg_rate,
             volume=vol, pos=s_calc,
             touschek_coeffs=touschek_coeffs)
-        return dit
+        if A_spin_num is not None:
+            A_spin_num = _np.array(A_spin_num).ravel()
+            avg_A_spin_num = _np.trapz(A_spin_num, x=s_calc)
+            avg_A_spin_den = _np.trapz(A_spin_den, x=s_calc)
+            avg_A_spin = avg_A_spin_num/avg_A_spin_den
+            touschek_dict['avg_A_spin'] = avg_A_spin
+        return touschek_dict
 
     @staticmethod
     def f_function_arg(kappa, kappam, b1_, b2_):
@@ -895,6 +939,21 @@ def get_touschek_integration_table(cls, ksi_ini=None, ksi_end=None):
             cls._calc_d_touschek_table(ksi_ini, ksi_end)
         return cls._KSI_TABLE, cls._D_TABLE
 
+    @staticmethod
+    def _calc_spin_funcs(ksi):
+        if _integrate is None:
+            raise ImportError('scipy library not available')
+        lim = 1000
+        int1, _ = _integrate.quad(
+            lambda x: _np.exp(-x)/x, ksi, _np.inf, limit=lim)
+        int2, _ = _integrate.quad(
+            lambda x: _np.exp(-x)*_np.log(x/ksi)/x**2, ksi, _np.inf, limit=lim)
+        int3, _ = _integrate.quad(
+            lambda x: _np.exp(-x)/x**2, ksi, _np.inf, limit=lim)
+        f_val = ksi*int2/2
+        c_val = int1 - f_val - ksi*int3
+        return c_val, f_val
+
     # ----- private methods -----
 
     def _process_eqparams(self, eqparams):

pyaccel/optics/acceptances.py

@@ -129,7 +129,7 @@ def calc_touschek_energy_acceptance(
 
     # ############ Check tunes ############
     # Make sure tunes don't cross int and half-int resonances
-    # Must be symmetric due to syncrhotron oscillations
+    # Must be symmetric due to synchrotron oscillations
     ap_tune = _np.full(energy_offsets.shape, _np.inf)
     if check_tune:
         tune0_int = _np.floor(tune0)

pyaccel/optics/beam_envelope.py

@@ -54,6 +54,7 @@ def __init__(self, accelerator, energy_offset=0.0):
         self._damping_numbers = _np.zeros(3)
         self._tunes = _np.zeros(3)
         self._fixed_point = None
+        self._U0 = None
         self.accelerator = accelerator
 
     def __str__(self):
@@ -306,8 +307,14 @@ def tilt_xyplane(self):
     @property
     def U0(self):
         """."""
-        res = _calc_U0(self._acc.energy, self.energy_offset, self._integral2)
-        return res
+        if self._U0 is None:
+            return _calc_U0(
+                self._acc.energy, self.energy_offset, self._integral2)
+        return self._U0
+
+    @U0.setter
+    def U0(self, value):
+        self._U0 = value
 
     @property
     def overvoltage(self):

pyaccel/optics/rad_integrals.py

@@ -34,6 +34,7 @@ def __init__(self, accelerator, energy_offset=0.0):
         self._alpha = 0.0
         self._integralsx = _np.zeros(6)
         self._integralsy = _np.zeros(6)
+        self._rho3_integral = None
         self.accelerator = accelerator
 
     def __str__(self):
@@ -327,6 +328,20 @@ def sigma_py(self):
         espread0 = self.espread0
         return _np.sqrt(emity*self._twi.gammay + (espread0*self._twi.etapy)**2)
 
+    @property
+    def spin_polarization_time(self):
+        """."""
+        num = 5*_np.sqrt(3)/8
+        Cs = _mp.constants.light_speed
+        compton = _mp.constants.reduced_planck_constant
+        compton /= _mp.constants.electron_mass
+        compton /= _mp.constants.light_speed
+        Cs *= compton
+        Cs *= _mp.constants.electron_radius
+        Cs *= self._acc.gamma_factor**5
+        l0_ = self._acc.length
+        return 1/(num * Cs * self._integralsx[2]/l0_)
+
     def as_dict(self):
         """."""
         dic = {par: getattr(self, par) for par in self.PARAMETERS}