added DelayCurrentConverter

devices/delay_current.py

@@ -0,0 +1,89 @@
+import numpy as np
+import pint
+from scipy import optimize
+
+
+unit = pint.UnitRegistry()
+MOMENTUM_SCALE = (1 * unit.GeV / unit.c).to_base_units().magnitude
+
+N_MAGNETS = 4
+N_DRIFTS = 2
+
+MAGNET_LENGTH = 0.3109 * unit.m
+DRIFT_LENGTH = 0.42 * unit.m
+DISTANCE_BETWEEN_MAGNETS_2_AND_3 = 0.16 * unit.m
+
+TOTAL_LENGTH = N_MAGNETS * MAGNET_LENGTH + N_DRIFTS * DRIFT_LENGTH + DISTANCE_BETWEEN_MAGNETS_2_AND_3
+
+# unclear what these magic numbers are:
+MAGNUM1 = 0.2998
+MAGNUM2 = 0.001934
+
+
+class DelayCurrentConverter:
+
+    def __init__(self, I_guess=30):
+        """
+        I_guess in ampere, must be dimensionless for minimizer
+        """
+        self.I_guess = I_guess
+        self.vectorized_current_A = np.vectorize(self.current_A)
+        self.vectorized_delay_fs = np.vectorize(self.delay_fs)
+
+
+    # dimensionless functions
+
+    def current_A(self, delay_fs, beam_energy_MeV):
+        delay = delay_fs * unit.fs
+        beam_momentum = beam_energy_MeV * unit.MeV / unit.c
+        return self.current(delay, beam_momentum).to("A").magnitude
+
+    def delay_fs(self, I_ampere, beam_energy_MeV):
+        I = I_ampere * unit.A
+        beam_momentum = beam_energy_MeV * unit.MeV / unit.c
+        return self.delay(I, beam_momentum).to("fs").magnitude
+
+
+    # functions below expect pint
+
+    def current(self, delay, beam_momentum):
+        res = optimize.root_scalar(self.delay_error, args=(delay, beam_momentum), x0=self.I_guess, x1=2*self.I_guess)
+        return res.root * unit.A
+
+    def delay(self, I, beam_momentum):
+        delta = self.path_length_difference(I, beam_momentum) / unit.c
+        return delta.to("fs")
+
+
+    def delay_error(self, I_ampere, delay_target, beam_momentum):
+        I = I_ampere * unit.A
+        result_delay = self.delay(I, beam_momentum)
+        distance = delay_target - result_delay
+        return distance.to("fs").magnitude
+
+    def path_length_difference(self, I, beam_momentum):
+        length = self.path_length(I, beam_momentum)
+        return length - TOTAL_LENGTH
+
+    def path_length(self, I, beam_momentum):
+        angle = self.deflection_angle(I, beam_momentum)
+        magnet = N_MAGNETS * self.orbit_radius_magnet(I, beam_momentum) * angle
+        drift = N_DRIFTS * DRIFT_LENGTH / np.cos(angle)
+        total = magnet + drift + DISTANCE_BETWEEN_MAGNETS_2_AND_3
+        return total
+
+    def orbit_radius_magnet(self, I, beam_momentum):
+        angle = self.deflection_angle(I, beam_momentum)
+        return MAGNET_LENGTH / np.sin(angle)
+
+    def deflection_angle(self, I, beam_momentum):
+        field = self.integrate_magnet_field(I)
+        scaled = beam_momentum / MOMENTUM_SCALE
+        argument = field / scaled * MAGNUM1 * unit.C
+        return np.arcsin(argument)
+
+    def integrate_magnet_field(self, I):
+        return I * MAGNUM2 * unit.m * unit.T / unit.A
+
+
+