updating

code_new/Sumit/rdown.py

@@ -28,13 +28,14 @@ import numpy as np
 import qnm
 import os
 
-
 f_fpars= [[2.95845, -2.58697, 0.0533469], [2.12539, -1.78054, 0.0865503], [1.74755, -1.44776, 0.123666], [1.78287, -1.53203, 0.129475], [2.04028, -1.83224, 0.112497]]
 q_fpars=[[0.584077, 1.52053, -0.480658], [0.00561441, 0.630715, -0.432664], [-0.197965, 0.515956, -0.369706], [-0.275097, 0.455691, -0.331543], [-0.287596, 0.398514, -0.309799]]
 
+
+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):
+    def __init__(self,mf,af,l,m,n=4,s=-2,time=[],fixed=False,qnm_model='berti'):
         self.mf = mf
         self.af = af
         self.l = l
@@ -44,7 +45,10 @@ class Ringdown_Spectrum:
         self.grav_220 = [qnm.modes_cache(s=s,l=self.l,m=self.m,n=i) for i in range (0,self.n+1)]
         self.dim = self.n+1
         self.fixed = fixed
-        
+        self.qnm_model = qnm_model
+        dict_omega = {'berti': self.QNM_Berti , 'qnm': self.QNM_spectrum}
+        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:
             self.time = np.arange(0,100,0.1)
             
@@ -62,7 +66,22 @@ class Ringdown_Spectrum:
         tau_m_a=-1/(np.imag(omegas_new))*self.mf
     
         return (w_m_a, tau_m_a)
+    
+    def QNM_Berti(self,rdowndata):
+        """ It computes the RD frequencies and damping times in NR units.
+        """     
+        position=np.argmax(rdowndata[0,0] >= (self.af))
+        #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
 
+        return w_m_a, tau_ma_a
+    
 
     def w_fpars_Berti(self,n):
         return f_fpars[n]
@@ -89,7 +108,7 @@ class Ringdown_Spectrum:
             q1,q2,q3 = tau_qpars_Berti(n)
             res[n]=(f1 + f2*(2**(-1/q3)*((-2*q1 + w[n]*tau[n])/q2)**(1/q3))**f3)/w[n]
         return res
-
+    
     def spin_from_wtau_loop(self,w,tau,l,m):
         res=[None]*dim
         for n in range (0,dim):
@@ -115,6 +134,34 @@ class Ringdown_Spectrum:
             # -1j to agree with SXS convention
         return ansatz
     
+    def rd_model_wtau_m_af(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]
+
+        w_m_a , tau_m_a = dict_omega[self.qnm_model](mass_vars,spin_vars,2,2)
+
+        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))
+        # -1j to agree with SXS convention
+        return ansatz
+    
+    def rd_model_wtau_fixed(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)]
+
+        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))
+        # -1j to agree with SXS convention
+        return ansatz
+    
     def rd_model_wq(self,theta):
         """RD model parametrized with the quality factor q.
         """  
@@ -152,14 +199,10 @@ class Ringdown_Spectrum:
         mass_vars = theta[-2]
         spin_vars = theta[-1]
 
-        w_m_a , tau_m_a = QNM_spectrum
+        w_m_a , tau_m_a = QNM_spectrum()
 
         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))
         # -1j to agree with SXS convention
-        return ansatz
-
-
-
-
+        return ansatz
\ No newline at end of file