changed binary.py in a first version to derive the fisher matrix with respect...

changed binary.py in a first version to derive the fisher matrix with respect to the parameters ln(A), position of the source, angular momentum of the source, planetary parameters

binary.py

@@ -101,6 +101,8 @@ class binary:
         self.K = (2*pi*G/self.P)**(1/3)*self.mP/(self.m1+self.m2+self.mP)**(2/3)*np.sin(theta_P)/c
         self.T_obs = T_obs*yr #s
         
+        self.key = key
+        
         self.chirp_mass = (self.m1*self.m2)**(3/5)/((self.m1+self.m2))**(1/5) #kg
         self.f_binary = freq/2
         self.f_GW = freq
@@ -201,7 +203,9 @@ class binary:
     def strain(self,t):
         return sqth*self.A(t)*np.cos(2*pi*self.freq_int(t)+self.polarization_phase(t)+self.doppler_phase(t))
 
-    def d_freq(self,t,i):
+    def d_freq(self,t,i): # f_0, K, P, phi_0
+        if i == 0:
+            return (1-self.K*np.cos(2*pi*t/self.P+self.phi_0))
         if i == 1:
             return -self.f_GW*np.cos(2*pi*t/self.P+self.phi_0)
         if i == 2:
@@ -224,23 +228,56 @@ class binary:
             print('Error in d_freq_int')
             return self.strain(t)
         
-    def d_strain(self,t,i):
-        return -sqth*self.A(t)*(2*pi*self.d_freq_int(t,i)+2*pi*self.d_freq(t,i)/c*AU*np.sin(self.theta_S_Ec)*np.cos(omega_E*t-self.phi_S_Ec))*np.sin(2*pi*self.freq_int(t)+self.polarization_phase(t)+self.doppler_phase(t))
+    def d_strain(self,t,i,bina=None,diff=1.): # 0: lnA, 1: f_0, 2: theta_S_Ec, 3: phi_S_Ec, 4: theta_L, 5: phi_L, 6: K, 7: P, 8: phi_0
+        if bina is not None:
+            return (bina.strain(t)-self.strain(t))/diff
+        elif i >= 6:
+            i -= 5
+            return -sqth*self.A(t)*(2*pi*self.d_freq_int(t,i)+2*pi*self.d_freq(t,i)/c*AU*np.sin(self.theta_S_Ec)*np.cos(omega_E*t-self.phi_S_Ec))*np.sin(2*pi*self.freq_int(t)+self.polarization_phase(t)+self.doppler_phase(t))
+        elif i == 0:
+            return self.strain(t)
+        elif i == 1:
+            return -sqth*self.A(t)*(2*pi*self.T_obs+2*pi*self.d_freq(t,0)/c*AU*np.sin(self.theta_S_Ec)*np.cos(omega_E*t-self.phi_S_Ec))*np.sin(2*pi*self.freq_int(t)+self.polarization_phase(t)+self.doppler_phase(t))
+        else:
+            raise ValueError
+        
+    def h_ixh_j(self,i,j,diff=1e-7):
+        B = None
+        C = None
+        if i == 2:
+            B = binary(self.theta_S_Ec+diff,self.phi_S_Ec,self.dist/pc,self.theta_L,self.phi_L,self.m1/M_S,self.m2/M_S,self.f_GW,self.mP/M_J,self.P/yr,self.theta_P,self.phi_0,self.T_obs/yr,self.key)
+        elif i == 3:
+            B = binary(self.theta_S_Ec,self.phi_S_Ec+diff,self.dist/pc,self.theta_L,self.phi_L,self.m1/M_S,self.m2/M_S,self.f_GW,self.mP/M_J,self.P/yr,self.theta_P,self.phi_0,self.T_obs/yr,self.key)
+        elif i == 4:
+            B = binary(self.theta_S_Ec,self.phi_S_Ec,self.dist/pc,self.theta_L+diff,self.phi_L,self.m1/M_S,self.m2/M_S,self.f_GW,self.mP/M_J,self.P/yr,self.theta_P,self.phi_0,self.T_obs/yr,self.key)
+        elif i == 5:
+            B = binary(self.theta_S_Ec,self.phi_S_Ec,self.dist/pc,self.theta_L,self.phi_L+diff,self.m1/M_S,self.m2/M_S,self.f_GW,self.mP/M_J,self.P/yr,self.theta_P,self.phi_0,self.T_obs/yr,self.key)
+        
+        if j == 2:
+            C = binary(self.theta_S_Ec+diff,self.phi_S_Ec,self.dist/pc,self.theta_L,self.phi_L,self.m1/M_S,self.m2/M_S,self.f_GW,self.mP/M_J,self.P/yr,self.theta_P,self.phi_0,self.T_obs/yr,self.key)
+        elif j == 3:
+            C = binary(self.theta_S_Ec,self.phi_S_Ec+diff,self.dist/pc,self.theta_L,self.phi_L,self.m1/M_S,self.m2/M_S,self.f_GW,self.mP/M_J,self.P/yr,self.theta_P,self.phi_0,self.T_obs/yr,self.key)
+        elif j == 4:
+            C = binary(self.theta_S_Ec,self.phi_S_Ec,self.dist/pc,self.theta_L+diff,self.phi_L,self.m1/M_S,self.m2/M_S,self.f_GW,self.mP/M_J,self.P/yr,self.theta_P,self.phi_0,self.T_obs/yr,self.key)
+        elif j == 5:
+            C = binary(self.theta_S_Ec,self.phi_S_Ec,self.dist/pc,self.theta_L,self.phi_L+diff,self.m1/M_S,self.m2/M_S,self.f_GW,self.mP/M_J,self.P/yr,self.theta_P,self.phi_0,self.T_obs/yr,self.key)
+        
+        return integrate.quad(lambda t: self.d_strain(t,i,B,diff)*self.d_strain(t,j,C,diff),0,self.T_obs,limit=int(self.T_obs*self.f_GW*15))[0]
 
