Skip to content

GitLab

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

finished work on binaries and reproduced the Tamanini plot by decreasing... 

finished work on binaries and reproduced the Tamanini plot by decreasing epsabs to 1.49e-3 (3..4 h per data point on Laptop)
added sckycover to see where we expect most binary systems in the milky way
started work on ig_binary to look for exoplanets via icegiantmission
parent ac6b99be
Expand all Hide whitespace changes
Inline Side-by-side
Figures/Correlation_Mat.png
View replaced file @ ac6b99be View file @ ccd1b36d
Figures/Correlation_Mat.png

12.7 KB | W: | H:

Figures/Correlation_Mat.png

12.2 KB | W: | H:

Figures/Correlation_Mat.png
Figures/Correlation_Mat.png
Figures/Correlation_Mat.png
Figures/Correlation_Mat.png
  • 2-up
  • Swipe
  • Onion skin
Figures/Polarization_Phase.png
View replaced file @ ac6b99be View file @ ccd1b36d
Figures/Polarization_Phase.png

110 KB | W: | H:

Figures/Polarization_Phase.png

93.2 KB | W: | H:

Figures/Polarization_Phase.png
Figures/Polarization_Phase.png
Figures/Polarization_Phase.png
Figures/Polarization_Phase.png
  • 2-up
  • Swipe
  • Onion skin
