Commit 1ba4ef0a authored by Marcus Haberland's avatar Marcus Haberland
Browse files

Better integration algorithm + different sizes of the Fisher matrix +...

Better integration algorithm + different sizes of the Fisher matrix + comparison of these sizes in the txt
parent 2f586195
B = binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=5,P=2,theta_P=pi/2,phi_P=pi/2,T_obs=4)
# 0: K, 1: P, 2: phi_0, 3: f_0, 4: ln(A), 5: theta_S_Ec, 6: phi_S_Ec, 7: theta_L, 8: phi_L
Fisher mat (l):
[[ 6.84162681e+15, 1.27514411e+02, -2.16485628e+08, 2.97509370e+18, 2.86264707e+06, -1.91990430e+10, -1.26610993e+10, 1.47461753e+10, 7.06489260e+09],
[ 1.27514411e+02, 1.46076660e-11, -1.00147822e-04, 9.21983563e+04, -3.92704967e-07, -4.42718335e-04, -6.62314061e-04, 3.05810785e-04, 1.46567797e-04],
[-2.16485628e+08, -1.00147822e-04, 1.00392578e+03, -3.21649452e+11, -1.57469501e-01, -9.99287111e+01, 1.09711875e+01, 4.67096748e+02, -8.03315798e+00],
[ 2.97509370e+18, 9.21983563e+04, -3.21649452e+11, 2.43245074e+21, 5.94288319e+09, -8.04259145e+12, -1.01007573e+13, 9.19115263e+12, 4.29120108e+12],
[ 2.86264707e+06, -3.92704967e-07, -1.57469501e-01, 5.94288319e+09, 1.17907708e+04, 7.16663868e+03, -1.94768798e+04, -6.05861551e+03, 1.26426796e+04],
[-1.91990430e+10, -4.42718335e-04, -9.99287111e+01, -8.04259145e+12, 7.16663868e+03, 7.10045005e+05, -4.54637146e+04, -7.67761154e+04, -3.21179034e+04],
[-1.26610993e+10, -6.62314061e-04, 1.09711875e+01, -1.01007573e+13, -1.94768798e+04, -4.54637146e+04, 3.74173021e+06, -5.22526619e+04, -6.58838947e+03],
[ 1.47461753e+10, 3.05810785e-04, 4.67096748e+02, 9.19115263e+12, -6.05861551e+03, -7.67761154e+04, -5.22526619e+04, 5.74558345e+04, 1.55489246e+04],
[ 7.06489260e+09, 1.46567797e-04, -8.03315798e+00, 4.29120108e+12, 1.26426796e+04, -3.21179034e+04, -6.58838947e+03, 1.55489246e+04, 2.91226477e+04]]
inv(F):
[[ 4.12692901e-16, -1.64253770e-03, -1.02962982e-10, -1.64360181e-19, 2.67954170e-11, -1.14006081e-13, -3.84497921e-14, -5.36682629e-11, -5.07718078e-11],
[-1.64253770e-03, 4.41495350e+11, 4.53595145e+04, 1.51935035e-06, -9.16734218e+02, -2.37250400e+01, 3.52649459e+01, -2.55814546e+03, -2.89238048e+02],
[-1.02962982e-10, 4.53595145e+04, 5.91470638e-03, 7.41384119e-13, -1.22324054e-04, -8.30372239e-06, 3.42616148e-06, -3.86347865e-04, -5.99211969e-05],
[-1.64360181e-19, 1.51935035e-06, 7.41384119e-13, 1.64355905e-21, 2.70244566e-14, -1.39338900e-14, 8.73622896e-16, -2.17158689e-13, -1.20704565e-13],
[ 2.67954170e-11, -9.16734218e+02, -1.22324054e-04, 2.70244566e-14, 3.77464899e-04, -3.19088474e-06, 2.83430777e-06, 9.29412077e-05, -2.22267141e-04],
[-1.14006081e-13, -2.37250400e+01, -8.30372239e-06, -1.39338900e-14, -3.19088474e-06, 1.81709240e-06, 1.95472244e-08, 3.56307270e-06, 3.68917533e-06],
[-3.84497921e-14, 3.52649459e+01, 3.42616148e-06, 8.73622896e-16, 2.83430777e-06, 1.95472244e-08, 2.97912252e-07, 7.47490706e-07, -1.83650108e-06],
[-5.36682629e-11, -2.55814546e+03, -3.86347865e-04, -2.17158689e-13, 9.29412077e-05, 3.56307270e-06, 7.47490706e-07, 1.07640084e-04, -3.59335464e-05],
[-5.07718078e-11, -2.89238048e+02, -5.99211969e-05, -1.20704565e-13, -2.22267141e-04, 3.68917533e-06, -1.83650108e-06, -3.59335464e-05, 1.85207947e-04]]
corr = inv(F) / np.outer(best_fit,best_fit)
[[ 1.04938198e-03, -4.14987183e-05, -1.04523642e-04, -2.62089676e-11, -8.63712760e-07, -5.78670634e-08, -9.75814863e-09, -2.72408696e-05, -1.28853453e-05],
[-4.14987183e-05, 1.10830223e-04, 4.57524714e-04, 2.40726536e-12, 2.93606046e-07, -1.19652869e-07, 8.89261296e-08, -1.29015354e-05, -7.29359411e-07],
[-1.04523642e-04, 4.57524714e-04, 2.39714021e-03, 4.71979789e-11, 1.57415504e-06, -1.68268596e-06, 3.47142736e-07, -7.82904460e-05, -6.07128660e-06],
[-2.62089676e-11, 2.40726536e-12, 4.71979789e-11, 1.64355905e-17, -5.46276459e-14, -4.43529493e-13, 1.39041402e-14, -6.91237576e-12, -1.92107282e-12],
[-8.63712760e-07, 2.93606046e-07, 1.57415504e-06, -5.46276459e-14, 1.54236656e-07, 2.05313147e-08, -9.11848427e-09, -5.98017584e-07, 7.15073872e-07],
[-5.78670634e-08, -1.19652869e-07, -1.68268596e-06, -4.43529493e-13, 2.05313147e-08, 1.84109953e-07, 9.90273954e-10, 3.61014744e-07, 1.86895806e-07],
[-9.75814863e-09, 8.89261296e-08, 3.47142736e-07, 1.39041402e-14, -9.11848427e-09, 9.90273954e-10, 7.54620550e-09, 3.78683215e-08, -4.65191159e-08],
[-2.72408696e-05, -1.29015354e-05, -7.82904460e-05, -6.91237576e-12, -5.98017584e-07, 3.61014744e-07, 3.78683215e-08, 1.09062208e-05, -1.82041472e-06],
[-1.28853453e-05, -7.29359411e-07, -6.07128660e-06, -1.92107282e-12, 7.15073872e-07, 1.86895806e-07, -4.65191159e-08, -1.82041472e-06, 4.69137210e-06]]
unc:
[2.03148443e-08, 6.64451165e+05, 7.69071283e-02, 4.05408318e-11, 1.94284559e-02, 1.34799570e-03, 5.45813386e-04, 1.03749739e-02, 1.36091126e-02]
rel_unc:
[ 3.23941658e-02, 1.05275934e-02, 7.69071283e-02, 4.05408318e-09, -3.92729750e-04, 4.29080357e-04, 8.68688983e-05, 3.30245678e-03, 2.16595755e-03]
---
reduced corr (m):
[[ 1.00171945e-03, -3.57984667e-05, -1.85756475e-05, -6.47519878e-07, -4.35798030e-08, -4.20964220e-05, -1.15747675e-05],
[-3.57984667e-05, 1.09887156e-04, 4.47222798e-04, -9.20004356e-08, 1.04651218e-07, -1.06733315e-05, -1.83544851e-06],
[-1.85756475e-05, 4.47222798e-04, 2.24215149e-03, -6.23139338e-07, 4.09154221e-07, -5.14624834e-05, -8.51768396e-06],
[-6.47519878e-07, -9.20004356e-08, -6.23139338e-07, 1.69783470e-07, 2.48776061e-09, 2.51296825e-07, 4.73866910e-08],
[-4.35798030e-08, 1.04651218e-07, 4.09154221e-07, 2.48776061e-09, 7.00017887e-09, 7.14592156e-09, -3.15940719e-09],
[-4.20964220e-05, -1.06733315e-05, -5.14624834e-05, 2.51296825e-07, 7.14592156e-09, 5.49585121e-06, 2.28340733e-07],
[-1.15747675e-05, -1.83544851e-06, -8.51768396e-06, 4.73866910e-08, -3.15940719e-09, 2.28340733e-07, 1.20669704e-06]]
rel_unc:
[3.16499519e-02, 1.04827075e-02, 7.43793457e-02, 4.12047898e-04, 8.36670716e-05, 2.34432319e-03, 1.09849763e-03]
---
reduced corr (s):
[[ 5.82505837e-04, -1.32769186e-04, -7.57159933e-04],
[-1.32769186e-04, 8.46298498e-05, 5.14886332e-04],
[-7.57159933e-04, 5.14886332e-04, 4.13549118e-03]]
rel_unc:
[0.02413516, 0.00919945, 0.06430778]
---
Tamanini:
~ [0.02547771, 0.01019108, 0.07643312]
l:
[ 3.23941658e-02, 1.05275934e-02, 7.69071283e-02]
m:
[3.16499519e-02, 1.04827075e-02, 7.43793457e-02]
s:
[0.02413516, 0.00919945, 0.06430778]
......@@ -23,6 +23,8 @@ sqth = np.sqrt(3)/2
omega_E = 2*pi/yr
arr_yr = np.linspace(0,yr,1000)
labels=['K','P','phi_0', 'theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L', 'f_0', 'ln(A)']
keys = enumerate(labels)
class binary:
'''
......@@ -215,9 +217,9 @@ class binary:
return (bina.strain(t)-self.strain(t))/diff
if i <= 2:
return -sqth*self.A(t)*(2*pi*self.d_freq_int(t,i+1)+2*pi*self.d_freq(t,i+1)/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))
if i == 3:
if i == 7:
return -sqth*self.A(t)*(2*pi*self.d_freq_int(t,0)+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))
if i == 4:
if i == 8:
return self.strain(t)
print('Error in d_strain')
raise ValueError
......@@ -227,60 +229,115 @@ class binary:
return lambda t: (bina.strain(t)-self.strain(t))/diff
if i <= 2:
return lambda t: -sqth*self.A(t)*(2*pi*self.d_freq_int(t,i+1)+2*pi*self.d_freq(t,i+1)/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))
if i == 3:
if i == 7:
return lambda t: -sqth*self.A(t)*(2*pi*self.d_freq_int(t,0)+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))
if i == 4:
if i == 8:
return lambda t: self.strain(t)
print('Error in d_strain')
raise ValueError
def h_ixh_j(self,i,j,diff=1e-6):
funcA = self.h_i(i,diff)
funcB = self.h_i(j,diff)
return integrate.quad(lambda t: funcA(t)*1e20*funcB(t)*1e20,0,self.T_obs,limit=int(self.T_obs*self.f_GW*20),epsrel=1.49e-2,epsabs=1.49e-2)[0]*1e-40
def h_ixh_j(self,i,j,diff=1e-6,epsrel=1.49e-2,epsabs=1.49e-2):
print('Integrating dh/d{} * dh/d{} ...\n'.format(labels[i],labels[j]))
if i == j:
funcA = self.h_i(i,diff)
funcB = self.h_i(j,diff)
return integrate.quad(lambda t: funcA(t)*1e20*funcB(t)*1e20,0,self.T_obs,limit=int(self.T_obs*self.f_GW*20),epsrel=epsrel,epsabs=0)[0]*1e-40
if i != j:
funcA = self.h_i(i,diff)
funcB = self.h_i(j,diff)
return integrate.quad(lambda t: funcA(t)*1e20*funcB(t)*1e20,0,self.T_obs,limit=int(self.T_obs*self.f_GW*20),epsrel=epsrel,epsabs=np.sqrt(self.Fishers[self.key-1,i,i]*1e40*self.Fishers[self.key-1,j,j]*1e40)*epsabs)[0]*1e-40
def h_i(self,i,diff=1e-6):
B = None
if i == 5:
if i == 3:
if self.theta_S_Ec+diff > pi:
diff *= -1
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 == 6:
elif i == 4:
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 == 7:
elif i == 5:
if self.theta_L+diff > pi:
diff *= -1
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 == 8:
elif i == 6:
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)
return self.d_strain2(i,B,diff)
def Fisher_mat(self,mode=1,diff=1e-6):
if mode == 1:
mat = np.zeros((9,9),np.float64)
for key in [1,2]:
self.key = key
for i in np.arange(9): # 0: K, 1: P, 2: phi_0, 3: f_0, 4: ln(A), 5: theta_S_Ec, 6: phi_S_Ec, 7: theta_L, 8: phi_L
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]
def Fisher_mat(self,mode='l',diff=1e-6,epsrel=1.49e-2,epsabs=1.49e-2,both=False):
if mode == 'l':
self.Fishers = np.zeros((2,9,9),np.float64)
if both:
for key in [1,2]:
self.key = key
for i in np.arange(9): # 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.Fishers[key-1][i][i] = self.h_ixh_j(i, i,diff)
for i in np.arange(9):
for j in np.arange(i+1,9):
self.Fishers[key-1][i][j] = self.h_ixh_j(i, j,diff)
self.Fishers[key-1][j][i] = self.Fishers[key-1][i][j]
else:
self.key = 1
for i in np.arange(9): # 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.Fishers[0][i][i] = self.h_ixh_j(i, i,diff)
for i in np.arange(9):
for j in np.arange(i+1,9):
self.Fishers[0][i][j] = self.h_ixh_j(i, j,diff)
self.Fishers[0][j][i] = self.Fishers[0][i][j]
elif mode == 'm': # without f_0
self.Fishers = np.zeros((2,7,7),np.float64)
if both:
for key in [1,2]:
self.key = key
for i in np.arange(7): # 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.Fishers[key-1][i][i] = self.h_ixh_j(i, i,diff)
for i in np.arange(7):
for j in np.arange(i+1,7):
self.Fishers[key-1][i][j] = self.h_ixh_j(i, j,diff)
self.Fishers[key-1][j][i] = self.Fishers[key-1][i][j]
else:
self.key = 1
for i in np.arange(7): # 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.Fishers[0][i][i] = self.h_ixh_j(i, i,diff)
for i in np.arange(7):
for j in np.arange(i+1,7):
self.Fishers[0][i][j] = self.h_ixh_j(i, j,diff)
self.Fishers[0][j][i] = self.Fishers[0][i][j]
elif mode == 's':
self.Fishers = np.zeros((2,3,3),np.float64)
if both:
for key in [1,2]:
self.key = key
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.Fishers[key-1][i][i] = self.h_ixh_j(i, i,diff)
for i in np.arange(3):
for j in np.arange(i+1,3):
self.Fishers[key-1][i][j] = self.h_ixh_j(i, j,diff)
self.Fishers[key-1][j][i] = self.Fishers[key-1][i][j]
else:
self.key = 1
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.Fishers[0][i][i] = self.h_ixh_j(i, i,diff)
for i in np.arange(3):
for j in np.arange(i+1,3):
self.Fishers[0][i][j] = self.h_ixh_j(i, j,diff)
self.Fishers[0][j][i] = self.Fishers[0][i][j]
else:
mat = np.zeros((3,3),np.float64)
for key in [1,2]:
self.key = key
for i in np.arange(3): # 0: K, 1: P, 2: phi_0, 3: f_0, 4: ln(A), 5: theta_S_Ec, 6: phi_S_Ec, 7: theta_L, 8: phi_L
for j in np.arange(i,3):
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() #2*mat/S_n
print('Please try mode == l for 9x9 matrix, == m for relevant YxY matrix, or == s for 3x3 matrix')
self.Fisher = self.Fishers[0] + self.Fishers[1]
self.Fisher *= 2./self.S_n() #2*mat/S_n
return self.Fisher
def rel_uncertainty(self,mode=1,diff=1e-6):
unc = np.sqrt(np.abs(np.diag(np.linalg.inv(self.Fisher_mat(mode,diff)))[:3])) #np.linalg.eigvals(self.Fisher_mat()) for diagonalized
return unc/np.array([self.K,self.P,1.])
def rel_uncertainty(self,mode='l',diff=1e-6,epsrel=1.49e-2,epsabs=1.49e-2,both=False):
try:
unc = np.sqrt(np.abs(np.diag(np.linalg.inv(self.Fisher_mat(mode,diff)))[:3])) #np.linalg.eigvals(self.Fisher_mat()) for diagonalized
return unc/np.array([self.K,self.P,1.])
except:
print('Something bad happened :c')
return [np.nan,np.nan,np.nan]
def hxh(self):
res = 0.
......
......@@ -3,6 +3,7 @@ import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize
import LISA
import time
pi = np.pi
AU = 1495978707e2 # m
......@@ -19,31 +20,52 @@ omega_E = 2*pi/yr
arr_yr = np.linspace(0,yr,1000)
from binary import binary
keys = enumerate(['K','P','phi_0', 'theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L', 'f_0', 'ln(A)'])
B = binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=5,P=2,theta_P=pi/2,phi_P=pi/2,T_obs=4)
#B.yr_plot('amp')
print(B.f_binary)
#print(B.f_binary)
#B.yr_plot(B.strain,np.linspace(0,yr/1000,1000),1)
#B.yr_plot(B.strain,np.linspace(0,yr/1000,1000),2)
#plt.show()
print(B.signal_to_noise())
print(B.rel_uncertainty())
fig = 'Figures/'
#print(B.signal_to_noise())
start = time.time()
print(B.rel_uncertainty('l'))
end = time.time()
print((end-start)/60)
# fig = 'Figures/'
#print(B.h_ixh_j(6, 6, 1e-6))
def sanity_plots():
add = yr/2
hour = np.linspace(0,5*60,1000) + add
plt.figure(dpi=300)
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)
for key in [1,2]:
B.key=key
C.key=key
plt.plot(hour - add,[B.strain(t) for t in hour],'-')
plt.plot(hour - add,[C.strain(t) for t in hour],':')
plt.legend(['$h_I$ w/ exoplanet','$h_I$ w/o exoplanet','$h_{II}$ w/ exoplanet','$h_{II}$ w/o exoplanet'])
plt.xlabel('$t$ in s for five minutes')
plt.ylabel('Strain $h$ dimensionless')
plt.title(r'For $P=2$ yr and $t_0=0.5$ yr')
plt.savefig(fig+'Strain_2.png')
plt.plot(arr_yr,[B.A(t) for t in arr_yr])
plt.legend(['$A_I$','$A_{II}$'])
plt.xlabel('$t$ in s for one yr')
plt.ylabel('Strain amplitude')
plt.savefig(fig+'Strain_Amplitude.png')
plt.show()
plt.figure(dpi=300)
for key in [1,2]:
B.key=key
plt.plot(arr_yr,[B.polarization_phase(t) for t in arr_yr])
plt.legend([r'$\varphi^{p}_I$',r'$\varphi^{p}_{II}$'])
plt.xlabel(r'$t$ in s for one yr')
plt.ylabel(r'Polarization Phase $\varphi_p$ in rad')
plt.savefig(fig+r'Polarization_Phase.png')
plt.show()
plt.figure(dpi=300)
for key in [1,2]:
B.key=key
plt.plot(arr_yr,[B.doppler_phase(t) for t in arr_yr])
plt.legend([r'$\varphi^{p}_I$',r'$\varphi^{p}_{II}$'])
plt.xlabel(r'$t$ in s for one yr')
plt.ylabel(r'Polarization Phase $\varphi_p$ in rad')
plt.savefig(fig+r'Polarization_Phase.png')
plt.show()
add = 0.
......@@ -62,20 +84,59 @@ def sanity_plots():
plt.savefig(fig+'Strain_1.png')
plt.show()
for i in np.arange(9):
plt.plot(hour,[C.h_i(t,i,B,1e-6) for t in hour])
plt.title(i)
for i, st in keys:
f = C.h_i(i,1e-6)
plt.plot(hour,[f(t) for t in hour+.8*yr])
plt.title(st)
plt.show()
def yr_d_strain(diff=1e-6):
arr_yr = np.linspace(0,1,1000)*yr
for i, st in enumerate(['K','P','phi_0', 'f_0', 'ln(A)', 'theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L']):
plt.figure(dpi=300)
f = B.h_i(i,diff)
plt.plot(arr_yr/yr,[f(t) for t in arr_yr])
#plt.legend([r'$\varphi^{p}_I$',r'$\varphi^{p}_{II}$'])
plt.title('One year strain for $\lambda_i=${}'.format(st))
plt.xlabel(r'$t$ in yr')
plt.ylabel(r'$\frac{\partial h}{\partial \lambda_i}$')
#plt.savefig(fig+'\d_strain_yr\d_strain_{}.png'.format(st))
plt.show()
for i, st in enumerate(['K','P','phi_0', 'f_0', 'ln(A)', 'theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L']):
plt.figure(dpi=300)
f = B.h_i(i,diff)
plt.plot(arr_yr/yr,[f(t)**2 for t in arr_yr])
#plt.legend([r'$\varphi^{p}_I$',r'$\varphi^{p}_{II}$'])
plt.title('One year strain for $\lambda_i=${}'.format(st))
plt.xlabel(r'$t$ in yr')
plt.ylabel(r'$\left(\frac{\partial h}{\partial \lambda_i}\right)^2$')
#plt.savefig(fig+'\d_strain_yr\d_strain^2_{}.png'.format(st))
plt.show()
def uncertainty_plot(n=100,mP=1,a=-2,b=1,mode=2):
Ps = np.logspace(a,b,n)
uncs = np.zeros((n,3),np.float64)
uncs2 = np.zeros((n,3),np.float64)
def uncertainty_plot(n0=100,mP=1,a=-2,b=1,mode='m'):
Ps = np.logspace(a,b,n0)
uncs = np.zeros((n0,3),np.float64)
uncs2 = np.zeros((n0,3),np.float64)
times = np.zeros((2,n0))
bs1 = []
bs2 = []
for n,P0 in enumerate(Ps):
uncs[n] = binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=mP,P=P0,theta_P=pi/2,phi_P=pi/2,T_obs=4).rel_uncertainty(mode)
uncs2[n] = binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=1e-3,mP=mP,P=P0,theta_P=pi/2,phi_P=pi/2,T_obs=4).rel_uncertainty(mode)
uncs*=np.sqrt(B.S_n()/2)*binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=mP,P=2,theta_P=0,phi_P=pi/2,T_obs=4).signal_to_noise()
uncs2*=np.sqrt(B.S_n()/2)*binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=1e-3,mP=mP,P=2,theta_P=0,phi_P=pi/2,T_obs=4).signal_to_noise()
print('Now working on binary #',n+1,' / ',n0,' (P={:.2f} yr) \n'.format(P0))
start = time.time()
B1 = binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=mP,P=P0,theta_P=pi/2,phi_P=pi/2,T_obs=4)
uncs[n] = B1.rel_uncertainty(mode)
end = time.time()
times[0,n] = (end-start)/60
start = time.time()
B2 = binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=1e-3,mP=mP,P=P0,theta_P=pi/2,phi_P=pi/2,T_obs=4)
uncs2[n] = B2.rel_uncertainty(mode)
end = time.time()
times[1,n] = (end-start)/60
bs1.append(B1)
bs2.append(B2)
uncs*=binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=mP,P=2,theta_P=0,phi_P=pi/2,T_obs=4).signal_to_noise()
uncs2*=binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=1e-3,mP=mP,P=2,theta_P=0,phi_P=pi/2,T_obs=4).signal_to_noise()
uncs*=mP
uncs2*=mP
plt.figure(dpi=300)
......@@ -96,7 +157,7 @@ def uncertainty_plot(n=100,mP=1,a=-2,b=1,mode=2):
plt.xlabel(r'$P$ in yr')
plt.savefig(fig+'Relative_Uncertainties_Zoom.png')
plt.show()
return uncs, uncs2
return uncs, uncs2, bs1, bs2, times
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)
......@@ -119,18 +180,20 @@ def one_year_degeneracy():
def deriv_check():
for add in [0,.5,1]:
B = binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=5,P=1.5,theta_P=pi/2,phi_P=pi/2,T_obs=4)
hour = np.linspace(0,5*60,1000) + add
for i in [1,2,3]:
plt.plot(hour,[B.d_freq(t, i) for t in hour],label='{}'.format(i))
B = binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=5,P=5,theta_P=pi/2,phi_P=pi/2,T_obs=4)
hour = np.linspace(0,5*60,1000) + add*yr
for i in [0,1,2]:
f = B.h_i(i,1e-6)
plt.plot(hour,[f(t) for t in hour],label='{}'.format(i))
#plt.plot(hour,[B.f_0(t) for t in hour],label='OG')
plt.legend()
plt.show()
for i in [1,2,3]:
plt.plot(hour,[B.d_freq_int(t, i) for t in hour],label='{}'.format(i))
for i in [3,4,5,6,7,8]:
f = B.h_i(i,1e-6)
plt.plot(hour,[f(t) for t in hour],label='{}'.format(i))
#plt.plot(hour,[B.freq_int(t) for t in hour],label='OG')
plt.legend()
plt.show()
plt.legend()
plt.show()
def pos_dep():
rel_unc = np.zeros((20,20,3))
......@@ -153,9 +216,9 @@ def test(n=100):
uncs = np.zeros((n,9),np.float64)
#uncs2 = np.zeros((n,3),np.float64)
for n,P0 in enumerate(Ps):
uncs[n] = np.diag(binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=10,P=P0,theta_P=pi/2,phi_P=pi/2,T_obs=4).Fisher_mat())
uncs[n] = [binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=10,P=P0,theta_P=pi/2,phi_P=pi/2,T_obs=4).h_ixh_j(i,i,1e-6) for i in np.arange(9)]
#uncs2[n] = binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=1e-3,mP=mP,P=P0,theta_P=pi/2,phi_P=pi/2,T_obs=4).rel_uncertainty()
uncs/=np.diag(B.Fisher)
#uncs/=np.diag(B.Fisher)
#uncs2*=mP
plt.figure(dpi=300)
plt.tight_layout()
......@@ -181,10 +244,17 @@ def test(n=100):
return uncs
#uncs = test()
def test2(n=20,i=5,j=5):
Ps = np.logspace(-2,1,n)
plt.loglog(Ps,[binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=10,P=P0,theta_P=pi/2,phi_P=pi/2,T_obs=4).h_ixh_j(i,i,1e-6) for P0 in Ps])
plt.show()
#test2()
#uncs = test(20)
#pos_dep()
#sanity_plots()
uncs, uncs2 = uncertainty_plot(100,mP=10,a=-2,b=1,mode=2)
#yr_d_strain(1e-6)
uncs, uncs2, bs1, bs2, times = uncertainty_plot(10,mP=10,a=-2,b=1,mode='m')
#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