diff --git a/tracking/optics.py b/tracking/optics.py
index 089b271adddd11dcfa063f9c5d4d93e19d9f25dd..2482a373478ba532fb17b22db7cf4286324d0606 100644
--- a/tracking/optics.py
+++ b/tracking/optics.py
@@ -122,4 +122,71 @@ class Optics:
                           self.dispY(position), self.disppY(position)]
             return np.array(dispersion)
     
+class PhyisicalModel:
+    """Store lattice physical parameters"""
+    def __init__(self, ring, a0=5e-3, b0=5e-3, shape='circ', rho=1.7e-8,
+                 n_points=1e4):
+        
+        self.n_points = int(n_points)
+        self.position = np.linspace(0, ring.L, self.n_points)
+        self.x0 = np.ones_like(self.position)*a0 # Horizontal half aperture [m]
+        self.y0 = np.ones_like(self.position)*b0 # Vertical half aperture [m]
+        
+        self.length = self.position[1:] - self.position[:-1]
+        self.center = (self.position[1:] + self.position[:-1])/2
+        self.rho = np.ones_like(self.center)*rho # Resistivity of the chamber material [ohm*m]
+        
+        self.shape = np.repeat(np.array([shape]), self.n_points-1) # Shape of the chamber cross section (elli/circ/rect)
 
+    def resistive_wall_coef(self):
+        L = self.end_position[-1]
+        sigma = 1/self.rho
+        beta_H_star = 1/L*(self.length*self.beta[0,:]).sum()
+        beta_V_star = 1/L*(self.length*self.beta[1,:]).sum()
+        sigma_star = 1/L*(self.length*sigma).sum()
+        
+        a1_H = (((self.length*self.beta[0,:]/(np.sqrt(sigma)*(self.a0)**1)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/1)
+        a2_H = (((self.length*self.beta[0,:]/(np.sqrt(sigma)*(self.a0)**2)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/2)
+
+        a1_V = (((self.length*self.beta[1,:]/(np.sqrt(sigma)*(self.b0)**1)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/1)
+        a2_V = (((self.length*self.beta[1,:]/(np.sqrt(sigma)*(self.b0)**2)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/2)
+        
+        a3_H = (((self.length*self.beta[0,:]/(np.sqrt(sigma)*(self.a0)**3)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/3)
+        a4_H = (((self.length*self.beta[0,:]/(np.sqrt(sigma)*(self.a0)**4)).sum())**(-1)*L*beta_H_star/np.sqrt(sigma_star))**(1/4)
+        
+        a3_V = (((self.length*self.beta[1,:]/(np.sqrt(sigma)*(self.b0)**3)).sum())**(-1)*L*beta_V_star/np.sqrt(sigma_star))**(1/3)
+        a4_V = (((self.length*self.beta[1,:]/(np.sqrt(sigma)*(self.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))
+        
+        print("rhorw = " + str(1/sigma_star))
+        print("beffL1 = " + str(a1_L))
+        print("beffL2 = " + str(a2_L))
+        print("beffH3 = " + str(a3_H))
+        print("beffH4 = " + str(a4_H))
+        print("beffV3 = " + str(a3_V))
+        print("beffV4 = " + str(a4_V))
+        
+    def change_values(self, pos1, pos2, a0, b0, shape, rho):
+        ind = (self.position > pos1) & (self.position < pos2)
+        self.x0[ind] = a0
+        self.y0[ind] = b0
+        ind2 = (self.position[:-1] > pos1) & (self.position[1:] < pos2)
+        self.rho[ind2] = rho
+        self.shape[ind2] = shape
+        
+    def taper(self, pos1, pos2, x0_start, y0_start, x0_end, y0_end, shape, rho):
+        ind = (self.position > pos1) & (self.position < pos2)
+        self.x0[ind] = np.linspace(x0_start, x0_end, sum(ind))
+        self.y0[ind] = np.linspace(y0_start, y0_end, sum(ind))
+        ind2 = (self.position[:-1] > pos1) & (self.position[1:] < pos2)
+        self.rho[ind2] = rho
+        self.shape[ind2] = shape
+        
+    def plot_aperture(self):
+        plt.plot(self.position,self.x0*1e3)
+        plt.plot(self.position,self.y0*1e3)
+        plt.xlabel("Longitudinal position [m]")
+        plt.ylabel("H/V aperture [mm]")
+        plt.legend(["H","V"])