binary.py
View file @ ccd1b36d
......@@ -29,6 +29,9 @@ keys = enumerate(labels)
ArcTan = lambda x,y: np.arctan2(y,x)
file = 'dict_binaries.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 binary:
'''
Parameters
......@@ -45,8 +48,8 @@ class binary:
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
sep : float
The seperation of the binary partners in AU
freq : float
The frequency of the GW emitted
mP : float
The mass of the single exoplanet in M_J
sepP : float
......@@ -55,11 +58,14 @@ class binary:
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
ts : array, optional
An array of the timescales of interest, only needed for plots, leave None else
key :
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=1,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',key=1):
'''
Parameters
----------
......@@ -173,12 +179,27 @@ class binary:
return np.arctan2((Lz - Ln*self.cos_theta_S(t)),nLxz)
def F_LISA(self,t=0):
'''
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
'''
return self.F_LISA_det_frame(self.cos_theta_S(t),self.phi_S(t)-pi/4*(self.key-1),self.psi_S(t))
def A(self,t=0):
'''
Returns
-------
A(t) : The amplitude of our binary in LISA in the common amplitude-phase-form
'''
return np.linalg.norm(self.gw_amplitude*self.F_LISA(t))
def polarization_phase(self,t=0):
'''
Returns
-------
phi_p(t) : The polarization
'''
vec = self.gw_amplitude*self.F_LISA(t)
return np.arctan2(-vec[1],vec[0])
......@@ -239,7 +260,7 @@ 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-2):
def h_ixh_j(self,i,j,diff=1e-6,epsrel=1.49e-2,epsabs=1.49e-3):
print('Integrating dh/d{} * dh/d{} ...\n'.format(labels[i],labels[j]))
if i == j:
funcA = self.h_i(i,diff)
......@@ -256,19 +277,20 @@ class binary:
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)
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,'m',self.key)
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)
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,'m',self.key)
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)
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,'m',self.key)
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)
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,'m',self.key)
return self.d_strain(i,B,diff)
def Fisher_mat(self,diff=1e-6,epsrel=1.49e-2,epsabs=1.49e-2,both=False):
def Fisher_mat(self,diff=1e-6,epsrel=1.49e-2,epsabs=1.49e-3,both=False):
sim, js = self.look_for_similar()
if self.mode == 'l':
self.Fishers = np.zeros((2,10,10),np.float64)
if both:
......@@ -283,11 +305,18 @@ class binary:
else:
self.key = 1
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)
if sim and i < 3:
self.Fishers[0][i][i] = js['Fisher'][i,i]
else:
self.Fishers[0][i][i] = self.h_ixh_j(i, i,diff)
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]
if sim and i < 3 and j < 3:
self.Fishers[0][i][j] = js['Fisher'][i,j]
self.Fishers[0][i][j] = js['Fisher'][i,j]
else:
self.Fishers[0][i][j] = self.h_ixh_j(i, j,diff)
self.Fishers[0][j][i] = self.Fishers[0][i][j]
elif self.mode == 'm': # without f_0
self.Fishers = np.zeros((2,9,9),np.float64)
......@@ -302,14 +331,23 @@ class binary:
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):
if sim and i < 3:
self.Fishers[0][i][i] = js['Fisher'][i,i]
else:# 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]
if sim and i < 3 and j < 3:
self.Fishers[0][i][j] = js['Fisher'][i,j]
self.Fishers[0][i][j] = js['Fisher'][i,j]
else:
self.Fishers[0][i][j] = self.h_ixh_j(i, j,diff)
self.Fishers[0][j][i] = self.Fishers[0][i][j]
elif self.mode == 's':
if sim:
return js['Fisher'][:3,:3]
self.Fishers = np.zeros((2,3,3),np.float64)
if both:
for key in [1,2]:
......@@ -328,13 +366,14 @@ class binary:
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]
self.Fishers[0] *= 2 #Account for two measurements
else:
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,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-3,both=False):
in_dict, js = self.in_dict_bin()
if in_dict:
return js['Tamanini_plot']
......@@ -354,7 +393,8 @@ 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=1e20*1.49e-2)[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=1.49e-2,epsabs=0)[0]
res *= np.sqrt(2)
return res
def signal_to_noise(self):
......@@ -366,8 +406,8 @@ class binary:
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]
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':
......@@ -377,9 +417,11 @@ class binary:
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],
'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']
return self.js
......@@ -406,12 +448,29 @@ class binary:
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']:
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']):
return True, i
return False, None
def look_for_similar(self):
if self.mode in ['m','l']:
B2 = self.copy()
B2.mode = 's'
return B2.in_dict_bin()
if self.mode == 's':
B2 = self.copy()
B2.mode = 'm'
return B2.in_dict_bin()
def create_from_json(json,mode='m'):
pars = json['pars']
B = binary(theta_S_Ec=pars['theta_S_Ec'],phi_S_Ec=pars['phi_S_Ec'],dist=pars['dist'],theta_L=pars['theta_L'],phi_L=pars['phi_L'],m1=.23,m2=.23,freq=pars['f_0'],mP=pars['M_P'],P=pars['P'],theta_P=pi/2,phi_P=pars['phi_0'],T_obs=4,mode=mode,key=1)
B.snr = json['SNR']*B.a0/np.exp(pars['ln(A)'])
B.Fisher = json['Fisher']
return B
def copy(self):
return binary(theta_S_Ec=self.theta_S_Ec,phi_S_Ec=self.phi_S_Ec,dist=self.dist/pc,theta_L=self.theta_L,phi_L=self.phi_L,m1=self.m1/M_S,m2=self.m2/M_S,freq=self.f_GW,mP=self.mP/M_J,P=self.P/yr,theta_P=self.theta_P,phi_P=self.phi_0,T_obs=self.T_obs/yr,mode=self.mode,key=self.key)
# Analyical stuff
......
binary_test.py
View file @ ccd1b36d
......@@ -33,11 +33,12 @@ 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')
start = time.time()
print(B.rel_uncertainty())
end = time.time()
print((end-start)/60)
B.add_json()
# 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))
......@@ -272,7 +273,7 @@ def correlation_mat():
plt.show()
pass
def look_at_binary():
def look_at_binary(B):
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)
......@@ -287,16 +288,17 @@ def look_at_binary():
plt.savefig(fig+'Correlation_Mat.png')
plt.show()
def plot_all_Tamanini(f0=10e-3,mode='s'):
def plot_all_Tamanini(f0=10e-3,mode='m'):
P = []
f = []
tam = []
for i in binaries[mode]:
P.append(i['pars','P'])
f.append(i['pars','f_0'])
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 = tam[f == f0]
tam2 = np.array(tam[f == f0,:])
plt.figure(dpi=300)
plt.tight_layout()
......@@ -310,6 +312,14 @@ def plot_all_Tamanini(f0=10e-3,mode='s'):
plt.ylabel(r'$\sqrt{(\Gamma^{-1})_{ii}}/\lambda_i\cdot$SNR$\cdot M_P/M_J$')
plt.xlabel(r'$P$ in yr')
def look_at_unc():
infile = open('dict_binaries.txt','rb')
binaries = pickle.load(infile)
infile.close()
for json in binaries['m']:
look_at_binary(json['binary'])
plt.show()
#test2()
#uncs = test(20)
#pos_dep()
......@@ -320,4 +330,6 @@ def plot_all_Tamanini(f0=10e-3,mode='s'):
#one_year_degeneracy()
#deriv_check()
#correlation_mat()
#is_analytical_truly_better()
\ No newline at end of file
#is_analytical_truly_better()
#plot_all_Tamanini(mode='m')
look_at_unc()
\ No newline at end of file
dict_binaries.txt
View file @ ccd1b36d
No preview for this file type
ig_binary.py 0 → 100644
View file @ ccd1b36d
This diff is collapsed.
reset_dict.py
View file @ ccd1b36d
......@@ -5,15 +5,31 @@ Created on Wed Dec 8 13:50:28 2021
@author: marcu
"""
import pickle
import numpy as np
from binary import binary
file = 'dict_binaries.txt'
yr = 24*3600*365.25 # s
if False:
infile = open('dict_binaries.txt','rb')
binaries = pickle.load(infile)
infile.close()
if True:
outfile = open(file,'wb')
pickle.dump({'s': [], 'm': [],'l':[]},outfile)
outfile.close()
infile = open('dict_binaries.txt','rb')
binaries = pickle.load(infile)
infile.close()
if True:
#for mode in ['s','m','l']:
mode = 'm'
bins = binaries[mode][:2] + binaries[mode][-2:]
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()
print(binaries)
\ No newline at end of file
skycover.csv 0 → 100644
View file @ ccd1b36d
This source diff could not be displayed because it is too large. You can view the blob instead.
skycover.py 0 → 100644
View file @ ccd1b36d
# -*- coding: utf-8 -*-
"""
Created on Fri Dec 10 13:36:35 2021
@author: marcu
"""
import astropy.units as u
import astropy.coordinates as coord
import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt
import pandas as pd
fig = 'Figures/'
c = coord.SkyCoord(lat=np.arccos(.3) * u.rad,
lon=5. * u.rad,
distance=0. * u.kpc,
frame='barycentricmeanecliptic').transform_to(coord.Galactocentric)
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):
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))
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)
pointer = np.array([new_point.x.value,new_point.y.value,new_point.z.value])
delta = pointer - Sun
path = lambda s: Sun + s*delta
return quad(lambda s: density_galactic(path(s)),0,np.inf)[0]
def generate_skycover(n=100):
thetas = np.arcsin(np.linspace(1,-1,n))
phis = np.linspace(0,2*np.pi,n)
vals = np.zeros((n,n))
for i, th in enumerate(thetas):
for j, ph in enumerate(phis):
vals[i,j] = projected_milky_way(th,ph)
V = np.sum(vals)
plt.imshow(vals/V)
np.savetxt("skycover.csv", vals/V, delimiter=",")
return vals, V
#generate_skycover()
Sky = pd.read_csv("skycover.csv",header=None).to_numpy()
plt.figure(dpi=300)
plt.imshow(np.log(Sky*.9+.1),extent=[0,2*np.pi,np.pi/2,-np.pi/2])
#plt.colorbar()
plt.title(r'log($\rho$) projected')
plt.xlabel(r'$\phi_S$')
plt.ylabel(r'$\theta_S$ uniform in sin$(\theta_S)$')
plt.xticks()
plt.savefig(fig+'SkyCover.png')
\ No newline at end of file
uncertainty_plot_Tamanini.py
View file @ ccd1b36d
......@@ -11,6 +11,7 @@ import matplotlib.pyplot as plt
from scipy import integrate, optimize
import LISA
import time
import pickle
pi = np.pi
AU = 1495978707e2 # m
......@@ -41,28 +42,22 @@ 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)
B1 = binary(freq=10e-3,mP=mP,P=P0,theta_P=pi/2,phi_P=pi/2,T_obs=4,mode=mode)
uncs[n] = B1.add_json()['Tamanini_plot']
end = time.time()
times[0,n] = (end-start)/60
print('Finished 10 mHz after {} min'.formar(times[0,n]))
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'.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
print('Finished 1 mHz after {} min'.format(times[1,n]))
'''
plt.figure(dpi=300)
plt.tight_layout()
......@@ -70,14 +65,15 @@ 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.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