90 lines
2.7 KiB
Python
90 lines
2.7 KiB
Python
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
|
|
|
|
|
|
|