Commit 487af732 authored by Marcus Haberland's avatar Marcus Haberland
Browse files

Changed the fisher matrix to include computation of the uncertainty in...

Changed the fisher matrix to include computation of the uncertainty in nuisance parameters position, angular momentum, GW frequency
parent 896829ae
......@@ -93,7 +93,6 @@ class binary:
self.phi_L = phi_L
self.m1 = m1*M_S #kg
self.m2 = m2*M_S #kg
self.sep = (G*(self.m1+self.m2)/(pi*freq)**2)**(1/6) #m
self.mP = mP*M_J #kg
self.P = P*yr #s
self.theta_P = theta_P
......@@ -110,6 +109,9 @@ class binary:
self.cos_inclination = np.cos(self.theta_L)*np.cos(self.theta_S_Ec)+np.sin(self.theta_L)*np.sin(self.theta_S_Ec)*np.cos(self.phi_L-self.phi_S_Ec) # = L*n
self.gw_amplitude = 4/self.dist*(G*self.chirp_mass/c**2)**(5/3)*(pi*self.f_GW/c)**(2/3)*np.array([(1.+self.cos_inclination**2)/2, -self.cos_inclination]) # [A_plus, A_cross]
def sep(self):
return (G*(self.m1+self.m2)/(pi*freq)**2)**(1/6) #m
def r_s(self,m):
'''
Returns
......@@ -157,19 +159,8 @@ class binary:
nLxz = .5*np.sin(self.theta_L)*np.sin(self.theta_S_Ec)*np.sin(self.phi_L-self.phi_S_Ec) - sqth*np.cos(omega_E*t)*(np.cos(self.theta_L)*np.sin(self.theta_S_Ec)*np.sin(self.phi_S_Ec)-np.cos(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_L)) - sqth*np.sin(omega_E*t)*(np.cos(self.theta_S_Ec)*np.sin(self.theta_L)*np.cos(self.phi_L) - np.cos(self.theta_L)*np.sin(self.theta_S_Ec)*np.cos(self.phi_S_Ec))
return np.arctan2((Lz - Ln*self.cos_theta_S(t)),nLxz)
def F_LISA_I(self,t=0):
return self.F_LISA_det_frame(self.cos_theta_S(t),self.phi_S(t),self.psi_S(t))
def F_LISA_II(self,t=0):
return self.F_LISA_det_frame(self.cos_theta_S(t),self.phi_S(t)-pi/4,self.psi_S(t))
def F_LISA(self,t=0):
if self.key == 1:
return self.F_LISA_I(t)
if self.key == 2:
return self.F_LISA_II(t)
else:
raise ValueError
return self.F_LISA_det_frame(self.cos_theta_S(t),self.phi_S(t)-pi/4*(self.key-1),self.psi_S(t))
def A(self,t=0):
return np.linalg.norm(self.gw_amplitude*self.F_LISA(t))
......@@ -236,9 +227,10 @@ class binary:
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:
elif i == 1: # TODO
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:
print('Error in d_strain')
raise ValueError
def h_ixh_j(self,i,j,diff=1e-7):
......@@ -262,9 +254,9 @@ class binary:
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]
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*20),points=self.T_obs*np.linspace(0,1,int(self.T_obs*self.f_GW*2)))[0]
def Fisher_mat(self,diff=1e-7):
def Fisher_mat(self,diff=1e-6):
mat = np.zeros((9,9),np.float64)
for key in [1,2]:
self.key = key
......@@ -280,9 +272,12 @@ class binary:
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 hxh(self):
res = 0.
for key in [1,2]:
self.key = key
res += integrate.quad(lambda t: self.strain(t)**2,0,self.T_obs,limit=int(self.T_obs*self.f_GW*20),points=self.T_obs*np.linspace(0,1,int(self.T_obs*self.f_GW*2)))[0]
return res
def signal_to_noise(self):
self.key = 1
Term_1, err_1 = integrate.quad(lambda t: self.strain(t)**2,0,self.T_obs,limit=int(self.T_obs*self.f_GW*5))
self.key = 2
Term_2, err_2 = integrate.quad(lambda t: self.strain(t)**2,0,self.T_obs,limit=int(self.T_obs*self.f_GW*5))
return np.sqrt(2*(Term_1+Term_2)/self.S_n)
\ No newline at end of file
return np.sqrt(2*self.hxh()/self.S_n)
\ No newline at end of file
......@@ -3,7 +3,6 @@ import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize
import LISA
from mpmath import quadosc
pi = np.pi
AU = 1495978707e2 # m
......@@ -127,7 +126,7 @@ def sanity_plots():
plt.show()
def uncertainty_plot(n=100,mP=1):
Ps = np.logspace(-2,1.1,n)
Ps = np.logspace(-2,1,n)
uncs = np.zeros((n,3),np.float64)
uncs2 = np.zeros((n,3),np.float64)
for n,P0 in enumerate(Ps):
......@@ -151,8 +150,9 @@ def uncertainty_plot(n=100,mP=1):
plt.legend()
plt.ylabel(r'$\sqrt{(\Gamma^{-1})_{ii}}/\lambda_i\cdot$SNR$\cdot M_P/M_J$')
plt.xlabel(r'$P$ in yr')
plt.savefig(fig+'Relative_Uncertainties.png')
plt.savefig(fig+'Relative_Uncertainties_Zoom.png')
plt.show()
return uncs, uncs2
def one_year_degeneracy():
C = binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=0,P=100,theta_P=0,phi_P=0,T_obs=4)
......@@ -189,6 +189,7 @@ def deriv_check():
plt.show()
#sanity_plots()
uncertainty_plot(50,mP=1)
uncs, uncs2 = uncertainty_plot(15,mP=10)
#A = np.power(uncs[:-2,:]*uncs[1:-1,:]*uncs[2:,:],1/3)[::3,:]
#one_year_degeneracy()
#deriv_check()
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment