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

@@ -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
+        return aperture
\ No newline at end of file
+    
+    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