Commit 74cae427 authored by Marcus Haberland's avatar Marcus Haberland
Browse files

add binary_test.py a playground for the binary.py class

parent c19ac911
from scipy.constants import *
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize
import LISA
from mpmath import quadosc
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
arr_yr = np.linspace(0,yr,1000)
from binary import binary
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.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())
print(B.rel_uncertainty())
fig = 'Figures/'
def sanity_plots():
plt.figure(dpi=300)
for key in [1,2]:
B.key=key
plt.plot(arr_yr,[B.A(t) for t in arr_yr])
plt.legend(['$A_I$','$A_{II}$'])
plt.xlabel('$t$ in s for one yr')
plt.ylabel('Strain amplitude')
plt.savefig(fig+'Strain_Amplitude.png')
plt.show()
plt.figure(dpi=300)
for key in [1,2]:
B.key=key
plt.plot(arr_yr,[B.polarization_phase(t) for t in arr_yr])
plt.legend([r'$\varphi^{p}_I$',r'$\varphi^{p}_{II}$'])
plt.xlabel(r'$t$ in s for one yr')
plt.ylabel(r'Polarization Phase $\varphi_p$ in rad')
plt.savefig(fig+r'Polarization_Phase.png')
plt.show()
plt.figure(dpi=300)
for key in [1,2]:
B.key=key
plt.plot(arr_yr,[B.doppler_phase(t) for t in arr_yr])
plt.legend([r'$\varphi^{D}_I$',r'$\varphi^{D}_{II}$'])
plt.xlabel(r'$t$ in s for one yr')
plt.ylabel(r'Doppler Phase $\varphi_D$ in rad')
plt.savefig(fig+'Doppler_Phase.png')
plt.show()
add = yr/2
hour = np.linspace(0,5*60,1000) + add
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)
for key in [1,2]:
B.key=key
C.key=key
plt.plot(hour - add,[B.strain(t) for t in hour],'-')
plt.plot(hour - add,[C.strain(t) for t in hour],':')
plt.legend(['$h_I$ w/ exoplanet','$h_I$ w/o exoplanet','$h_{II}$ w/ exoplanet','$h_{II}$ w/o exoplanet'])
plt.xlabel('$t$ in s for five minutes')
plt.ylabel('Strain $h$ dimensionless')
plt.title(r'For $P=2$ yr and $t_0=0.5$ yr')
plt.savefig(fig+'Strain_2.png')
plt.show()
add = 0.
hour = np.linspace(0,5*60,1000) + add
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)
for key in [1,2]:
B.key=key
C.key=key
plt.plot(hour - add,[B.strain(t) for t in hour],'-')
plt.plot(hour - add,[C.strain(t) for t in hour],':')
plt.legend(['$h_I$ w/ exoplanet','$h_I$ w/o exoplanet','$h_{II}$ w/ exoplanet','$h_{II}$ w/o exoplanet'])
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()
add = yr
hour = np.linspace(0,5*60,1000) + add
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)
for key in [1,2]:
B.key=key
C.key=key
plt.plot(hour - add,[B.strain(t) for t in hour],'-')
plt.plot(hour - add,[C.strain(t) for t in hour],':')
plt.legend(['$h_I$ w/ exoplanet','$h_I$ w/o exoplanet','$h_{II}$ w/ exoplanet','$h_{II}$ w/o exoplanet'])
plt.xlabel('$t$ in s for five minutes')
plt.ylabel('Strain $h$ dimensionless')
plt.title(r'For $P=2$ yr and $t_0=1$ yr')
plt.savefig(fig+'Strain_3.png')
plt.show()
add = 2*yr
hour = np.linspace(0,5*60,1000) + add
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)
for key in [1,2]:
B.key=key
C.key=key
plt.plot(hour - add,[B.strain(t) for t in hour],'-')
plt.plot(hour - add,[C.strain(t) for t in hour],':')
plt.legend(['$h_I$ w/ exoplanet','$h_I$ w/o exoplanet','$h_{II}$ w/ exoplanet','$h_{II}$ w/o exoplanet'])
plt.xlabel('$t$ in s for five minutes')
plt.ylabel('Strain $h$ dimensionless')
plt.title(r'For $P=2$ yr and $t_0=2$ yr')
plt.savefig(fig+'Strain_4.png')
plt.show()
def uncertainty_plot(n=100,mP=1):
Ps = np.logspace(-2,1.1,n)
uncs = np.zeros((n,3),np.float64)
uncs2 = np.zeros((n,3),np.float64)
for n,P0 in enumerate(Ps):
uncs[n] = 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).rel_uncertainty()
uncs2[n] = 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).rel_uncertainty()
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.png')
plt.show()
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)
for i,adds in enumerate([0,.25,.5,.75,1]):
add = adds*yr
hour = np.linspace(0,5*60,1000) + add
plt.figure(dpi=300)
for key in [1,2]:
B.key=key
C.key=key
plt.plot(hour - add,[B.strain(t) for t in hour],'-')
plt.plot(hour - add,[C.strain(t) for t in hour],':')
plt.legend(['$h_I$ w/ exoplanet','$h_I$ w/o exoplanet','$h_{II}$ w/ exoplanet','$h_{II}$ w/o exoplanet'])
plt.xlabel('$t$ in s for five minutes')
plt.ylabel('Strain $h$ dimensionless')
plt.title(r'For $P=1$ yr and $t_0={}$ yr'.format(adds))
plt.savefig(fig+'OneYearPlanet/Strain_{}.png'.format(i))
plt.show()
def deriv_check():
for add in [0,.5,1]:
B = binary(np.arccos(.3),5.,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=10e-3,mP=5,P=1.5,theta_P=pi/2,phi_P=pi/2,T_obs=4)
hour = np.linspace(0,5*60,1000) + add
for i in [1,2,3]:
plt.plot(hour,[B.d_freq(t, i) for t in hour],label='{}'.format(i))
#plt.plot(hour,[B.f_0(t) for t in hour],label='OG')
plt.legend()
plt.show()
for i in [1,2,3]:
plt.plot(hour,[B.d_freq_int(t, i) for t in hour],label='{}'.format(i))
#plt.plot(hour,[B.freq_int(t) for t in hour],label='OG')
plt.legend()
plt.show()
#sanity_plots()
uncertainty_plot(50,mP=1)
#one_year_degeneracy()
#deriv_check()
\ 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