-    def Fisher_mat(self):
-        mat = np.zeros(6,np.float64)
+    def Fisher_mat(self,diff=1e-7):
+        mat = np.zeros((9,9),np.float64)
         for key in [1,2]:
             self.key = key
-            for n, ind in enumerate([[1,1],[1,2],[1,3],[2,2],[2,3],[3,3]]):
-                i, j = ind
-                mat[n] += integrate.quad(lambda t: self.d_strain(t,i)*self.d_strain(t,j),0,self.T_obs,limit=int(self.T_obs*self.f_GW*5))[0]
-        self.Fisher = 2.*np.array([[mat[0],mat[1],mat[2]],
-                                   [mat[1],mat[3],mat[4]],
-                                   [mat[2],mat[4],mat[5]]])/self.S_n
+            for i in np.arange(9): # 0: lnA, 1: f_0, 2: theta_S_Ec, 3: phi_S_Ec, 4: theta_L, 5: phi_L, 6: K, 7: P, 8: phi_0
+                for j in np.arange(i,9):
+                    mat[i][j] += self.h_ixh_j(i, j,diff)
+                    if j != i:
+                        mat[j][i] = mat[i][j]
+        self.Fisher = 2.*mat/self.S_n
         return self.Fisher
     
     def rel_uncertainty(self):
-        unc = np.sqrt(np.diag(np.linalg.inv(self.Fisher_mat()))) #np.linalg.eigvals(self.Fisher_mat()) for diagonalized
+        unc = np.sqrt(np.abs(np.diag(np.linalg.inv(self.Fisher_mat()))[6:])) #np.linalg.eigvals(self.Fisher_mat()) for diagonalized
         return unc/np.array([self.K,self.P,1.])*self.signal_to_noise()
     
     def signal_to_noise(self):