Skip to content

GitLab

Commit ac6b99be authored by Marcus Haberland's avatar Marcus Haberland
Browse files

added pickle support to save binaries I already computed and save time

parent aabb1d28
Expand all Hide whitespace changes
Inline Side-by-side
Differentiation.nb
View file @ ac6b99be
This diff is collapsed.
binary.py
View file @ ac6b99be
......@@ -9,6 +9,7 @@ import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize
import LISA
import pickle
pi = np.pi
AU = 1495978707e2 # m
......@@ -23,9 +24,10 @@ 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', 'ln(A)', 'f_0']
labels=['K','P','phi_0', 'theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L', 'ln(A)', 'f_1', 'f_0']
keys = enumerate(labels)
ArcTan = lambda x,y: np.arctan2(y,x)
file = 'dict_binaries.txt'
class binary:
'''
......@@ -57,7 +59,7 @@ class binary:
An array of the timescales of interest, only needed for plots, leave None else
key :
'''
def __init__(self,theta_S_Ec,phi_S_Ec,dist,theta_L,phi_L,m1=.23,m2=.23,freq=2.6e-3,mP=1,P=2,theta_P=pi/2,phi_P=1,T_obs=4,key=1):
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=1,T_obs=4,mode='m',key=1):
'''
Parameters
----------
......@@ -85,6 +87,8 @@ class binary:
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)
'''
......@@ -102,11 +106,14 @@ 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
assert mode in ['s','m','l']
self.mode = mode
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
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.a0 = 4/self.dist*(G*self.chirp_mass/c**2)**(5/3)*(pi*self.f_GW/c)**(2/3)
......@@ -179,7 +186,7 @@ class binary:
return 2*pi*self.freq(t)/c*AU*np.sin(self.theta_S_Ec)*np.cos(omega_E*t-self.phi_S_Ec)
def freq(self,t):
return self.f_GW*(1-self.K*np.cos(2*pi*t/self.P+self.phi_0))
return (self.f_GW + self.f_1*t)*(1 - self.K*np.cos(2*pi*t/self.P + self.phi_0))
def freq_int(self,t):
return self.f_GW*(t-self.K*self.P/(2*pi)*(np.sin(2*pi*t/self.P+self.phi_0)-np.sin(self.phi_0)))
......@@ -189,52 +196,46 @@ class binary:
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))
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)
return -((self.f_GW + self.f_1*t)*np.cos((2*pi*t)/self.P + self.phi_0))
if i == 2:
return -self.f_GW*self.K*2*pi*t/(self.P**2)*np.sin(2*pi*t/self.P+self.phi_0)
return (-2*self.K*pi*t*(self.f_GW + self.f_1*t)*np.sin((2*pi*t)/self.P + self.phi_0))/self.P**2
if i == 3:
return +self.f_GW*self.K*np.sin(2*pi*t/self.P+self.phi_0)
return self.K*(self.f_GW + self.f_1*t)*np.sin((2*pi*t)/self.P + self.phi_0)
if i == 4:
return self.K*(self.f_GW + self.f_1*t)*np.sin((2*pi*t)/self.P + self.phi_0)
else:
print('Error in d_freq')
return self.strain(t)
def d_freq_int(self,t,i):
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((pi*t)/self.P + self.phi_0)*np.sin((pi*t)/self.P))/pi
if i == 1:
return -((self.f_GW*self.P*np.cos((pi*t)/self.P + self.phi_0)*np.sin((pi*t)/self.P))/pi)
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)
if i == 2:
#return +self.f_GW*self.K/self.P*t*np.cos(2*pi/self.P*t+self.phi_0) + self.f_GW*self.K/(2*pi)*(np.sin(2*pi/self.P*t+self.phi_0)-np.sin(self.phi_0))
return self.f_GW*self.K*((t*np.cos((2*pi*t)/self.P + self.phi_0))/self.P - (np.cos((pi*t)/self.P + self.phi_0)*np.sin((pi*t)/self.P))/pi)
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)
if i == 3:
return (self.f_GW*self.K*self.P*np.sin((pi*t)/self.P)*np.sin((pi*t)/self.P + self.phi_0))/pi
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)
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)
else:
print('Error in d_freq_int')
return self.strain(t)
def d_strain(self,t,i,bina=None,diff=1e-6): # 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
if bina is not None:
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 == 8:
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 == 7:
return self.strain(t)
print('Error in d_strain')
raise ValueError
def d_strain2(self,i,bina=None,diff=1e-6): # 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
def d_strain(self,i,bina=None,diff=1e-6): # 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
if bina is not None:
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 == 8:
if i == 9:
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 == 7:
return lambda t: self.strain(t)
if i == 8:
return lambda t: -sqth*self.A(t)*(2*pi*self.d_freq_int(t,4)+2*pi*self.d_freq(t,4)/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))
print('Error in d_strain')
raise ValueError
......@@ -265,50 +266,50 @@ class binary:
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)
return self.d_strain(i,B,diff)
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)
def Fisher_mat(self,diff=1e-6,epsrel=1.49e-2,epsabs=1.49e-2,both=False):
if self.mode == 'l':
self.Fishers = np.zeros((2,10,10),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)
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)
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):
for i in np.arange(10):
for j in np.arange(i+1,10):
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)
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)
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):
for i in np.arange(10):
for j in np.arange(i+1,10):
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)
elif self.mode == 'm': # without f_0
self.Fishers = np.zeros((2,9,9),np.float64)
if both:
for key in [1,2]:
self.key = key
for i in np.arange(8): # 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)
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(8):
for j in np.arange(i+1,8):
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(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)
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(7):
for j in np.arange(i+1,7):
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 == 's':
elif self.mode == 's':
self.Fishers = np.zeros((2,3,3),np.float64)
if both:
for key in [1,2]:
......@@ -328,28 +329,91 @@ class binary:
self.Fishers[0][i][j] = self.h_ixh_j(i, j,diff)
self.Fishers[0][j][i] = self.Fishers[0][i][j]
else:
print('Please try mode == l for 9x9 matrix, == m for relevant 8x8 matrix, or == s for 3x3 matrix')
print('Please try mode == l for 10x10 matrix, == m for relevant 9x9 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='l',diff=1e-6,epsrel=1.49e-2,epsabs=1.49e-2,both=False):
def rel_uncertainty(self,diff=1e-6,epsrel=1.49e-2,epsabs=1.49e-2,both=False):
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_mat(mode,diff)))[:3])) #np.linalg.eigvals(self.Fisher_mat()) for diagonalized
unc = np.sqrt(np.abs(np.diag(np.linalg.inv(self.Fisher))[:3])) #np.linalg.eigvals(self.Fisher_mat()) for diagonalized
self.add_json()
return unc/np.array([self.K,self.P,1.])
except:
print('Something bad happened :c')
return [np.nan,np.nan,np.nan]
unc = np.sqrt(np.abs(np.diag(np.linalg.inv(self.Fisher_mat(diff)))[:3])) #np.linalg.eigvals(self.Fisher_mat()) for diagonalized
self.add_json()
return unc/np.array([self.K,self.P,1.])
def hxh(self):
def hxh(self,both=False):
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))[0]
if both:
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))[0]
else:
res += 1e-20*integrate.quad(lambda t: (self.strain(t)*1e10)**2,0,self.T_obs,limit=int(self.T_obs*self.f_GW*20),epsrel=1.49e-2,epsabs=1e20*1.49e-2)[0]
return res
def signal_to_noise(self):
return np.sqrt(2*self.hxh()/self.S_n())
if not hasattr(self, 'snr'):
self.snr = np.sqrt(2*self.hxh()/self.S_n())
return self.snr
def json(self):
if not hasattr(self, 'Fisher'):
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']
vals = [self.K,self.P,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]
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_0]
rel_unc = err / np.outer(vals2,vals2)
self.js = {'Fisher' : self.Fisher,'Error' : err,
'pars' : {label : vals[i] for i,label in enumerate(labels)},
'rel_unc' : rel_unc,
'exo_rel_unc' : np.diag(rel_unc)[:3],
'SNR' : self.signal_to_noise(),
}
self.js['Tamanini_plot'] = self.js['SNR']*self.js['pars']['M_P']*self.js['exo_rel_unc']
return self.js
def add_json(self):
in_dict, js = self.in_dict_bin()
if not in_dict:
infile = open(file,'rb')
dict_binaries = pickle.load(infile)
infile.close()
dict_binaries[self.mode].append(self.json())
outfile = open(file,'wb')
pickle.dump(dict_binaries,outfile)
outfile.close()
return self.js
else:
return js
def in_dict_bin(self):
infile = open('dict_binaries.txt','rb')
dict_binaries = pickle.load(infile)[self.mode]
infile.close()
for i in dict_binaries:
pars = i['pars']
if self.P == pars['P'] and self.phi_0 == pars['phi_0'] and self.theta_S_Ec == pars['theta_S_Ec'] and self.phi_S_Ec == pars['phi_S_Ec'] and self.theta_L == pars['theta_L'] and self.phi_L == pars['phi_L'] and self.f_GW == pars['f_0'] and self.mP/M_J == pars['M_P'] and self.chirp_mass == pars['M_c']:
return True, i
return False, None
# Analyical stuff
def strain_analytical(self,t):
return (np.sqrt(3)*np.cos(2*pi*self.freq_int(t) + ArcTan((self.a0*(1 + (np.cos(self.theta_S_Ec)*np.cos(self.theta_L) + np.cos(self.phi_S_Ec - self.phi_L)*np.sin(self.theta_S_Ec)*np.sin(self.theta_L))**2)*
......
binary_test.py
View file @ ac6b99be
......@@ -4,6 +4,8 @@ import matplotlib.pyplot as plt
from scipy import integrate, optimize
import LISA
import time
import astropy
import pickle
pi = np.pi
AU = 1495978707e2 # m
......@@ -21,20 +23,21 @@ fig = 'Figures/'
arr_yr = np.linspace(0,yr,1000)
from binary import binary
labels = ['K','P','phi_0', 'theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L', 'f_0', 'ln(A)']
infile = open('dict_binaries.txt','rb')
binaries = pickle.load(infile)
infile.close()
labels = ['K','P','phi_0', 'theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L', 'ln(A)', 'f_1', 'f_0']
keys = enumerate(labels)
B = binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=5,P=1,theta_P=pi/2,phi_P=pi/2,T_obs=4)
#B.yr_plot('amp')
#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())
# start = time.time()
#print(B.rel_uncertainty('m'))
# end = time.time()
# print((end-start)/60)
file = 'dict_binaries.txt'
B = binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=10,P=1,theta_P=pi/2,phi_P=pi/2,T_obs=4,mode='m')
start = time.time()
print(B.rel_uncertainty())
end = time.time()
print((end-start)/60)
B.add_json()
#print(B.h_ixh_j(6, 6, 1e-6))
......@@ -116,53 +119,6 @@ def yr_d_strain(diff=1e-6):
#plt.savefig(fig+'\d_strain_yr\d_strain^2_{}.png'.format(st))
plt.show()
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):
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
print('Finished 10 mHz after {} min'.formar(times[0,n]))
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
print('Finished 1 mHz after {} min'.formar(times[1,n]))
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)
plt.tight_layout()
plt.loglog(Ps,uncs[:,0],'r-',label=r'$\sigma_K/K$')
plt.loglog(Ps,uncs[:,1],'b-',label=r'$\sigma_P/P$')
plt.loglog(Ps,uncs[:,2],'g-',label=r'$\sigma_\varphi$')
plt.loglog(Ps,uncs2[:,0],'r--')
plt.loglog(Ps,uncs2[:,1],'b--')
plt.loglog(Ps,uncs2[:,2],'g--')
plt.title(r'Positions as in [1], $M_{b1,2}=0.23M_\odot,r=1kpc,M_P=1M_J, T_{obs}=4 yr$')
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_Zoom.png')
plt.show()
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)
B = binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=5,P=1,theta_P=pi/2,phi_P=pi/2,T_obs=4)
......@@ -316,6 +272,44 @@ def correlation_mat():
plt.show()
pass
def look_at_binary():
labels=[r'K',r'P',r'$\phi_0$', r'$\theta_S$', '$\phi_S$', r'$\theta_L$', '$\phi_L$','ln(A)',r'$f_1$']
Fish = B.Fisher
Fish_inv = np.linalg.inv(Fish)
vals = np.array([B.K,B.P,1.,.5,1.,.5,1.,np.log(B.a0),B.f_1])
corr = Fish_inv / np.outer(vals,vals)
rel_unc = np.log10(np.abs(corr*(B.signal_to_noise()*5)**2))
plt.matshow(rel_unc)
plt.xticks(np.arange(9),labels=labels,rotation=40)
plt.yticks(np.arange(9),labels=labels)
plt.colorbar()
plt.title(r'log$_{10}($Correlation matrix$\times$SNR$^2\times (M_P / M_J)^2)$')
plt.savefig(fig+'Correlation_Mat.png')
plt.show()
def plot_all_Tamanini(f0=10e-3,mode='s'):
P = []
f = []
tam = []
for i in binaries[mode]:
P.append(i['pars','P'])
f.append(i['pars','f_0'])
tam.append(i['Tamanini_plot'])
P2 = P[f == f0]
tam2 = tam[f == f0]
plt.figure(dpi=300)
plt.tight_layout()
plt.loglog(P2,tam2[:,0],'r-',label=r'$\sigma_K/K$')
plt.loglog(P2,tam2[:,1],'b-',label=r'$\sigma_P/P$')
plt.loglog(P2,tam2[:,2],'g-',label=r'$\sigma_\varphi$')
plt.title(r'Positions as in [1], $M_{b1,2}=0.23M_\odot,r=1kpc,M_P=1M_J, T_{obs}=4 yr$')
plt.legend()
plt.ylabel(r'$\sqrt{(\Gamma^{-1})_{ii}}/\lambda_i\cdot$SNR$\cdot M_P/M_J$')
plt.xlabel(r'$P$ in yr')
#test2()
#uncs = test(20)
#pos_dep()
......@@ -326,4 +320,4 @@ def correlation_mat():
#one_year_degeneracy()
#deriv_check()
#correlation_mat()
is_analytical_truly_better()
\ No newline at end of file
#is_analytical_truly_better()
\ No newline at end of file
reset_dict.py 0 → 100644
View file @ ac6b99be
# -*- coding: utf-8 -*-
"""
Created on Wed Dec 8 13:50:28 2021
@author: marcu
"""
import pickle
file = 'dict_binaries.txt'
if False:
outfile = open(file,'wb')
pickle.dump({'s': [], 'm': [],'l':[]},outfile)
outfile.close()
infile = open('dict_binaries.txt','rb')
binaries = pickle.load(infile)
infile.close()
print(binaries)
\ 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