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

Did some computations and played around with the data analysis after some long computations

parent 1ba4ef0a
This diff is collapsed.
......@@ -108,7 +108,8 @@ class binary:
self.f_GW = freq
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.gw_amplitude = 4/self.dist*(G*self.chirp_mass/c**2)**(5/3)*(pi*self.f_GW/c)**(2/3)*np.array([(1.+self.cos_inclination**2)/2, -self.cos_inclination]) # [A_plus, A_cross]
self.a0 = 4/self.dist*(G*self.chirp_mass/c**2)**(5/3)*(pi*self.f_GW/c)**(2/3)
self.gw_amplitude = self.a0*np.array([(1.+self.cos_inclination**2)/2, -self.cos_inclination]) # [A_plus, A_cross]
def sep(self):
return (G*(self.m1+self.m2)/(pi*freq)**2)**(1/6) #m
......@@ -348,3 +349,60 @@ class binary:
def signal_to_noise(self):
return np.sqrt(2*self.hxh()/self.S_n())
def strain_analytical(self,t):
sqrtth = np.sqrt(3)
return (sqrtth*np.cos(2*pi*self.freq_int(t) + np.arctan(((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))*
(np.cos(2*np.arctan((np.cos(self.theta_L)/2. - (sqrtth*np.cos(self.phi_L - t*omega_E)*np.sin(self.theta_L))/2. -
(np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)*(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)))/
(-0.5*(np.sin(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_S_Ec - self.phi_L)) - (sqrtth*np.cos(t*omega_E)*(np.cos(self.theta_L)*np.sin(self.theta_S_Ec)*np.sin(self.phi_S_Ec) - np.cos(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_L)))/2. -
(sqrtth*(-(np.cos(self.theta_L)*np.cos(self.phi_S_Ec)*np.sin(self.theta_S_Ec)) + np.cos(self.theta_S_Ec)*np.cos(self.phi_L)*np.sin(self.theta_L))*np.sin(t*omega_E))/2.)))*
(np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)*
np.sin(2*(t*omega_E - np.arctan((1./np.sin(self.theta_S_Ec)*1./np.sin(self.phi_S_Ec - t*omega_E)*(sqrtth*np.cos(self.theta_S_Ec) + np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec)))/2.))) +
(np.cos(2*(t*omega_E - np.arctan((1./np.sin(self.theta_S_Ec)*1./np.sin(self.phi_S_Ec - t*omega_E)*(sqrtth*np.cos(self.theta_S_Ec) + np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec)))/2.)))*
(1 + (np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)**2)*
np.sin(2*np.arctan((np.cos(self.theta_L)/2. - (sqrtth*np.cos(self.phi_L - t*omega_E)*np.sin(self.theta_L))/2. -
(np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)*(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)))/
(-0.5*(np.sin(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_S_Ec - self.phi_L)) -
(sqrtth*np.cos(t*omega_E)*(np.cos(self.theta_L)*np.sin(self.theta_S_Ec)*np.sin(self.phi_S_Ec) - np.cos(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_L)))/2. -
(sqrtth*(-(np.cos(self.theta_L)*np.cos(self.phi_S_Ec)*np.sin(self.theta_S_Ec)) + np.cos(self.theta_S_Ec)*np.cos(self.phi_L)*np.sin(self.theta_L))*np.sin(t*omega_E))/2.))))/2.))/
((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)*
((np.cos(2*(t*omega_E - np.arctan((1./np.sin(self.theta_S_Ec)*1./np.sin(self.phi_S_Ec - t*omega_E)*(sqrtth*np.cos(self.theta_S_Ec) + np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec)))/2.)))*
np.cos(2*np.arctan((np.cos(self.theta_L)/2. - (sqrtth*np.cos(self.phi_L - t*omega_E)*np.sin(self.theta_L))/2. -
(np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)*(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)))/
(-0.5*(np.sin(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_S_Ec - self.phi_L)) -
(sqrtth*np.cos(t*omega_E)*(np.cos(self.theta_L)*np.sin(self.theta_S_Ec)*np.sin(self.phi_S_Ec) - np.cos(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_L)))/2. -
(sqrtth*(-(np.cos(self.theta_L)*np.cos(self.phi_S_Ec)*np.sin(self.theta_S_Ec)) + np.cos(self.theta_S_Ec)*np.cos(self.phi_L)*np.sin(self.theta_L))*np.sin(t*omega_E))/2.)))*
(1 + (np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)**2))/2. -
(np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)*
np.sin(2*(t*omega_E - np.arctan((1./np.sin(self.theta_S_Ec)*1./np.sin(self.phi_S_Ec - t*omega_E)*(sqrtth*np.cos(self.theta_S_Ec) + np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec)))/2.)))*
np.sin(2*np.arctan((np.cos(self.theta_L)/2. - (sqrtth*np.cos(self.phi_L - t*omega_E)*np.sin(self.theta_L))/2. -
(np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)*(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)))/
(-0.5*(np.sin(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_S_Ec - self.phi_L)) - (sqrtth*np.cos(t*omega_E)*(np.cos(self.theta_L)*np.sin(self.theta_S_Ec)*np.sin(self.phi_S_Ec) - np.cos(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_L)))/2. -
(sqrtth*(-(np.cos(self.theta_L)*np.cos(self.phi_S_Ec)*np.sin(self.theta_S_Ec)) + np.cos(self.theta_S_Ec)*np.cos(self.phi_L)*np.sin(self.theta_L))*np.sin(t*omega_E))/2.)))))) +
(2*AU*self.freq(t)*pi*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/c)*np.sqrt(self.a0**2*(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*
(np.cos(2*np.arctan((np.cos(self.theta_L)/2. - (sqrtth*np.cos(self.phi_L - t*omega_E)*np.sin(self.theta_L))/2. -
(np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)*(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)))/
(-0.5*(np.sin(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_S_Ec - self.phi_L)) - (sqrtth*np.cos(t*omega_E)*(np.cos(self.theta_L)*np.sin(self.theta_S_Ec)*np.sin(self.phi_S_Ec) - np.cos(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_L)))/2. -
(sqrtth*(-(np.cos(self.theta_L)*np.cos(self.phi_S_Ec)*np.sin(self.theta_S_Ec)) + np.cos(self.theta_S_Ec)*np.cos(self.phi_L)*np.sin(self.theta_L))*np.sin(t*omega_E))/2.)))*
(np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)*
np.sin(2*(t*omega_E - np.arctan((1./np.sin(self.theta_S_Ec)*1./np.sin(self.phi_S_Ec - t*omega_E)*(sqrtth*np.cos(self.theta_S_Ec) + np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec)))/2.))) +
(np.cos(2*(t*omega_E - np.arctan((1./np.sin(self.theta_S_Ec)*1./np.sin(self.phi_S_Ec - t*omega_E)*(sqrtth*np.cos(self.theta_S_Ec) + np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec)))/2.)))*
(1 + (np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)**2)*
np.sin(2*np.arctan((np.cos(self.theta_L)/2. - (sqrtth*np.cos(self.phi_L - t*omega_E)*np.sin(self.theta_L))/2. -
(np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)*(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)))/
(-0.5*(np.sin(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_S_Ec - self.phi_L)) - (sqrtth*np.cos(t*omega_E)*(np.cos(self.theta_L)*np.sin(self.theta_S_Ec)*np.sin(self.phi_S_Ec) - np.cos(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_L)))/2. -
(sqrtth*(-(np.cos(self.theta_L)*np.cos(self.phi_S_Ec)*np.sin(self.theta_S_Ec)) + np.cos(self.theta_S_Ec)*np.cos(self.phi_L)*np.sin(self.theta_L))*np.sin(t*omega_E))/2.))))/2.)**2 +
self.a0**2*(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)**2*
((np.cos(2*(t*omega_E - np.arctan((1./np.sin(self.theta_S_Ec)*1./np.sin(self.phi_S_Ec - t*omega_E)*(sqrtth*np.cos(self.theta_S_Ec) + np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec)))/2.)))*
np.cos(2*np.arctan((np.cos(self.theta_L)/2. - (sqrtth*np.cos(self.phi_L - t*omega_E)*np.sin(self.theta_L))/2. -
(np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)*(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)))/
(-0.5*(np.sin(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_S_Ec - self.phi_L)) - (sqrtth*np.cos(t*omega_E)*(np.cos(self.theta_L)*np.sin(self.theta_S_Ec)*np.sin(self.phi_S_Ec) - np.cos(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_L)))/2. -
(sqrtth*(-(np.cos(self.theta_L)*np.cos(self.phi_S_Ec)*np.sin(self.theta_S_Ec)) + np.cos(self.theta_S_Ec)*np.cos(self.phi_L)*np.sin(self.theta_L))*np.sin(t*omega_E))/2.)))*
(1 + (np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)**2))/2. -
(np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)*
np.sin(2*(t*omega_E - np.arctan((1./np.sin(self.theta_S_Ec)*1./np.sin(self.phi_S_Ec - t*omega_E)*(sqrtth*np.cos(self.theta_S_Ec) + np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec)))/2.)))*
np.sin(2*np.arctan((np.cos(self.theta_L)/2. - (sqrtth*np.cos(self.phi_L - t*omega_E)*np.sin(self.theta_L))/2. -
(np.cos(self.theta_S_Ec)/2. - (sqrtth*np.cos(self.phi_S_Ec - t*omega_E)*np.sin(self.theta_S_Ec))/2.)*(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)))/
(-0.5*(np.sin(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_S_Ec - self.phi_L)) - (sqrtth*np.cos(t*omega_E)*(np.cos(self.theta_L)*np.sin(self.theta_S_Ec)*np.sin(self.phi_S_Ec) - np.cos(self.theta_S_Ec)*np.sin(self.theta_L)*np.sin(self.phi_L)))/2. -
(sqrtth*(-(np.cos(self.theta_L)*np.cos(self.phi_S_Ec)*np.sin(self.theta_S_Ec)) + np.cos(self.theta_S_Ec)*np.cos(self.phi_L)*np.sin(self.theta_L))*np.sin(t*omega_E))/2.))))**2))/2.
\ No newline at end of file
......@@ -16,23 +16,25 @@ 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 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)
labels = ['K','P','phi_0', 'theta_S_Ec', 'phi_S_Ec', 'theta_L', 'phi_L', 'f_0', 'ln(A)']
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('l'))
end = time.time()
print((end-start)/60)
# fig = 'Figures/'
# start = time.time()
#print(B.rel_uncertainty('m'))
# end = time.time()
# print((end-start)/60)
#print(B.h_ixh_j(6, 6, 1e-6))
......@@ -128,11 +130,13 @@ def uncertainty_plot(n0=100,mP=1,a=-2,b=1,mode='m'):
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()
......@@ -249,12 +253,62 @@ def test2(n=20,i=5,j=5):
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()
def is_analytical_truly_better():
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)
arr_yr = np.linspace(0,1,1000)*yr
plt.figure(dpi=300)
plt.plot(arr_yr/yr,[B.strain(t) for t in arr_yr],'-',label='numerical')
plt.plot(arr_yr/yr,[B.strain_analytical(t) for t in arr_yr],':',label='analytic')
plt.legend()
plt.title('One year strain')
plt.xlabel(r'$t$ in yr')
plt.ylabel(r'h(t)')
#plt.savefig(fig+'\d_strain_yr\d_strain_{}.png'.format(st))
plt.show()
hour = np.linspace(0,5*60,1000)
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)
plt.plot(hour,[B.strain(t) for t in hour],'-')
plt.plot(hour,[B.strain_analytical(t) for t in hour],'-')
plt.plot(hour,[C.strain(t) for t in hour],':')
plt.plot(hour,[C.strain_analytical(t) for t in hour],':')
plt.legend(['$h_I$ w/ exoplanet - num','$h_I$ w/ exoplanet - ana','$h_I$ w/o exoplanet - num','$h_I$ w/o exoplanet - ana'])
plt.xlabel('$t$ in s for five minutes')
plt.ylabel('Strain $h$ dimensionless')
plt.title(r'For $P=2$ yr and $t_0=0$ yr')
plt.savefig(fig+'Strain_1.png')
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]]
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)
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()
pass
#test2()
#uncs = test(20)
#pos_dep()
#sanity_plots()
#yr_d_strain(1e-6)
uncs, uncs2, bs1, bs2, times = uncertainty_plot(10,mP=10,a=-2,b=1,mode='m')
#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
#deriv_check()
is_analytical_truly_better()
\ No newline at end of file
# -*- coding: utf-8 -*-
"""
Created on Sun Dec 5 13:19:29 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
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 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)
n0 = 10
mP = 10
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()
\ 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