fixed bugs on QNM_Berti function in rdown.py

code_new/Sumit/rdown.py

@@ -35,7 +35,7 @@ q_fpars=[[0.584077, 1.52053, -0.480658], [0.00561441, 0.630715, -0.432664], [-0.
 c=2.99792458*10**8;G=6.67259*10**(-11);MS=1.9885*10**30;
 class Ringdown_Spectrum:
     """RDown model generator"""
-    def __init__(self,mf,af,l,m,n=4,s=-2,time=[],fixed=False,qnm_model='berti'):
+    def __init__(self,mf,af,l,m,n=4,s=-2,time=[],fixed=False,qnm_model='berti',rdowndata=[]):
         self.mf = mf
         self.af = af
         self.l = l
@@ -46,7 +46,8 @@ class Ringdown_Spectrum:
         self.dim = self.n+1
         self.fixed = fixed
         self.qnm_model = qnm_model
-        dict_omega = {'berti': self.QNM_Berti , 'qnm': self.QNM_spectrum}
+        self.dict_omega = {'berti': self.QNM_Berti , 'qnm': self.QNM_spectrum}
+        self.rdowndata = rdowndata
         dic = {'w-tau':self.rd_model_wtau , 'w-q': self.rd_model_wq, 'w-tau-fixed':self.rd_model_wtau_fixed,'w-tau-fixed-m-af': self.rd_model_wtau_m_af}
       
         if len(self.time)==0:
@@ -67,18 +68,18 @@ class Ringdown_Spectrum:
     
         return (w_m_a, tau_m_a)
     
-    def QNM_Berti(self,rdowndata):
+    def QNM_Berti(self,mass,spin):
         """ It computes the RD frequencies and damping times in NR units.
         """     
-        position=np.argmax(rdowndata[0,0] >= (self.af))
+        position=np.argmax(self.rdowndata[0,0] >= (spin))
         #w_m_a=f1+f2*(1-af)**f3
         w_m_a=[None]*(self.n+1)
         tau_ma_a=[None]*(self.n+1)
     
         for i in range(self.n+1):
-            qnm=rdowndata[i,1:3,position]
-            w_m_a[i] = qnm[0]/self.mf
-            tau_ma_a[i] = -1/(qnm[1])*self.mf
+            qnm=self.rdowndata[i,1:3,position]
+            w_m_a[i] = qnm[0]/mass
+            tau_ma_a[i] = -1/(qnm[1])*mass
 
         return w_m_a, tau_ma_a
     
@@ -134,31 +135,31 @@ class Ringdown_Spectrum:
             # -1j to agree with SXS convention
         return ansatz
     
-    def rd_model_wtau_m_af(theta):
+    def rd_model_wtau_m_af(self,theta):
         """RD model parametrized with the damping time tau and with the QNM spectrum fixd to GR. The QNM spectrum is given from the mass and spin.
         """ 
-        xvars = theta[ : (dim)]
-        yvars = theta[(dim) : 2*(dim)]
-        mass_vars = theta[index_mass]
-        spin_vars = theta[index_spin]
+        xvars = theta[ : (self.dim)]
+        yvars = theta[(self.dim) : 2*(self.dim)]
+        mass_vars = theta[-2]
+        spin_vars = theta[-1]
 
-        w_m_a , tau_m_a = dict_omega[self.qnm_model](mass_vars,spin_vars,2,2)
+        w_m_a , tau_m_a = self.dict_omega[self.qnm_model](mass_vars,spin_vars)
 
         ansatz = 0
-        for i in range (0,dim):
-            ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-timesrd_final_tsh/tau_m_a[i]) * (np.cos(w_m_a[i]*timesrd_final_tsh)-1j*np.sin(w_m_a[i]*timesrd_final_tsh))
+        for i in range (0,self.dim):
+            ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-self.time/tau_m_a[i]) * (np.cos(w_m_a[i]*self.time)-1j*np.sin(w_m_a[i]*self.time))
         # -1j to agree with SXS convention
         return ansatz
     
-    def rd_model_wtau_fixed(theta):
+    def rd_model_wtau_fixed(self,theta):
         """RD model parametrized with the damping time tau and with the QNM spectrum fixd to GR. 
         """ 
-        xvars = theta[ : (dim)]
-        yvars = theta[(dim) : 2*(dim)]
+        xvars = theta[ : (self.dim)]
+        yvars = theta[(self.dim) : 2*(self.dim)]
 
         ansatz = 0
-        for i in range (0,dim):
-            ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-timesrd_final_tsh/tau[i]) * (np.cos(w[i]*timesrd_final_tsh)-1j*np.sin(w[i]*timesrd_final_tsh))
+        for i in range (0,self.dim):
+            ansatz += (xvars[i]*np.exp(1j*yvars[i]))*np.exp(-self.time/tau[i]) * (np.cos(w[i]*self.time)-1j*np.sin(w[i]*self.time))
         # -1j to agree with SXS convention
         return ansatz