Commit 80f4c847 authored by Marcus Haberland's avatar Marcus Haberland
Browse files

Commented all of the functions in a hand-wavy manor and fixed some bugs

parent 46193617
......@@ -25,6 +25,7 @@ sqth = np.sqrt(3)/2
labels=['K','P','phi_0', 'theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L', 'ln(A)', 'f_1', 'f_0']
keys = enumerate(labels)
mode_size = {'s' : 3, 'm' : 9, 'l' : 10}
file = r'C:\Users\marcu\Documents\Studium\Kurse\Sem_9\Semester_Project\icegiantexoplanets\dict_binaries_ig.txt'
......@@ -32,39 +33,8 @@ def isclose(a, b, rel_tol=1e-04, abs_tol=0.0):
return abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)
class ig_binary:
'''
Parameters
----------
theta_S_Ec : float in [0,pi]
Inclination angle of the source in ecliptic coordinates in rad
phi_S_Ec : float in [0,2*pi]
Azimuth angle of the source in the ecliptic coordinates in rad
dist : float
Distance from us to binary system in pc
theta_L : float in [0,pi]
The inclination angle of the normal to the binaries orbital plane in ecliptic coords in rad
phi_L : float in [0,2*pi]
The azimuth angle of the normal to the binaries orbital plane in ecliptic coords
m1, m2 : float
The mass of individual binary partners in M_S
freq : float
The frequency of the GW emitted
mP : float
The mass of the single exoplanet in M_J
sepP : float
The seperation from the centre of mass of the binaries to the planet in AU
theta_P : float in [0,pi]
The inclination angle of the normal to the planets orbital plane in ecliptic coords in rad
phi_P : float in [0,2*pi]
The azimuth angle of the normal to the planets orbital plane in ecliptic coords
T_obs : float
Mission duration in years
mode : 's', 'm' or 'l'
Specifies the dimensions of the Fisher-mat: 3x3, 9x9, or 10x10
key : int in [1,2,3]
Used for plots etc: Indication if we should look at length difference I (1), legth difference II (2) or both (3)
'''
def __init__(self,theta_S_Ec=np.arccos(.3),phi_S_Ec=5.,dist=8e3,theta_L=np.arccos(-.2),phi_L=4.,m1=.23,m2=.23,freq=2.6e-3,mP=1,P=2,theta_P=pi/2,phi_P=pi/2,mode='m',epsabs=1.49e-2,epsrel=1.49e-2):
def __init__(self,theta_S_Ec=np.arccos(.3),phi_S_Ec=5.,dist=8e3,theta_L=np.arccos(-.2),phi_L=4.,m1=.23,m2=.23,freq=2.6e-3,mP=1,P=2,theta_P=pi/2,phi_P=pi/2,mode='m',epsabs=1e-2,epsrel=1.49e-2):
'''
Parameters
----------
......@@ -97,6 +67,8 @@ class ig_binary:
key : int in [1,2,3]
Used for plots etc: Indication if we should look at length difference I (1), legth difference II (2) or both (3)
'''
# Save the parameters of the binary system + exoplanet
self.theta_S_Ec = theta_S_Ec
self.phi_S_Ec = phi_S_Ec
self.dist = dist*pc #m
......@@ -110,6 +82,7 @@ class ig_binary:
self.phi_0 = phi_P
self.K = (2*pi*G/self.P)**(1/3)*self.mP/(self.m1+self.m2+self.mP)**(2/3)*np.sin(theta_P)/c
# Save the direction specifications of our detector
self.ig_direction = 0. # TODO: Direction earth-neptune
self.T_2 = 30000 # TODO: Change to a list
......@@ -118,32 +91,43 @@ class ig_binary:
assert mode in ['s','m','l']
self.mode = mode
# Compute relevant parameters of out binary system
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
self.f_1 = 96/5*pi**(8/3)*freq**(11/3)*(G*self.chirp_mass/c**3)**(5/3)
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.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.a0 = 4/self.dist*(G*self.chirp_mass/c**2)**(5/3)*(pi*self.f_GW/c)**(2/3)
self.gw_amplitude = self.a0*np.array([(1.+self.cos_inclination**2)/2, self.cos_inclination]) # [A_plus, A_cross]
# Compute the beam pattern function
#self.strain_amplitude = self.A(t)
self.F = self.F_ig()
# Compute the scalar product k * n = cos(angle los, detector) = mu of detector and line of sight
self.kn = np.sin(theta_S_Ec)*(np.cos(self.ig_direction)*np.cos(phi_S_Ec)+np.sin(self.ig_direction)*np.sin(phi_S_Ec))
print('mu={}'.format(self.kn))
#assert self.kn != 1
assert self.kn != 1
# Set the parameters for numerical integration
self.epsrel = epsrel
self.epsabs = epsabs
# Compute the relevant phases for the GW signal in the detector
self.polarization_phase()
self.A()
def sep(self):
'''
Returns the seperation for a given binary system in m
'''
return (G*(self.m1+self.m2)/(pi*freq)**2)**(1/6) #m
def S_n(self):
'''
Compute the strain sensitivity - currently left out because we report per SNR
'''
return 1. # irrelevant anyways, as we report in SNR
def r_s(self,m):
......@@ -158,7 +142,7 @@ class ig_binary:
'''
Returns
-------
psi_S : The polarization angle psi_S of the wavefront in the detector frame
psi_S : The polarization angle psi_S of the wavefront in the detector frame used to bring the GW to amplitude-phase-form - stolen from Cutler
'''
Lz = .5*np.cos(self.theta_L)-sqth*np.sin(self.theta_L)*np.cos(self.ig_direction-self.phi_L)
Ln = self.cos_inclination
......@@ -169,80 +153,166 @@ class ig_binary:
'''
Returns
-------
[F_plus(t), F_cross(t)] : The orbital frequency of the binary from newtonian physics in Hz, taking into account the two linearly independant signals in arm/key=1, 2
[F_plus(t), F_cross(t)] : The beam pattern function of the specified geometry, already transformed, s.t. we can get to amplitude-phase-form
'''
# beam pattern function for one arm detector, see Maggiore
F = [np.cos(self.theta_S_Ec)**2*np.cos(self.phi_rel)**2-np.sin(self.phi_rel)**2,2*np.cos(self.theta_S_Ec)*np.sin(self.phi_rel)*np.cos(self.phi_rel)]
# beam pattern function for h_xx
return np.array([F[0]*np.cos(2*self.psi_S())-F[1]*np.sin(2*self.psi_S()),
F[0]*np.sin(2*self.psi_S())+F[1]*np.cos(2*self.psi_S())])
# add the beam pattern rotation by polarization angle
return np.array([F[0]*np.cos(2*self.psi_S()) - F[1]*np.sin(2*self.psi_S()),
F[0]*np.sin(2*self.psi_S()) + F[1]*np.cos(2*self.psi_S())])
def nxhxn(self,t=0):
'''
Returns
-------
A(t) : The amplitude of our binary in LISA in the common amplitude-phase-form
h_nn : The strain in direction of the detector
'''
return self.amp*np.cos(2*pi*self.freq_int(t) + self.polarization)
def psi_bar(self,t):
'''
Returns
------
psi_bar : The strain as reported by Armstrong
'''
return self.nxhxn(t)/(1-(self.kn)**2)
def d_psi_bar(self,t):
'''
Returns
-------
A derivative d psi_bar / d lambda = d_psi_bar * d (2*pi*freq_int + polarization) / d lambda - if there is no dependance of the amplitude on lambda
'''
return -self.amp*np.sin(2*pi*self.freq_int(t) + self.polarization)/(1-(self.kn)**2)
def freq(self,t):
'''
Returns
-------
The frequency of the GW emitted by the binary with an exoplanetary doppler shift
'''
return (self.f_GW + self.f_1*t)*(1 - self.K*np.cos(2*pi*t/self.P + self.phi_0))
def polarization_phase(self):
'''
Returns
-------
phi_p(t) : The polarization
polarization : The residual phase shift of the strain in amplitude-phase form
'''
vec = self.gw_amplitude*self.F_ig()
self.polarization = np.arctan2(-vec[1],vec[0])
return self.polarization_phase
return self.polarization
def A(self):
'''
Returns
-------
'''
self.amp = np.linalg.norm(self.gw_amplitude*self.F_ig())
return self.amp
def freq_int(self,t):
return (2*t*(2*self.f_GW + self.f_1*t) + (self.K*self.P*(self.f_1*self.P*np.cos(self.phi_0) - self.f_1*self.P*np.cos((2*pi*t)/self.P + self.phi_0) + 2*self.f_GW*pi*np.sin(self.phi_0) - 2*pi*(self.f_GW + self.f_1*t)*np.sin((2*pi*t)/self.P + self.phi_0)))/pi**2)/4.
'''
Returns
-------
Analytical Integral 0 to t f(t') dt' with f(t) = self.freq(t)
'''
return (2*t*(2*self.f_GW + self.f_1*t) + (self.K*self.P*(self.f_1*self.P*np.cos(self.phi_0) - self.f_1*self.P*np.cos((2*np.pi*t)/self.P + self.phi_0) + 2*self.f_GW*np.pi*np.sin(self.phi_0) - 2*np.pi*(self.f_GW + self.f_1*t)*np.sin((2*np.pi*t)/self.P + self.phi_0)))/np.pi**2)/4.
def strain(self,t):
'''
Returns
-------
y(t) as specified by Armstrong
'''
return (self.kn-1)/2*self.psi_bar(t) - self.kn*self.psi_bar(t-(1+self.kn)/2*self.T_2) + (1+self.kn)/2*self.psi_bar(t-self.T_2)
def d_freq_int(self,t,i):
'''
Returns
-------
The derivative d freq_int / d lambda_i where:
i | lambda_i
------------
0 | f_0
1 | f_1
2 | K
3 | P
4 | phi_0
Verified with Differentiation2.nb and Analytical_Frequencies_IG.txt
'''
if i == 0:
return t - (self.K*self.P*np.cos((pi*t)/self.P + self.phi_0)*np.sin((pi*t)/self.P))/pi
return t - (self.K*self.P*np.cos((np.pi*t)/self.P + self.phi_0)*np.sin((np.pi*t)/self.P))/np.pi
if i == 1:
return (self.P*(self.f_1*self.P*np.cos(self.phi_0) - self.f_1*self.P*np.cos((2*pi*t)/self.P + self.phi_0) + 2*self.f_GW*pi*np.sin(self.phi_0) - 2*pi*(self.f_GW + self.f_1*t)*np.sin((2*pi*t)/self.P + self.phi_0)))/(4.*pi**2)
return t**2/2. - (self.K*self.P*(-(self.P*np.cos(self.phi_0)) + self.P*np.cos((2*np.pi*t)/self.P + self.phi_0) + 2*np.pi*t*np.sin((2*np.pi*t)/self.P + self.phi_0)))/(4.*np.pi**2)
if i == 2:
return (self.K*(self.f_1*self.P**2*np.cos(self.phi_0) + (-(self.f_1*self.P**2) + 2*self.f_GW*pi**2*t + 2*self.f_1*pi**2*t**2)*np.cos((2*pi*t)/self.P + self.phi_0) + self.P*pi*(self.f_GW*np.sin(self.phi_0) - (self.f_GW + 2*self.f_1*t)*np.sin((2*pi*t)/self.P + self.phi_0))))/(2.*self.P*pi**2)
return (self.P*(self.f_1*self.P*np.cos(self.phi_0) - self.f_1*self.P*np.cos((2*np.pi*t)/self.P + self.phi_0) + 2*self.f_GW*np.pi*np.sin(self.phi_0) - 2*np.pi*(self.f_GW + self.f_1*t)*np.sin((2*np.pi*t)/self.P + self.phi_0)))/(4.*np.pi**2)
if i == 3:
return (self.K*self.P*(2*self.f_GW*pi*np.cos(self.phi_0) - 2*pi*(self.f_GW + self.f_1*t)*np.cos((2*pi*t)/self.P + self.phi_0) + self.f_1*self.P*(-np.sin(self.phi_0) + np.sin((2*pi*t)/self.P + self.phi_0))))/(4.*pi**2)
return (self.K*(self.f_1*self.P**2*np.cos(self.phi_0) + (-(self.f_1*self.P**2) + 2*self.f_GW*np.pi**2*t + 2*self.f_1*np.pi**2*t**2)*np.cos((2*np.pi*t)/self.P + self.phi_0) + self.P*np.pi*(self.f_GW*np.sin(self.phi_0) - (self.f_GW + 2*self.f_1*t)*np.sin((2*np.pi*t)/self.P + self.phi_0))))/(2.*self.P*np.pi**2)
if i == 4:
return t**2/2. - (self.K*self.P*(-(self.P*np.cos(self.phi_0)) + self.P*np.cos((2*pi*t)/self.P + self.phi_0) + 2*pi*t*np.sin((2*pi*t)/self.P + self.phi_0)))/(4.*pi**2)
return (self.K*self.P*(2*self.f_GW*np.pi*np.cos(self.phi_0) - 2*np.pi*(self.f_GW + self.f_1*t)*np.cos((2*np.pi*t)/self.P + self.phi_0) + self.f_1*self.P*(-np.sin(self.phi_0) + np.sin((2*np.pi*t)/self.P + self.phi_0))))/(4.*np.pi**2)
else:
print('Error in d_freq_int')
return self.strain(t)
def d_strain(self,i,bina=None,diff=1e-6): # 0: K, 1: P, 2: phi_0, 3: theta_S_Ec, 4: phi_S_Ec, 5: theta_L, 6: phi_L, 7: ln(A), 8: f_1, 9: f_0
'''
Returns
-------
d y / d lambda_i, where:
i | lambda_i
------------
0 | K
1 | P
2 | phi_0
3 | theta_S_Ec
4 | phi__S_Ec
5 | theta_L
6 | phi_L
7 | ln(A)
8 | f_1
9 | f_0
and derivatives 3 .. 6 have to be performed numerically via an instance of bina, which can be set up via self.h_i()
'''
if bina is not None:
return lambda t: (bina.strain(t)-self.strain(t))/diff
if i <= 2:
return lambda t: (self.kn-1)/2*2*pi*self.d_freq_int(t,i)*self.d_psi_bar(t) - self.kn*2*pi*self.d_freq_int(t-(self.kn+1)/2*self.T_2,i)*self.d_psi_bar(t-(self.kn+1)/2*self.T_2) + (self.kn+1)/2*2*pi*self.d_freq_int(t-self.T_2,i)*self.d_psi_bar(t-self.T_2)
return lambda t: (self.kn-1)/2*2*pi*self.d_freq_int(t,i+2)*self.d_psi_bar(t) - self.kn*2*pi*self.d_freq_int(t-(self.kn+1)/2*self.T_2,i+2)*self.d_psi_bar(t-(self.kn+1)/2*self.T_2) + (self.kn+1)/2*2*pi*self.d_freq_int(t-self.T_2,i+2)*self.d_psi_bar(t-self.T_2)
if i == 9:
return lambda t: (self.kn-1)/2*2*pi*self.d_freq_int(t,0)*self.d_psi_bar(t) - self.kn*2*pi*self.d_freq_int(t-(self.kn+1)/2*self.T_2,0)*self.d_psi_bar(t-(self.kn+1)/2*self.T_2) + (self.kn+1)/2*2*pi*self.d_freq_int(t-self.T_2,0)*self.d_psi_bar(t-self.T_2)
if i == 7:
return lambda t: self.strain(t)
if i == 8:
return lambda t: (self.kn-1)/2*2*pi*self.d_freq_int(t,4)*self.d_psi_bar(t) - self.kn*2*pi*self.d_freq_int(t-(self.kn+1)/2*self.T_2,4)*self.d_psi_bar(t-(self.kn+1)/2*self.T_2) + (self.kn+1)/2*2*pi*self.d_freq_int(t-self.T_2,4)*self.d_psi_bar(t-self.T_2)
return lambda t: (self.kn-1)/2*2*pi*self.d_freq_int(t,1)*self.d_psi_bar(t) - self.kn*2*pi*self.d_freq_int(t-(self.kn+1)/2*self.T_2,1)*self.d_psi_bar(t-(self.kn+1)/2*self.T_2) + (self.kn+1)/2*2*pi*self.d_freq_int(t-self.T_2,1)*self.d_psi_bar(t-self.T_2)
print('Error in d_strain')
raise ValueError
def h_ixh_j(self,i,j,diff=1e-6):
'''
Returns
-------
The relevant integral for the Fisher matrix: $\I dy/dlambda_i(t) \cdot dy/dlambda_j(t) dt$, where
i | lambda_i
------------
0 | K
1 | P
2 | phi_0
3 | theta_S_Ec
4 | phi__S_Ec
5 | theta_L
6 | phi_L
7 | ln(A)
8 | f_1
9 | f_0
'''
print('Integrating dh/d{} * dh/d{} ...\n'.format(labels[i],labels[j]))
if i == j:
funcA = self.h_i(i,diff)
......@@ -251,33 +321,69 @@ class ig_binary:
if i != j:
funcA = self.h_i(i,diff)
funcB = self.h_i(j,diff)
return self.integral_over_obs(funcA,funcB,np.sqrt(self.Fisher[i,i]*1e40*self.Fisher[j,j]*1e40))
return self.integral_over_obs(funcA,funcB,min([self.Fisher[i,i],self.Fisher[j,j]]))
def integral_over_obs(self,funcA,funcB,comp):
'''
Returns
-------
The integral over the observational time for the ice giant mission (40 days each year for 10 years) of one or two functions:
\I funcA * funcB dt or \I funcA^2 dt
Parameters
----------
funcA : function (t in sec)
funcB : function (t in sec) or None
If None, compute \I funcA^2 dt, else compute \I funcA * funcB
Comp : float
A scaling for the epsabs if we are on the off-diagonal of the Fisher-matrix, where sometimes parameters are weekly correlated, i.e. we expect Int ~ 0
'''
res = 0.
if funcB is not None:
for i, start in enumerate(np.linspace(0,10,10)*yr):
print('-')
res += integrate.quad(lambda t: funcA(t)*1e20*funcB(t)*1e20,start,start + 40*day,limit=int(40*day*self.f_GW*5),epsrel=self.epsrel,epsabs=comp*self.epsabs/10)[0]*1e-40
res += integrate.quad(lambda t: funcA(t)*funcB(t),start,start + 40*day,limit=int(40*day*self.f_GW*5),epsrel=self.epsrel,epsabs=comp*self.epsabs/10)[0]
print('\n')
return res
else:
for i, start in enumerate(np.linspace(0,10,10)*yr):
print('-')
res += integrate.quad(lambda t: (funcA(t)*1e20)**2,start,start + 40*day,limit=int(40*day*self.f_GW*5),epsrel=self.epsrel,epsabs=comp*self.epsabs/10)[0]*1e-40
res += integrate.quad(lambda t: (funcA(t))**2,start,start + 40*day,limit=int(40*day*self.f_GW*5),epsrel=self.epsrel,epsabs=0)[0]
print('\n')
return res
def h_i(self,i,diff=1e-6):
'''
Returns
-------
A function d y / d lambda_i(t), where:
i | lambda_i
------------
0 | K
1 | P
2 | phi_0
3 | theta_S_Ec
4 | phi_S_Ec
5 | theta_L
6 | phi_L
7 | ln(A)
8 | f_1
9 | f_0
and derivatives 3 .. 6 will be performed numerically with a setup of a close binary system in parameter space
'''
B = None
if i == 3:
if self.theta_S_Ec+diff > pi:
# we care for an overflow in theta_S_Ec and just swap to the left-sided derivative
diff *= -1
B = ig_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,'m')
elif i == 4:
B = ig_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,'m')
elif i == 5:
if self.theta_L+diff > pi:
# we care for an overflow in theta_S_Ec and just swap to the left-sided derivative
diff *= -1
B = ig_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,'m')
elif i == 6:
......@@ -286,111 +392,165 @@ class ig_binary:
return self.d_strain(i,B,diff)
def Fisher_mat(self,diff=1e-6):
sim, js = self.look_for_similar()
if self.mode == 'l':
self.Fisher = np.zeros((10,10),np.float64)
for i in np.arange(10): # 0: K, 1: P, 2: phi_0, 3: theta_S_Ec, 4: phi_S_Ec, 5: theta_L, 6: phi_L, 7: f_0, 8: ln(A)
if sim and i < 3:
self.Fisher[i][i] = js['Fisher'][i,i]
else:
self.Fisher[i][i] = self.h_ixh_j(i, i,diff)
for i in np.arange(10):
for j in np.arange(i+1,10):
if sim and i < 3 and j < 3:
self.Fisher[i][j] = js['Fisher'][i,j]
self.Fisher[i][j] = js['Fisher'][i,j]
else:
self.Fisher[i][j] = self.h_ixh_j(i, j,diff)
self.Fisher[j][i] = self.Fisher[i][j]
elif self.mode == 'm': # without f_0
self.Fisher = np.zeros((9,9),np.float64)
for i in np.arange(9):
if sim and i < 3:
self.Fisher[i][i] = js['Fisher'][i,i]
else:
self.Fisher[i][i] = self.h_ixh_j(i, i,diff)
for i in np.arange(9):
for j in np.arange(i+1,9):
if sim and i < 3 and j < 3:
self.Fisher[i][j] = js['Fisher'][i,j]
self.Fisher[i][j] = js['Fisher'][i,j]
else:
self.Fisher[i][j] = self.h_ixh_j(i, j,diff)
self.Fisher[j][i] = self.Fisher[i][j]
elif self.mode == 's':
'''
Returns
-------
The Fisher matrix Gamma[i,j] = $\I_{T_{obs}} dy/d\lambda_i cdot dy/d\lambda_j dt$ in dependance of the mode:
mode | Fisher matrix
---------------------
's' | (3,3) matrix of only exoplanet parameters: lambda = (K, P, phi_0)
'm' | (9,9) matrix of the uncertain parameters: lambda = (K, P, phi_0, theta_S_Ec, phi_S_Ec, theta_L, phi_L, ln(A), f_1)
'l' | (10,10) matrix of all parameters: lambda = (K, P, phi_0, theta_S_Ec, phi_S_Ec, theta_L, phi_L, ln(A), f_1, f_0)
Also writes the Fisher matrix into self.Fisher
'''
# First, if we didn't already compute the Fisher matrix in this binaries_ig instance, compute it
if not hasattr(self,'Fisher'):
# Then, look if we already performed this computation in dict_binaries_ig.txt
sim, js, mode_sim = self.look_for_similar()
if sim:
self.Fisher = js['Fisher'][:3,:3]
return js['Fisher'][:3,:3]
self.Fisher = np.zeros((3,3),np.float64)
for i in np.arange(3): # 0: K, 1: P, 2: phi_0, 3: theta_S_Ec, 4: phi_S_Ec, 5: theta_L, 6: phi_L, 7: f_0, 8: ln(A),
self.snr = js['SNR']
size = mode_size[self.mode]
self.Fisher = np.zeros((size,size),np.float64)
for i in np.arange(size):
if sim:
if i < mode_size[mode_sim]:
self.Fisher[i][i] = js['Fisher'][i,i]
continue
self.Fisher[i][i] = self.h_ixh_j(i, i,diff)
for i in np.arange(3):
for j in np.arange(i+1,3):
for i in np.arange(size):
for j in np.arange(i+1,size):
if sim:
if j < mode_size[mode_sim]:
self.Fisher[i][j] = js['Fisher'][i,j]
self.Fisher[j][i] = js['Fisher'][j,i]
continue
self.Fisher[i][j] = self.h_ixh_j(i, j,diff)
self.Fisher[j][i] = self.Fisher[i][j]
else:
print('Please try mode == l for 10x10 matrix, == m for relevant 9x9 matrix, or == s for 3x3 matrix')
self.Fisher *= 2/self.S_n() # for signal to noise with S_n = 1 as it is irrelevant
self.Fisher *= 2 #Account for two measurements
self.Fisher *= 2/self.S_n() # for signal to noise with S_n = 1 as it is irrelevant
self.Fisher *= 2 #Account for two measurements / two missions approximately
return self.Fisher
def rel_uncertainty(self,diff=1e-6,epsrel=1.49e-2,epsabs=1.49e-2,both=False):
def rel_uncertainty(self,diff=1e-6):
'''
Returns
-------
[sigma_K/K, sigma_P/P, sigma_phi]The relative uncertainties of the exoplanet parameters of interest, cf. Tamanini (2020)
Parameters
----------
diff : float
The difference for the numerical differentiation (y(par + diff) - y(par)) / diff
'''
# Check if already saved in dict_binaries_ig.txt
in_dict, js = self.in_dict_bin()
if in_dict:
return js['Tamanini_plot']
try:
unc = np.sqrt(np.abs(np.diag(np.linalg.inv(self.Fisher))[:3])) #np.linalg.eigvals(self.Fisher_mat()) for diagonalized
unc = np.sqrt(np.diag(np.linalg.inv(self.Fisher_mat(diff))[:3]))
self.add_json()
return unc/np.array([self.K,self.P,1.])
except:
unc = np.sqrt(np.abs(np.diag(np.linalg.inv(self.Fisher_mat(diff)))[:3])) #np.linalg.eigvals(self.Fisher_mat()) for diagonalized
print('Error in sqrt! Please try again with higher epsrel or epsabs')
unc = np.sqrt(np.abs(np.diag(np.linalg.inv(self.Fisher_mat(diff)))[:3]))
self.add_json()
return unc/np.array([self.K,self.P,1.])
def hxh(self):
'''
Returns
-------
The integrated response $\I_T_obs y \cdot y dt$
'''
return self.integral_over_obs(lambda t: self.strain(t),None,0)
def signal_to_noise(self):
'''
Returns
-------
The signal-to-noise S/N of the binary over the whole observational time
'''
if not hasattr(self, 'snr'):
self.snr = np.sqrt(2*self.hxh()/self.S_n()*2) # x2 for two spacecrafts
return self.snr
def json(self):
if not hasattr(self, 'Fisher'):
'''
Returns
-------
The important calculated properties of the binary system which are to be saved in dict_binaries_ig.txt
json = {'Fisher' : Fisher matrix (mode),
'Error' : Inverse of the Fisher matrix, i.e. the expected variances,
'pars' : {label : Relevant parameters of the binary system},
'correlation' : Correlation matrix of the parameters, i.e. Errors[i,j]/(lambda_i * lambda_j),
'rel_unc' : Relative Uncertainties on the parameters, i.e. diagonal(Correlation),
'exo_rel_unc' : Relative Uncertainties on the exoplanet parameters (K, P, phi_0), i.e. diagonal(Correlation)[:3],
'SNR' : S/N,
'binary' : self
'Tamanini_plot' : Rescaled Relative Uncertainties on the exoplanet parameters (K, P, phi_0) as in Tamanini (2020), i.e. Uncertainty * S/N * M_P/M_J
}
'''
if not hasattr(self,'js'): # Compute it!
self.Fisher_mat()
err = np.linalg.inv(self.Fisher)
labels = ['K','P','phi_0', 'theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L', 'ln(A)', 'f_1', 'f_0','M_P','M_c','dist']
vals = [self.K,self.P/yr,self.phi_0,self.theta_S_Ec,self.phi_S_Ec,self.theta_L,self.phi_L,np.log(self.a0),self.f_1,self.f_GW,self.mP/M_J,self.chirp_mass,self.dist/pc]
if self.mode == 'm':
vals2 = [self.K,self.P,1.,1.,1.,1.,1.,np.log(self.a0),self.f_1]
if self.mode == 's':
vals2 = [self.K,self.P,1.]
if self.mode == 'l':
vals2 = [self.K,self.P,1.,1.,1.,1.,1.,np.log(self.a0),self.f_1,self.f_GW]
rel_unc = err / np.outer(vals2,vals2)
self.js = {'Fisher' : self.Fisher,'Error' : err,
'pars' : {label : vals[i] for i,label in enumerate(labels)},
'correlation' : rel_unc,
'rel_unc' : np.sqrt(np.abs(np.diag(rel_unc))),
'exo_rel_unc' : np.sqrt(np.abs(np.diag(rel_unc)))[:3],
'SNR' : self.signal_to_noise(),
'binary' : self
}
self.js['Tamanini_plot'] = self.js['SNR']*self.js['pars']['M_P']*self.js['exo_rel_unc']
err = np.linalg.inv(self.Fisher)
labels = ['K','P','phi_0', 'theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L', 'ln(A)', 'f_1', 'f_0','M_P','M_c','dist']
vals = [self.K,self.P/yr,self.phi_0,self.theta_S_Ec,self.phi_S_Ec,self.theta_L,self.phi_L,np.log(self.a0),self.f_1,self.f_GW,self.mP/M_J,self.chirp_mass,self.dist/pc]
if self.mode == 's':
vals2 = [self.K,self.P,1.]
if self.mode == 'm':
vals2 = [self.K,self.P,1.,1.,1.,1.,1.,np.log(self.a0),self.f_1]
if self.mode == 'l':
vals2 = [self.K,self.P,1.,1.,1.,1.,1.,np.log(self.a0),self.f_1,self.f_GW]
rel_unc = err / np.outer(vals2,vals2)
self.js = {'Fisher' :