Commit 7522eed2 authored by Marcus Haberland's avatar Marcus Haberland
Browse files

added icegiant data analysis (ig_binary.py) and did some corrections for a final full calculation

parent ccd1b36d
This diff is collapsed.
......@@ -65,7 +65,7 @@ class 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)
'''
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,T_obs=4,mode='m',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=pi/2,T_obs=4,mode='m',epsrel=1.49e-2,epsabs=1.49e-2,key=1):
'''
Parameters
----------
......@@ -116,6 +116,9 @@ class binary:
self.mode = mode
self.key = key
self.epsrel=epsrel
self.epsabs=epsabs
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
......@@ -210,12 +213,12 @@ class binary:
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)))
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.
def strain(self,t):
return sqth*self.A(t)*np.cos(2*pi*self.freq_int(t)+self.polarization_phase(t)+self.doppler_phase(t))
def d_freq(self,t,i): # f_0, K, P, phi_0
def d_freq(self,t,i): # f_0, K, P, phi_0, f_1
if i == 0:
return 1 - self.K*np.cos((2*pi*t)/self.P + self.phi_0)
if i == 1:
......@@ -225,7 +228,7 @@ class binary:
if i == 3:
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)
return t*(1 - self.K*np.cos((2*pi*t)/self.P + self.phi_0))
else:
print('Error in d_freq')
return self.strain(t)
......@@ -260,17 +263,15 @@ class binary:
print('Error in d_strain')
raise ValueError
def h_ixh_j(self,i,j,diff=1e-6,epsrel=1.49e-2,epsabs=1.49e-3):
def h_ixh_j(self,i,j,diff=1e-6):
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
return integrate.quad(lambda t: (funcA(t)*1e20)**2,0,self.T_obs,limit=int(self.T_obs*self.f_GW*20),epsrel=self.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
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=self.epsrel,epsabs=np.sqrt(self.Fishers[self.key-1,i,i]*1e40*self.Fishers[self.key-1,j,j]*1e40)*self.epsabs)[0]*1e-40
def h_i(self,i,diff=1e-6):
B = None
......@@ -289,7 +290,7 @@ class binary:
return self.d_strain(i,B,diff)
def Fisher_mat(self,diff=1e-6,epsrel=1.49e-2,epsabs=1.49e-3,both=False):
def Fisher_mat(self,diff=1e-6,both=False):
sim, js = self.look_for_similar()
if self.mode == 'l':
self.Fishers = np.zeros((2,10,10),np.float64)
......@@ -373,7 +374,7 @@ class binary:
self.Fisher *= 2./self.S_n() #2*mat/S_n
return self.Fisher
def rel_uncertainty(self,diff=1e-6,epsrel=1.49e-2,epsabs=1.49e-3,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']
......@@ -393,7 +394,7 @@ class binary:
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=0)[0]
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=self.epsrel,epsabs=0)[0]
res *= np.sqrt(2)
return res
......@@ -418,8 +419,8 @@ class binary:
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],
'rel_unc' : np.sqrt(np.diag(rel_unc)),
'exo_rel_unc' : np.sqrt(np.diag(rel_unc))[:3],
'SNR' : self.signal_to_noise(),
'binary' : self
}
......@@ -449,6 +450,7 @@ class binary:
for i in dict_binaries:
pars = i['pars']
if isclose(self.P/yr, pars['P']) and isclose(self.phi_0, pars['phi_0']) and isclose(self.theta_S_Ec, pars['theta_S_Ec']) and isclose(self.phi_S_Ec, pars['phi_S_Ec']) and isclose(self.theta_L, pars['theta_L']) and isclose(self.phi_L, pars['phi_L']) and isclose(self.f_GW, pars['f_0']) and isclose(self.mP/M_J, pars['M_P']) and isclose(self.chirp_mass, pars['M_c']):
print('Found pickled data in .txt')
return True, i
return False, None
......
......@@ -32,7 +32,7 @@ labels = ['K','P','phi_0', 'theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L', 'ln(A)'
keys = enumerate(labels)
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')
XRAY = binary(np.arccos(.3),5.,8.6e6,np.arccos(-.2),4.,m1=10,m2=100,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())
......@@ -294,18 +294,16 @@ def plot_all_Tamanini(f0=10e-3,mode='m'):
tam = []
for i in binaries[mode]:
P.append(i['pars']['P'])
f.append(i['pars']['f_0'])
tam.append(i['Tamanini_plot'])
tam = np.array(tam)
P2 = P[f == f0]
tam2 = np.array(tam[f == f0,:])
tam2 = np.array(tam)
P2 = np.array(P)
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.loglog(P2,tam2[:,0],'rx',label=r'$\sigma_K/K$')
plt.loglog(P2,tam2[:,1],'bx',label=r'$\sigma_P/P$')
plt.loglog(P2,tam2[:,2],'gx',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()
......@@ -332,4 +330,4 @@ def look_at_unc():
#correlation_mat()
#is_analytical_truly_better()
#plot_all_Tamanini(mode='m')
look_at_unc()
\ No newline at end of file
#look_at_unc()
\ No newline at end of file
No preview for this file type
This diff is collapsed.
# -*- coding: utf-8 -*-
"""
Created on Tue Dec 14 20:19:19 2021
@author: marcu
"""
import astropy.units as u
import astropy.coordinates as coord
# JPL Horizons: 2040-Dec-01, Neptune-Sun, icrf
ra=2+1/60*9+1/3600*44.81
dec=+11+1/60*13+1/3600*45.1
c = coord.SkyCoord(ra=ra * u.deg,
dec=dec * u.deg,
frame='icrs')
print(c.transform_to(coord.BarycentricMeanEcliptic))
\ No newline at end of file
......@@ -7,29 +7,30 @@ Created on Wed Dec 8 13:50:28 2021
import pickle
import numpy as np
from binary import binary
file = 'dict_binaries.txt'
file = 'dict_binaries_ig.txt'
yr = 24*3600*365.25 # s
infile = open('dict_binaries.txt','rb')
binaries = pickle.load(infile)
infile.close()
if False:
infile = open('dict_binaries.txt','rb')
binaries = pickle.load(infile)
infile.close()
print(binaries)
if True:
outfile = open(file,'wb')
pickle.dump({'s': [], 'm': [],'l':[]},outfile)
outfile.close()
if True:
if False:
#for mode in ['s','m','l']:
mode = 'm'
bins = binaries[mode][:2] + binaries[mode][-2:]
binaries[mode].pop(9)
binaries[mode].pop(6)
bins = binaries[mode]
for bina in bins:
vals = bina['pars']
#print(vals)
B = binary(freq=vals['f_0'],mP=vals['M_P'],P=vals['P'],mode=mode)
B.Fisher = bina['Fisher']
B.snr = bina['SNR']
print(np.log(B.a0) == vals['ln(A)'])
B.add_json()
\ No newline at end of file
print(binaries)
\ No newline at end of file
This diff is collapsed.
......@@ -21,12 +21,12 @@ c = coord.SkyCoord(lat=np.arccos(.3) * u.rad,
Sun = np.array([c.x.value,c.y.value,c.z.value]) # kpc
Centre = coord.SkyCoord(x=0 * u.kpc, y=0* u.kpc, z=0* u.kpc, frame='galactocentric').transform_to(coord.BarycentricMeanEcliptic)
def density_galactic(arr,rho_0=1/794,alpha=.25,R_b2=.64,R_d=2.5,Z_d=.4):
def density_galactic(arr,rho_0=1/29378.,alpha=.25,R_b2=.64,R_d=2.5,Z_d=.4):
x,y,z = arr
u2 = x**2 + y**2
r2 = u2 + z**2
u = np.sqrt(u2)
return rho_0*(alpha*np.exp(-r2/R_b2) + (1-alpha)*np.exp(-u/R_d)/np.cosh(z/Z_d))
return rho_0*(alpha*np.exp(-r2/R_b2) + (1-alpha)*np.exp(-u/R_d)/np.cosh(z/Z_d)**2)
def projected_milky_way(theta_S,phi_S):
new_point = coord.SkyCoord(lat=theta_S * u.rad, lon= phi_S * u.rad, distance= 1. * u.kpc, frame='barycentricmeanecliptic').transform_to(coord.Galactocentric)
......@@ -35,7 +35,7 @@ def projected_milky_way(theta_S,phi_S):
path = lambda s: Sun + s*delta
return quad(lambda s: density_galactic(path(s)),0,np.inf)[0]
return quad(lambda s: s**2*density_galactic(path(s)),0,np.inf)[0]
def generate_skycover(n=100):
thetas = np.arcsin(np.linspace(1,-1,n))
......@@ -51,7 +51,7 @@ def generate_skycover(n=100):
np.savetxt("skycover.csv", vals/V, delimiter=",")
return vals, V
#generate_skycover()
#print(generate_skycover())
Sky = pd.read_csv("skycover.csv",header=None).to_numpy()
plt.figure(dpi=300)
......
......@@ -30,7 +30,7 @@ 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 = binary(np.arccos(.3),5.,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)
n0 = 10
mP = 10
......@@ -42,11 +42,13 @@ 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))
bina = []
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(freq=10e-3,mP=mP,P=P0,theta_P=pi/2,phi_P=pi/2,T_obs=4,mode=mode)
B1 = binary(freq=10e-3,mP=mP,P=P0,theta_P=pi/2,phi_P=pi/2,T_obs=4,mode=mode,epsabs=1.49e-4)
uncs[n] = B1.add_json()['Tamanini_plot']
bina.append(B1.add_json()['binary'])
end = time.time()
times[0,n] = (end-start)/60
print('Finished 10 mHz after {} min'.format(times[0,n]))
......@@ -61,9 +63,9 @@ for n,P0 in enumerate(Ps):
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,uncs[:,0],'rx',label=r'$\sigma_K/K$')
plt.loglog(Ps,uncs[:,1],'bx',label=r'$\sigma_P/P$')
plt.loglog(Ps,uncs[:,2],'gx',label=r'$\sigma_\varphi$')
'''
plt.loglog(Ps,uncs2[:,0],'r--')
......@@ -73,6 +75,7 @@ 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.grid()
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_s.png')
......
# -*- coding: utf-8 -*-
"""
Created on Tue Dec 14 22:04:52 2021
@author: marcu
"""
from scipy.constants import *
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize
import LISA
import time
import pickle
pi = np.pi
AU = 1495978707e2 # m
pc = 3.0857e16 # 2e5 AU
yr = 24*3600*365.25 # s
M_J = 1.898e27 # kg
M_S = 1.989e30 # kg
r_S = 2*G/c**2 # m
sqth = np.sqrt(3)/2
omega_E = 2*pi/yr
fig = 'Figures/'
arr_yr = np.linspace(0,yr,1000)
from ig_binary import ig_binary
keys = enumerate(['K','P','phi_0', 'theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L', 'f_0', 'ln(A)'])
n0 = 15
mP = 10
a = -2
b = 1
mode = 's'
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))
bina = []
for n,P0 in enumerate(Ps):
print('Now working on binary #',n+1,' / ',n0,' (P={:.2f} yr) \n'.format(P0))
start = time.time()
B1 = ig_binary(freq=10e-3,mP=mP,P=P0,dist=1e3,theta_P=pi/2,phi_P=pi/2,mode=mode)
uncs[n] = B1.add_json()['Tamanini_plot']
bina.append(B1.add_json()['binary'])
end = time.time()
times[0,n] = (end-start)/60
print('Finished 10 mHz after {} min'.format(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'.format(times[1,n]))
'''
plt.figure(dpi=300)
plt.tight_layout()
plt.loglog(Ps,uncs[:,0],'rx',label=r'$\sigma_K/K$')
plt.loglog(Ps,uncs[:,1],'bx',label=r'$\sigma_P/P$')
plt.loglog(Ps,uncs[:,2],'gx',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.grid()
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_s.png')
plt.show()
\ 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