Feature/electron cooler

ducktrack/elements.py

@@ -9,6 +9,7 @@
 from scipy.constants import mu_0
 from scipy.constants import m_e as me_kg
 from scipy.constants import e as qe
+from scipy.constants import physical_constants
 
 me = me_kg*clight**2/qe
 
@@ -840,6 +841,137 @@ def track(self,p):
         if self.length > 0.0:
             raise NotImplementedError('Radiation is not implemented')
 
+class ElectronCooler(Element):
+
+    _description =[
+        ("current","","",0.0),
+        ("length","","",0.0),
+        ("radius_e_beam","","",0.0),
+        ("temp_perp","","",0.0),
+        ("temp_long","","",0.0),
+        ("magnetic_field","","",0.0),
+
+        ("offset_x","","",0.0),
+        ("offset_px","","",0.0),
+        ("offset_y","","",0.0),
+        ("offset_py","","",0.0),
+        ("offset_energy","","",0.0),
+        ("magnetic_field_ratio","","",0.0),
+        ("space_charge_factor","","",0.0),
+
+        ("classical_e_radius","","",physical_constants['classical electron radius'][0]),
+        ("me_kg","","",me_kg),
+        ]
+
+    def force(self, p):
+        current=self.current
+        length=self.length
+        radius_e_beam=self.radius_e_beam
+        temp_perp=self.temp_perp
+        temp_long=self.temp_long
+        magnetic_field=self.magnetic_field
+        magnetic_field_ratio=self.magnetic_field_ratio
+        space_charge_factor=self.space_charge_factor
+        classical_e_radius=self.classical_e_radius
+        me_kg = self.me_kg    
+
+        # All parameters are taken relative to the electron beam
+        x     = p.x     - self.offset_x
+        px    = p.px    - self.offset_px
+        y     = p.y     - self.offset_y
+        py    = p.py    - self.offset_py
+        delta = p.delta # offset_energy is implemented when longitudinal velocity is computed    
+        x, px, y, py, delta = map(np.atleast_1d, (x, px, y, py, delta))
+
+        # Radial and angular coordinates
+        theta = np.arctan2(y, x)
+        radius = np.hypot(x,y)
+
+        # Particle beam parameters
+        beta0, gamma0, q0, p0c, mass0 = p.beta0, p.gamma0, p.q0, p.p0c, p.mass0
+        total_momentum = p0c * (1.0 + delta)
+        gamma = np.sqrt(1.0 + (total_momentum / mass0) ** 2)
+        beta = np.sqrt(1.0 - 1.0 / gamma**2)
+        beta_x = px * p0c / (mass0 * gamma)
+        beta_y = py * p0c / (mass0 * gamma)
+
+        # Compute electron density
+        volume_e_beam = np.pi*(radius_e_beam)**2*length #m3
+        num_e_per_s = current/qe # number of electrons per second
+        self.tau=length/(gamma0*beta0*clight) # time spent in the electron cooler
+        electron_density = num_e_per_s*self.tau/volume_e_beam # density of electrons     
+
+        # Electron beam properties
+        v_perp_temp = (qe*temp_perp/me_kg)**(1./2) # transverse electron rms velocity
+        v_long_temp = (qe*temp_long/me_kg)**(1./2) # longitudinal electron rms velocity
+        rho_larmor = me_kg * v_perp_temp / (qe * magnetic_field) # depends on transverse temperature, larmor radius
+        elec_plasma_frequency = np.sqrt(electron_density * qe**2 / (me_kg * epsilon_0))
+        v_rms_magnet = beta0 * gamma0 * clight * magnetic_field_ratio # velocity spread due to magnetic imperfections
+        #V_eff = np.sqrt(v_long_temp**2 + v_rms_magnet**2) # effective electron beam velocity spread
+        mass_electron_ev = me_kg * clight**2 / qe #eV
+        energy_electron_initial = (gamma0 - 1) * mass_electron_ev #eV
+        energy_e_total = energy_electron_initial + self.offset_energy
+
+        friction_coefficient = electron_density*q0**2*qe**4 /(4*me_kg*(np.pi*epsilon_0)**2) # Coefficient used for computation of friction force 
+
+        # compute angular frequency of rotation of e-beam due to space charge
+        omega_e_beam = space_charge_factor*1/(2*np.pi*epsilon_0*clight) * current/(radius_e_beam**2*beta0*gamma0*magnetic_field)
+
+        Fx = np.zeros_like(x)
+        Fy = np.zeros_like(y)
+        Fl = np.zeros_like(delta)        
+
+        # Radial_velocity_dependence due to space charge
+        #  -> from equation 100b in Helmut Poth: Electron cooling. page 186      
+        space_charge_coefficient = space_charge_factor *classical_e_radius / (qe * clight) * (gamma0 + 1) / (gamma0**2);# //used for computation of the space charge energy offset
+        dE_E = space_charge_coefficient * current * (radius / radius_e_beam)**2 / (beta0)**3
+        E_diff_space_charge = dE_E * energy_e_total        
+
+        E_kin_total = energy_e_total + E_diff_space_charge
+        gamma_total = 1 + (E_kin_total / mass_electron_ev)
+        beta_total  = np.sqrt(1 - 1 / (gamma_total**2))
+
+        # Velocity differences
+        dVz = beta * clight - beta_total* clight      # Longitudinal velocity difference              
+        dVx = beta_x * clight                         # Horizontal velocity difference            
+        dVy = beta_y * clight                         # Vertical velocity difference
+        dVx -= omega_e_beam * radius * -np.sin(theta) # Apply x-component of e-beam rotation
+        dVy -= omega_e_beam * radius * +np.cos(theta) # Apply y-component of e-beam rotation
+        dV_squared = dVx**2+dVy**2+dVz**2
+        V_tot = np.sqrt(dV_squared + v_long_temp**2 + v_rms_magnet**2) # Total velocity difference due to all effects
+
+        # Coulomb logarithm        
+        rho_min = q0 *qe**2/(4*np.pi*epsilon_0*me_kg*V_tot**2)       # Minimum impact parameter
+        rho_max_shielding = V_tot/(elec_plasma_frequency)            # Maximum impact parameter based on charge shielding
+        rho_max_interaction = V_tot*self.tau                         # Maximum impact parameter based on interaction time
+        rho_max = np.minimum(rho_max_shielding, rho_max_interaction) # Take the smaller of the two maximum impact parameters
+        log_coulomb = np.log((rho_max+rho_min+rho_larmor)/(rho_min+rho_larmor)) # Coulomb logarithm
+
+        friction_denominator = V_tot**3 # Compute this coefficient once because its going to be used three times
+
+        Fx = -friction_coefficient * dVx/friction_denominator * log_coulomb  # Newton
+        Fy = -friction_coefficient * dVy/friction_denominator * log_coulomb  # Newton
+        Fl = -friction_coefficient * dVz/friction_denominator * log_coulomb  # Newton
+        # If particle is outside electron beam, set cooling force to zero
+        outside_beam_indices = radius >= radius_e_beam
+        Fx[outside_beam_indices] = 0.0
+        Fy[outside_beam_indices] = 0.0
+        Fl[outside_beam_indices] = 0.0
+
+        Fx = Fx * 1/qe # convert to eV/m 
+        Fy = Fy * 1/qe # convert to eV/m 
+        Fl = Fl * 1/qe # convert to eV/m 
+        return Fx,Fy,Fl
+
+    def track(self, p):
+        Fx,Fy,Fl=self.force(p)
+        Fx = Fx * clight # convert to eV/c because p0c is also in eV/c
+        Fy = Fy * clight # convert to eV/c because p0c is also in eV/c
+        Fl = Fl * clight # convert to eV/c because p0c is also in eV/c
+        p.delta += np.squeeze( Fl*p.gamma0*self.tau/p.p0c)
+        p.px    += np.squeeze( Fx*p.gamma0*self.tau/p.p0c)
+        p.py    += np.squeeze( Fy*p.gamma0*self.tau/p.p0c)        
+
 __all__ = [
     "BeamBeam4D",
     "BeamBeam6D",
@@ -860,5 +992,6 @@ def track(self,p):
     "SRotation",
     "XYShift",
     "LinearTransferMatrix",
-    "FirstOrderTaylorMap"
+    "FirstOrderTaylorMap",
+    "ElectronCooler"
 ]

