Commit 2678307a authored by Marcus Haberland's avatar Marcus Haberland
Browse files

Update: Continued work in preperation of Draft v1

parent cba59b31
This diff is collapsed.
......@@ -9,6 +9,7 @@ from scipy.constants import *
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize
from tamanini_data import Tamanini
import LISA
import time
import astropy
......@@ -31,7 +32,7 @@ arr_yr = np.linspace(0,yr,1000)
from binary import binary
from ig_binary import ig_binary
infile = open('dict_binaries.txt','rb')
infile = open('dict_binaries_full.txt','rb')
binaries = pickle.load(infile)
infile.close()
......@@ -42,8 +43,15 @@ infile.close()
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)
def comparison(mode='m', weight=1e-4):
bins = binaries[mode]
def comparison_Tam(mode='l', weight=1e-4):
self.P = pd.read_csv(folder+'blue_'+label_csv+'.csv').to_numpy()
self.K = pd.read_csv(folder+'red_'+label_csv+'.csv').to_numpy()
self.phi = pd.read_csv(folder+'green_'+label_csv+'.csv').to_numpy()
self.Pm = np.stack([self.P[:,0],self.K[:,0],self.phi[:,0]]).mean(0)
self.tamm = np.stack([self.K[:,1],self.P[:,1],self.phi[:,1]])
ig_bins = binaries_ig[mode]
P = np.array([b['pars']['P'] for b in bins])
P_ig = np.array([b['pars']['P'] for b in ig_bins])
......@@ -56,6 +64,9 @@ def comparison(mode='m', weight=1e-4):
if mode == 'm':
b = bins[i]['binary']
vals = np.array([b.K,b.P,1.,1.,1.,1.,1.,np.log(b.a0),b.f_1])
if mode == 'l':
b = bins[i]['binary']
vals = np.array([b.K,b.P,1.,1.,1.,1.,1.,np.log(b.a0),b.f_1,b.f_GW])
for j, Pj in enumerate(P_ig):
if isclose(Pi,Pj):
P_fin.append(Pi)
......@@ -64,7 +75,7 @@ def comparison(mode='m', weight=1e-4):
F = bins[i]['Fisher']/bins[i]['SNR']**2/(1+weight*A/B)/bins[i]['pars']['M_P']**2 + weight*ig_bins[j]['Fisher']/2/(weight*A+B)/ig_bins[j]['pars']['M_P']**2
F_fin.append(F)
Tam_fin.append(np.sqrt(np.diag(np.linalg.inv(np.array(F_fin[-1]))/np.outer(vals,vals))))
#print('SNR_L = {:.3f}, SNR_IG = {:.3f}, total = {:.3f}\n'.format(bins[i]['SNR'],weight*bins[i]['SNR']+))
print('SNR_L = {:.3f}, SNR_IG = {:.3f}, total = {:.3f}\n'.format(bins[i]['SNR'],weight*bins[i]['SNR'](1+bins[i]['Fisher'][7,7]/ig_bins[j]['Fisher'][7,7]),bins[i]['SNR']**2 + (weight*bins[i]['SNR'](1+bins[i]['Fisher'][7,7]/ig_bins[j]['Fisher'][7,7]))**2))
Tam_fin = np.array(Tam_fin)
plt.figure(dpi=300)
......@@ -79,11 +90,75 @@ def comparison(mode='m', weight=1e-4):
plt.xlabel(r'$P$ in yr')
plt.ylabel(r'$\sqrt{(\Gamma^{-1})_{ii}}/\lambda_i\cdot$SNR$\cdot M_P/M_J$')
Tamanini(10)
plt.tight_layout()
plt.grid()
plt.savefig(fig+'rel_uncertainty_added.pdf')
def comparison(mode='l', weight=1e-4):
bins = binaries[mode]
ig_bins = binaries_ig[mode]
P = [b['pars']['P'] for b in bins]
P = np.array(P)
P_ig = np.array([b['pars']['P'] for b in ig_bins])
P_fin = []
F_fin = []
Tam_fin = []
Tam_L = []
for i, Pi in enumerate(P):
if mode == 's':
vals = np.array([bins[i]['pars']['K'],Pi*yr,1.])
if mode == 'm':
b = bins[i]['binary']
vals = np.array([b.K,b.P,1.,1.,1.,1.,1.,np.log(b.a0),b.f_1])
if mode == 'l':
b = bins[i]['binary']
vals = np.array([b.K,b.P,1.,1.,1.,1.,1.,np.log(b.a0),b.f_1,b.f_GW])
for j, Pj in enumerate(P_ig):
if isclose(Pi,Pj):
P_fin.append(Pi)
A = ig_bins[j]['SNR']**2/2
B = bins[i]['SNR']**2/2*bins[i]['binary'].S_n()
SNR = (1+weight**2)*bins[i]['SNR']**2
FIG = ig_bins[j]['Fisher']/2 * weight**2 * bins[i]['SNR']**2/A / ig_bins[j]['pars']['M_P']**2 / SNR
FL = bins[i]['Fisher']*bins[i]['binary'].S_n()/2*bins[j]['SNR']**2/B / bins[i]['pars']['M_P']**2 / SNR
F = FIG + FL
F_fin.append(F)
Tam_fin.append(np.sqrt(np.diag(np.linalg.inv(np.array(F_fin[-1]))/np.outer(vals,vals))))
Tam_L.append(bins[i]['Tamanini_plot'])
#print('SNR_L = {:.3f}, SNR_IG = {:.3f}, total = \n'.format(bins[i]['SNR'], weight/Binary.S_n()*bins[i]['Fisher'][7,7])) #weight*bins[i]['SNR']*(1+(bins[i]['Fisher'][7,7]*Binary.S_n()/2)/(ig_bins[j]['Fisher'][7,7]*Binary.S_n()/weight*2)))) #bins[i]['SNR']**2 + (weight*bins[i]['SNR'](1+bins[i]['Fisher'][7,7]/ig_bins[j]['Fisher'][7,7]))**2))
Tam_fin = np.array(Tam_fin)
Tam_L = np.array(Tam_L)
plt.figure(dpi=300)
plt.tight_layout()
# plt.title(r'LISA x IceGiant for $S_n^{LISA}(10 mHZ)/S_n^{IG}(10 mHz) = $'+'{:.1e}'.format(weight))
plt.loglog(P,Tam_L[:,0],'rx')
plt.loglog(P,Tam_L[:,1],'bx')
plt.loglog(P,Tam_L[:,2],'gx')
plt.loglog(P_fin,Tam_fin[:,0],'rs',label=r'$\sigma_K/K$')
plt.loglog(P_fin,Tam_fin[:,1],'bs',label=r'$\sigma_P/P$')
plt.loglog(P_fin,Tam_fin[:,2],'gs',label=r'$\sigma_\varphi$')
plt.xlabel(r'$P$ in yr')
plt.ylabel(r'$\sigma_i/\lambda_i\cdot S/N \cdot M_P/M_J$')
Tamanini(10)
plt.tight_layout()
plt.grid()
plt.savefig(fig+'rel_uncertainty_added.png')
plt.savefig(fig+'rel_uncertainty_added_100x.pdf')
comparison(weight=1e-4,mode='m')
fig, axes = plt.figure(dpi=300)
comparison(weight=1e-2,mode='l',axes.flatten[1])
plt.show()
\ No newline at end of file
This diff is collapsed.
......@@ -6,6 +6,7 @@ import LISA
import time
from astropy import units
import astropy.coordinates as coord
from tamanini_data import Tamanini
import pickle
pi = np.pi
......@@ -34,6 +35,11 @@ infile = open('dict_binaries_ig.txt','rb')
binaries_ig = pickle.load(infile)
infile.close()
infile = open('dict_binaries_full.txt','rb')
binaries_full = 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)
......@@ -277,38 +283,37 @@ def is_analytical_truly_better():
plt.savefig(fig+'\d_strain_yr\d_strain_comparison_5_min_{}e-1yr.png'.format(10*add))
plt.show()
def correlation_mat():
corr = [[ 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]]
corr = np.array(corr)
labels=[r'K',r'P',r'$\phi_0$','$f_0$', 'ln(A)',r'$\theta_S$', '$\phi_S$', r'$\theta_L$', '$\phi_L$']
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)
def correlation_mat(delete=True):
infile = open('dict_binaries_ig.txt','rb')
binaries_ig = pickle.load(infile)
infile.close()
latest = binaries_ig['l'][-1]
corr = latest['correlation']
labels=[r'K',r'P',r'$\phi_0$', r'$\theta_S$', '$\phi_S$', r'$\theta_L$', '$\phi_L$', 'ln(A)', '$f_1$', '$f_0$']
rel_unc = np.abs(corr)*latest['SNR']**2
plt.matshow(np.log10(rel_unc))
plt.xticks(np.arange(10),labels=labels,rotation=40)
plt.yticks(np.arange(10),labels=labels)
plt.colorbar()
plt.title(r'log$_{10}($Correlation matrix$\times$SNR$^2\times (M_P / M_J)^2)$')
plt.title('log$_{10}$(Correlation matrix) for IG, SNR = 1 and $M_P$ = '+str(int(latest['pars']['M_P']))+' $M_J$')
plt.savefig(fig+'Correlation_Mat.png')
plt.show()
labels=[r'K',r'P',r'$\phi_0$', 'ln(A)',r'$\theta_S$', '$\phi_S$', r'$\theta_L$', '$\phi_L$']
corr = np.delete(np.delete(corr,3,0),3,1)
rel_unc = np.log10(np.abs(corr*(B.signal_to_noise()*5)**2))
plt.matshow(rel_unc)
plt.xticks(np.arange(8),labels=labels,rotation=40)
plt.yticks(np.arange(8),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_2.png')
plt.show()
pass
if delete:
labels=[r'K',r'P',r'$\phi_0$', r'$\theta_S$', '$\phi_S$', r'$\theta_L$', '$\phi_L$', 'ln(A)', '$f_1$']
l2 = latest['binary'].copy()
l2.mode = 'm'
corr2 = l2.json()['correlation']
rel_unc2 = np.abs(corr2*(latest['SNR']))**2
plt.matshow((rel_unc2 - rel_unc[:-1,:-1]) / rel_unc[:-1,:-1])
plt.xticks(np.arange(9),labels=labels,rotation=40)
plt.yticks(np.arange(9),labels=labels)
plt.colorbar()
plt.title(r'Relative difference between mode l and m for IG$')
plt.savefig(fig+'Correlation_Mat_2.png')
plt.show()
pass
def look_at_binary(B,delete=[],txt='Correlation_Mat_new.png'):
labels=[r'K',r'P',r'$\phi_0$', r'$\theta_S$', '$\phi_S$', r'$\theta_L$', '$\phi_L$','ln(A)',r'$f_1$',r'$f_0$']
......@@ -353,6 +358,9 @@ def plot_all_Tamanini():
plt.loglog(Ps,tams[:,1],'b-')
plt.loglog(Ps,tams[:,2],'g-')
Tamanini(10)
Tamanini(1)
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.grid()
......@@ -445,7 +453,133 @@ def XRayBinary(f=10e-3,n=10):
plt.show()
return snr
def strain_png():
plt.style.use('seaborn-bright')
plt.figure(dpi=300)
plt.tight_layout()
B = binary()
arr_yr = np.linspace(0,1,10**4)
for key in [1,2]:
B.key=key
plt.plot(arr_yr,[B.strain(t) for t in arr_yr*yr],alpha=.6)
plt.legend([r'$\alpha=I$',r'$\alpha=II$'])
plt.xlabel(r'$t$ in yr')
plt.ylabel(r'Strain $h_\alpha (t)$')
plt.grid()
plt.savefig(fig+r'Strain.pdf')
plt.show()
def LISA_png():
plt.figure(dpi=300)
plt.tight_layout()
B = binary()
arr_yr = np.linspace(0,1,10**4)
for key in [1,2]:
B.key=key
plt.plot(arr_yr,[B.strain(t) for t in arr_yr*yr],alpha=.6)
plt.legend([r'$\alpha=I$',r'$\alpha=II$'])
plt.xlabel(r'$t$ in yr')
plt.ylabel(r'Strain $h_\alpha (t)$')
plt.grid()
plt.savefig(fig+r'Strain.pdf')
plt.show()
def integrand_ex_png():
plt.style.use('seaborn-bright')
plt.figure(dpi=300)
plt.tight_layout()
B = binary(P=8.5)
funcA = B.h_i(7)
funcB = B.h_i(1)
arr_yr = np.linspace(0,4,10**5)
for key in [1,2]:
B.key=key
plt.plot(arr_yr,[funcA(t)*funcB(t) for t in arr_yr*yr],alpha=.6)
plt.legend([r'$\alpha=I$',r'$\alpha=II$'])
plt.xlabel(r'$t$ in yr')
plt.ylabel(r'Integrand $\partial_P h_\alpha \times \partial_{\mathrm{ln}(A)} h_\alpha(t)$')
plt.grid()
plt.savefig(fig+r'Fisher_integrand.pdf')
plt.show()
def LISA_IG_s_png():
PI = []
PL = []
tamI = []
tamL = []
for i in binaries_full['l']:
B = i['binary']
PL.append(i['pars']['P'])
tamL.append(B.reduced_fisher_mat(['theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L', 'ln(A)', 'f_1', 'f_0'],False,True))
for iz, i in enumerate(binaries_ig['s']):
if iz < 3:
continue
B = i['binary']
if B.f_GW > 5e-3 and i['pars']['P'] > 0.03:
PI.append(i['pars']['P'])
tamI.append(i['Tamanini_plot'])
tamI = np.array(tamI)
tamL = np.array(tamL)
plt.figure(dpi=300)
plt.tight_layout()
res = [(PI[i],tamI[i]) for i in np.arange(len(PI))]
res.sort(key=lambda x: x[0])
res.pop(1)
PI = [i[0] for i in res]
tamI = np.array([i[1] for i in res])
plt.loglog(PI,tamI[:,0],'r-.',label=r'$\sigma_K/K$')
plt.loglog(PI,tamI[:,1],'b-.',label=r'$\sigma_P/P$')
plt.loglog(PI,tamI[:,2],'g-.',label=r'$\sigma_\varphi$')
plt.loglog(PI,tamI[:,0],'rx',markersize=5)
plt.loglog(PI,tamI[:,1],'bx',markersize=5)
plt.loglog(PI,tamI[:,2],'gx',markersize=5)
plt.loglog(PL,tamL[:,0],'r^')
plt.loglog(PL,tamL[:,1],'b^')
plt.loglog(PL,tamL[:,2],'g^')
Tamanini(10,False)
plt.legend()
plt.grid()
plt.ylabel(r'$\sigma_i/\lambda_i\cdot S/N \cdot M_P/M_J$')
plt.xlabel(r'$P$ in yr')
plt.savefig(fig+'Tamanini_comp.pdf')
def IG_1_png():
PI = []
tamI = []
for i in binaries_ig['s']:
B = i['binary']
if B.f_GW < 5e-3:
PI.append(i['pars']['P'])
tamI.append(i['Tamanini_plot'])
tamI = np.array(tamI)
plt.figure(dpi=300)
plt.tight_layout()
plt.loglog(PI,tamI[:,0],'rx',label=r'$\sigma_K/K$')
plt.loglog(PI,tamI[:,1],'bx',label=r'$\sigma_P/P$')
plt.loglog(PI,tamI[:,2],'gx',label=r'$\sigma_\varphi$')
Tamanini('all',False)
plt.legend()
plt.grid()
plt.ylabel(r'$\sigma_i/\lambda_i\cdot S/N \cdot M_P/M_J$')
plt.xlabel(r'$P$ in yr')
#plt.savefig(fig+'Tamanini_comp.pdf')
IG_1_png()
#LISA_IG_s_png()
#integrand_ex_png()
#strain_png()
#test2()
#uncs = test(20)
#pos_dep()
......
This diff is collapsed.
No preview for this file type
No preview for this file type
# -*- coding: utf-8 -*-
"""
Created on Sun Mar 13 12:42:40 2022
@author: marcu
"""
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
exos = pd.read_csv('data_tamanini/exoplanets.csv',header=0,skiprows=24,na_values='')
print(exos)
confirmed = exos['soltype'].to_numpy() == 'Published Confirmed'
Ps = (exos['pl_orbper'].to_numpy()/365.25)[confirmed]
Ms = exos['pl_bmassj'].to_numpy()[confirmed]
names = exos['pl_name'].to_numpy()[confirmed]
new = np.concatenate(([True],names[1:] != names[:-1]))
transit = exos['tran_flag'].to_numpy()
rv = (exos['rv_flag'].to_numpy() & (transit == 0))[confirmed]
transit = transit[confirmed]
micro = exos['micro_flag'].to_numpy()[confirmed]
imaging = exos['ima_flag'].to_numpy()[confirmed]
misc = (transit + micro + rv + imaging) == 0
plt.figure(dpi=300)
plt.tight_layout()
plt.style.use('seaborn-white')
for label in [[transit,'Transit'],[rv,'Radial Velocity'],[micro,'Microlensing'],[imaging,'Imaging'],[misc,'other']]:
plt.scatter(Ps[(label[0] != 0)],Ms[(label[0] != 0)],label=label[1],marker='.',alpha=.4)
print('{}: {}'.format(label[1],np.sum( label[0])))
plt.xscale('log')
plt.xlabel('Orbital Period $P$ in yr')
plt.yscale('log')
plt.ylabel('Planet mass $M_P$ or lower mass $\sin (i)M_P$ in $M_J$')
plt.grid()
plt.legend()
plt.savefig('Figures/exoplanets.pdf')
\ No newline at end of file
......@@ -5,12 +5,16 @@ Created on Fri Dec 10 15:43:10 2021
@author: marcu
"""
from tempfile import tempdir
from scipy.constants import *
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize
import LISA
import pickle
import mpmath
mpmath.mp.dps = 50
pi = np.pi
AU = 1495978707e2 # m
......@@ -27,14 +31,14 @@ labels=['K','P','phi_0', 'theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L', 'ln(A)',
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'
file = 'dict_binaries_ig.txt'
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:
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):
def __init__(self,theta_S_Ec=np.arccos(.3),phi_S_Ec=5.,dist=1e3,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=0,epsrel=1e-2):
'''
Parameters
----------
......@@ -64,8 +68,6 @@ class ig_binary:
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)
'''
# Save the parameters of the binary system + exoplanet
......@@ -84,7 +86,7 @@ class ig_binary:
# Save the direction specifications of our detector
self.ig_direction = 0. # TODO: Direction earth-neptune
self.T_2 = 30000 # TODO: Change to a list
self.T_2 = 30000 # TODO: Change to a list # Change to 15000 !!!
self.phi_rel = self.phi_S_Ec - self.ig_direction
......@@ -122,13 +124,13 @@ class ig_binary:
'''
Returns the seperation for a given binary system in m
'''
return (G*(self.m1+self.m2)/(pi*freq)**2)**(1/6) #m
return (G*(self.m1+self.m2)/(pi*self.f_GW)**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
return 1. # irrelevant atm, as we report in SNR
def r_s(self,m):
'''
......@@ -209,7 +211,7 @@ class ig_binary:
'''
Returns
-------
The amplitude of the GW signal
'''
self.amp = np.linalg.norm(self.gw_amplitude*self.F_ig())
return self.amp
......@@ -264,7 +266,7 @@ class ig_binary:
'''
Returns
-------
d y / d lambda_i, where:
The function d y(t) / d lambda_i, where:
i | lambda_i
------------
......@@ -285,12 +287,13 @@ class ig_binary:
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+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,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)
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)
print('Error in d_strain')
raise ValueError
......@@ -316,14 +319,13 @@ class ig_binary:
print('Integrating dh/d{} * dh/d{} ...\n'.format(labels[i],labels[j]))
if i == j:
funcA = self.h_i(i,diff)
funcB = None
return self.integral_over_obs(funcA,funcB,0)
return self.integral_over_obs(funcA,None)
if i != j:
funcA = self.h_i(i,diff)
funcB = self.h_i(j,diff)
return self.integral_over_obs(funcA,funcB,min([self.Fisher[i,i],self.Fisher[j,j]]))
return self.integral_over_obs(funcA,funcB)
def integral_over_obs(self,funcA,funcB,comp):
def integral_over_obs(self,funcA,funcB):
'''
Returns
-------
......@@ -334,21 +336,22 @@ class ig_binary:
----------
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
If None, compute \I funcA^2 dt, else compute \I funcA * funcB dt
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)*funcB(t),start,start + 40*day,limit=int(40*day*self.f_GW*5),epsrel=self.epsrel,epsabs=comp*self.epsabs/10)[0]
temp_res = integrate.quad(lambda t: funcA(t)*funcB(t),start,start + 40*day,limit=int(40*day*self.f_GW*50),epsrel=self.epsrel,epsabs=0)
res += temp_res[0]
print('Integral = {:.3e} +- {:.3e}, rel. error: {:.3f}'.format(temp_res[0],temp_res[1],abs(temp_res[1]/temp_res[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))**2,start,start + 40*day,limit=int(40*day*self.f_GW*5),epsrel=self.epsrel,epsabs=0)[0]
temp_res = integrate.quad(lambda t: (funcA(t))**2,start,start + 40*day,limit=int(40*day*self.f_GW*50),epsrel=self.epsrel,epsabs=0)
res += temp_res[0]
print('Integral = {:.3e} +- {:.3e}, rel. error: {:.3f}'.format(temp_res[0],temp_res[1],abs(temp_res[1]/temp_res[0])))
print('\n')
return res
......@@ -426,7 +429,7 @@ class ig_binary:
if sim:
if j < mode_size[mode_sim]:
self.Fisher[i][j] = js['Fisher'][i,j]
self.Fisher[j][i] = js['Fisher'][j,i]
self.Fisher[j][i] = self.Fisher[i][j]
continue
self.Fisher[i][j] = self.h_ixh_j(i, j,diff)
self.Fisher[j][i] = self.Fisher[i][j]
......@@ -437,9 +440,10 @@ class ig_binary: