From 8e4a8ff26e1780ef336c01043dfe8f3ab194dfbf Mon Sep 17 00:00:00 2001
From: Gamelin Alexis <alexis.gamelin@synchrotron-soleil.fr>
Date: Sat, 22 Apr 2023 15:09:50 +0200
Subject: [PATCH] Improve PhysicalModel class
Sections with rho=0 are not considered in resistive_wall_effective_radius any more.
Add a plot_resistivity method.
---
mbtrack2/utilities/optics.py | 54 ++++++++++++++++++++++++------------
1 file changed, 37 insertions(+), 17 deletions(-)
diff --git a/mbtrack2/utilities/optics.py b/mbtrack2/utilities/optics.py
index bcb3ed0..979ed94 100644
--- a/mbtrack2/utilities/optics.py
+++ b/mbtrack2/utilities/optics.py
@@ -481,17 +481,21 @@ class PhysicalModel:
"""
Return the effective radius of the chamber for resistive wall
calculations as defined in Eq. 27 of [1].
+
+ Section with rho = 0 ohm.m are not considered in this calculation.
Parameters
----------
optics : Optics object
x_right : bool, optional
If True, x_right is used, if Fasle, x_left is used.
- y_top : TYPE, optional
+ y_top : bool, optional
If True, y_top is used, if Fasle, y_bottom is used.
Returns
-------
+ L : float
+ Total length considered in [m].
rho_star : float
Effective resistivity of the chamber material in [ohm.m].
a1_L : float
@@ -520,6 +524,8 @@ class PhysicalModel:
Spectrometers, Detectors and Associated Equipment 806 (2016): 221-230.
"""
+ idx = self.rho != 0
+
if x_right is True:
a0 = (self.x_right[1:] + self.x_right[:-1])/2
else:
@@ -529,30 +535,35 @@ class PhysicalModel:
b0 = (self.y_top[1:] + self.y_top[:-1])/2
else:
b0 = np.abs((self.y_bottom[1:] + self.y_bottom[:-1])/2)
+
+ a0 = a0[idx]
+ b0 = b0[idx]
- beta = optics.beta(self.center)
- L = self.position[-1]
- sigma = 1/self.rho
- beta_H_star = 1/L*(self.length*beta[0,:]).sum()
- beta_V_star = 1/L*(self.length*beta[1,:]).sum()
- sigma_star = 1/L*(self.length*sigma).sum()
+ beta = optics.beta(self.center[idx])
+ length = self.length[idx]
- a1_H = (((self.length*beta[0,:]/(np.sqrt(sigma)*(a0)**1)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/1)
- a2_H = (((self.length*beta[0,:]/(np.sqrt(sigma)*(a0)**2)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/2)
+ L = length.sum()
+ sigma = 1/self.rho[idx]
+ beta_H_star = 1/L*(length*beta[0,:]).sum()
+ beta_V_star = 1/L*(length*beta[1,:]).sum()
+ sigma_star = 1/L*(length*sigma).sum()
+
+ a1_H = (((length*beta[0,:]/(np.sqrt(sigma)*(a0)**1)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/1)
+ a2_H = (((length*beta[0,:]/(np.sqrt(sigma)*(a0)**2)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/2)
- a1_V = (((self.length*beta[1,:]/(np.sqrt(sigma)*(b0)**1)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/1)
- a2_V = (((self.length*beta[1,:]/(np.sqrt(sigma)*(b0)**2)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/2)
+ a1_V = (((length*beta[1,:]/(np.sqrt(sigma)*(b0)**1)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/1)
+ a2_V = (((length*beta[1,:]/(np.sqrt(sigma)*(b0)**2)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/2)
- a3_H = (((self.length*beta[0,:]/(np.sqrt(sigma)*(a0)**3)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/3)
- a4_H = (((self.length*beta[0,:]/(np.sqrt(sigma)*(a0)**4)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/4)
+ a3_H = (((length*beta[0,:]/(np.sqrt(sigma)*(a0)**3)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/3)
+ a4_H = (((length*beta[0,:]/(np.sqrt(sigma)*(a0)**4)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/4)
- a3_V = (((self.length*beta[1,:]/(np.sqrt(sigma)*(b0)**3)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/3)
- a4_V = (((self.length*beta[1,:]/(np.sqrt(sigma)*(b0)**4)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/4)
+ a3_V = (((length*beta[1,:]/(np.sqrt(sigma)*(b0)**3)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/3)
+ a4_V = (((length*beta[1,:]/(np.sqrt(sigma)*(b0)**4)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/4)
a1_L = min((a1_H,a1_V))
a2_L = min((a2_H,a2_V))
- return (1/sigma_star, a1_L, a2_L, a3_H, a4_H, a3_V, a4_V)
+ return (L, 1/sigma_star, beta_H_star, beta_V_star, a1_L, a2_L, a3_H, a4_H, a3_V, a4_V)
def change_values(self, start_position, end_position, x_right=None,
y_top=None, shape=None, rho=None, x_left=None,
@@ -684,4 +695,13 @@ class PhysicalModel:
self.xm(s),
self.yp(s),
self.ym(s)])
- return aperture
\ No newline at end of file
+ return aperture
+
+ def plot_resistivity(self):
+ """Plot resistivity along the ring."""
+ fig, ax = plt.subplots(1)
+ ax.plot(self.position[1:], self.rho)
+ ax.set(xlabel="Longitudinal position [m]",
+ ylabel="Resistivity [ohm.m]")
+
+ return (fig, ax)
\ No newline at end of file
--
GitLab