examples/electron_cooler/Electron_cooler_IPAC2023_results.py

@@ -0,0 +1,148 @@
+# copyright ############################### #
+# This file is part of the Xtrack Package.  #
+# Copyright (c) CERN, 2021.                 #
+# ######################################### #
+
+# This file will generate and display the following IPAC results:
+# https://accelconf.web.cern.ch/ipac2023/doi_per_institute/tupm027/index.html
+# Specifications for a new electron cooler of the antiproton decelerator at CERN 
+
+
+import numpy as np
+import xtrack as xt
+import xpart as xp
+import matplotlib.pyplot as plt
+
+######################################
+# lattice parameters for AD at 300 MeV
+######################################
+
+#lattice parameters for AD at 300Mev/c
+# from https://acc-models.web.cern.ch/acc-models/ad/scenarios/lowenergy/lowenergy.tfs
+qx = 5.45020077392
+qy = 5.41919929346
+dqx=-20.10016919292
+dqy=-22.29552755573
+circumference = 182.43280000000 #m
+
+# relativistic factors
+gamma_rel = 1.04987215550 # at 300 MeV/c
+
+# optics at e-cooler (approximate), in m
+beta_x = 10 
+beta_y = 4
+D_x = 0
+
+# electron cooler parameters
+current = 2.4  # A current
+length = 1.5 # m cooler length
+radius_e_beam = 25*1e-3 #m radius of the electron beam
+temp_perp = 100e-3 # <E> [eV] = kb*T
+temp_long =  1e-3 # <E> [eV]
+magnetic_field = 0.060 # 100 Gauss in ELENA
+# idea is to study magnetic field imperfections
+magnetic_field_ratio_list = [0,1e-4,5e-4,1e-3] #Iterate over different values of the magnetic field quality to see effect on cooling performance.
+#magnetic_field_ratio is the ratio of transverse componenet of magnetic field and the longitudinal component. In the ideal case, the ratio is 0.
+
+# some initial beam parameters
+emittance = 35e-6
+dp_p = 2e-3 
+q0 = 1
+
+# simulation parameters: simulate 20 s of cooling, and take data once every 100 ms
+max_time_s = 20
+int_time_s = 0.01
+
+# some constants, and simple computations
+clight = 299792458.0
+mass0 = 938.27208816*1e6 #ev/c^2
+
+beta_rel = np.sqrt(gamma_rel**2 - 1)/gamma_rel
+p0c = mass0*beta_rel*gamma_rel #eV/c
+T_per_turn = circumference/(clight*beta_rel)
+
+# compute length of simulation, as well as sample interval, in turns
+num_turns = int(max_time_s/T_per_turn)
+save_interval = int(int_time_s/T_per_turn)
+
+# compute initial beam parameters
+x_init = np.sqrt(beta_x*emittance)
+y_init = np.sqrt(beta_y*emittance)
+
+particles0 = xp.Particles(
+            mass0=mass0,
+            p0c=p0c,
+            q0=q0,
+            x=x_init,
+            px=0,
+            y=0,
+            py=0,
+            delta=0,
+            zeta=0)        
+
+# arc to do linear tracking
+arc = xt.LineSegmentMap(
+        qx=qx, qy=qx,
+        dqx=dqx, dqy=dqy,
+        length=circumference,
+        betx=beta_x,
+        bety=beta_y,
+        dx=D_x)
+
+# compute length of simulation, as well as sample interval, in turns
+num_turns = int((max_time_s / T_per_turn).item())
+save_interval = int((int_time_s / T_per_turn).item())
+
+# create a monitor object, to reduce holded data
+monitor = xt.ParticlesMonitor(start_at_turn=0, stop_at_turn=1,
+                        n_repetitions=int(num_turns/save_interval),
+                        repetition_period=save_interval,
+                        num_particles=1)
+
+##############################
+# electron cooler parameters #
+##############################  
+electron_cooler = xt.ElectronCooler(
+                length=length,
+                radius_e_beam=radius_e_beam,
+                current=current,
+                temp_perp=temp_perp,
+                temp_long=temp_long,
+                magnetic_field=magnetic_field, 
+                magnetic_field_ratio=0,
+                space_charge_factor=0)
+
+line = xt.Line(elements=[monitor, electron_cooler, arc],element_names=['monitor','electron_cooler','arc'])
+line.particle_ref = xp.Particles(mass0=mass0, q0=q0, p0c=p0c)
+line.build_tracker()
+
+################
+# prepare plot #
+################
+
+plt.figure(figsize=(10, 5))
+plt.rcParams.update({'font.size': 14}) 
+
+##########################################
+# scan over magnetic field imperfections #
+##########################################
+
+for magnetic_field_ratio in magnetic_field_ratio_list:        
+
+        electron_cooler.magnetic_field_ratio=magnetic_field_ratio
+        particles=particles0.copy()
+        # track
+        line.track(particles, num_turns=num_turns,
+              turn_by_turn_monitor=False)
+        x = monitor.x[:,:,0]
+        px = monitor.px[:,:,0]        
+        time = monitor.at_turn[:, 0, 0] * T_per_turn
+        # compute action at the end for all turns
+        action_x = (x**2/beta_x + beta_x*px**2)
+
+        plt.plot(time,action_x*1e6,label='$B_{\\perp}/B_{\\parallel}$='f'{magnetic_field_ratio:.0e}')
+
+plt.xlabel('Time [s]')
+plt.ylabel('$J_x$ $[\\mu m]$')
+plt.legend()
+plt.show()