Commit 1cf15046 authored by Marcus Haberland's avatar Marcus Haberland
Browse files

save this version before I fuck it up with the analytic derivatives

parent d8847367
......@@ -105,13 +105,16 @@ class binary:
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
self.S_n = LISA.LISA().Sn(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]
def sep(self):
return (G*(self.m1+self.m2)/(pi*freq)**2)**(1/6) #m
def S_n(self):
return LISA.LISA().Sn(self.f_GW)
def r_s(self,m):
'''
Returns
......@@ -146,7 +149,7 @@ class binary:
phi_S_det : The azimuth angle phi_S of the source in the detector frame
'''
arctan_term = np.arctan2((np.sqrt(3)*np.cos(self.theta_S_Ec)+np.sin(self.theta_S_Ec)*np.cos(omega_E*t-self.phi_S_Ec)),(2*np.sin(self.theta_S_Ec)*np.sin(omega_E*t-self.phi_S_Ec)))
return (omega_E*t + arctan_term) % (2*pi)
return (omega_E*t - arctan_term) % (2*pi)
def psi_S(self,t=0):
'''
......@@ -165,19 +168,6 @@ class binary:
def A(self,t=0):
return np.linalg.norm(self.gw_amplitude*self.F_LISA(t))
def yr_plot(self,func,ts=np.linspace(0,yr,1000),key=None):
if key is not None:
self.key = key
if func == 'amp':
for key in [1,2]:
self.key = key
plt.plot(ts,[self.A(t) for t in ts])
plt.legend(['AI','AII'])
plt.show()
return
funcs = [func(t) for t in ts]
plt.plot(ts,funcs)
def polarization_phase(self,t=0):
vec = self.gw_amplitude*self.F_LISA(t)
return np.arctan2(-vec[1],vec[0])
......@@ -207,77 +197,95 @@ class binary:
print('Error in d_freq')
return self.strain(t)
def d_freq_int(self,t,i): # TODO: f_1
def d_freq_int(self,t,i):
if i == 0:
return t - (self.K*self.P*np.cos(pi*t)/self.P + self.phi_0)*np.sin((pi*t)/self.P)/pi
if i == 1:
return -self.f_GW/(2*pi)*self.P*(np.sin(2*pi*t/self.P+self.phi_0) - np.sin(self.phi_0))
return -((self.f_GW*self.P*np.cos((pi*t)/self.P + self.phi_0)*np.sin((pi*t)/self.P))/pi)
if i == 2:
#return +self.f_GW*self.K/self.P*t*np.cos(2*pi/self.P*t+self.phi_0) + self.f_GW*self.K/(2*pi)*(np.sin(2*pi/self.P*t+self.phi_0)-np.sin(self.phi_0))
return (self.K*self.f_GW*(self.P*(np.sin(self.phi_0)-np.sin(self.phi_0+2*pi*t/self.P))+2*np.pi*np.cos(self.phi_0+2*np.pi*t/self.P)*t))/(2*self.P*np.pi)
return self.f_GW*self.K*((t*np.cos((2*pi*t)/self.P + self.phi_0))/self.P - (np.cos((pi*t)/self.P + self.phi_0)*np.sin((pi*t)/self.P))/pi)
if i == 3:
return -self.f_GW/(2*pi)*self.P*self.K*(np.cos(2*pi*t/self.P+self.phi_0) - np.cos(self.phi_0))
return (self.f_GW*self.K*self.P*np.sin((pi*t)/self.P)*np.sin((pi*t)/self.P + self.phi_0))/pi
else:
print('Error in d_freq_int')
return self.strain(t)
def d_strain(self,t,i,bina=None,diff=1.): # 0: lnA, 1: f_0, 2: theta_S_Ec, 3: phi_S_Ec, 4: theta_L, 5: phi_L, 6: K, 7: P, 8: phi_0
def d_strain(self,t,i,bina=None,diff=1e-6): # 0: K, 1: P, 2: phi_0, 3: f_0, 4: ln(A), 5: theta_S_Ec, 6: phi_S_Ec, 7: theta_L, 8: phi_L
if bina is not None:
return (bina.strain(t)-self.strain(t))/diff
elif i >= 6:
i -= 5
return -sqth*self.A(t)*(2*pi*self.d_freq_int(t,i)+2*pi*self.d_freq(t,i)/c*AU*np.sin(self.theta_S_Ec)*np.cos(omega_E*t-self.phi_S_Ec))*np.sin(2*pi*self.freq_int(t)+self.polarization_phase(t)+self.doppler_phase(t))
elif i == 0:
if i <= 2:
return -sqth*self.A(t)*(2*pi*self.d_freq_int(t,i+1)+2*pi*self.d_freq(t,i+1)/c*AU*np.sin(self.theta_S_Ec)*np.cos(omega_E*t-self.phi_S_Ec))*np.sin(2*pi*self.freq_int(t)+self.polarization_phase(t)+self.doppler_phase(t))
if i == 3:
return -sqth*self.A(t)*(2*pi*self.d_freq_int(t,0)+2*pi*self.d_freq(t,0)/c*AU*np.sin(self.theta_S_Ec)*np.cos(omega_E*t-self.phi_S_Ec))*np.sin(2*pi*self.freq_int(t)+self.polarization_phase(t)+self.doppler_phase(t))
if i == 4:
return self.strain(t)
elif i == 1: # TODO
return -sqth*self.A(t)*(2*pi*self.T_obs+2*pi*self.d_freq(t,0)/c*AU*np.sin(self.theta_S_Ec)*np.cos(omega_E*t-self.phi_S_Ec))*np.sin(2*pi*self.freq_int(t)+self.polarization_phase(t)+self.doppler_phase(t))
else:
print('Error in d_strain')
raise ValueError
print('Error in d_strain')
raise ValueError
def d_strain2(self,i,bina=None,diff=1e-6): # 0: K, 1: P, 2: phi_0, 3: f_0, 4: ln(A), 5: theta_S_Ec, 6: phi_S_Ec, 7: theta_L, 8: phi_L
if bina is not None:
return lambda t: (bina.strain(t)-self.strain(t))/diff
if i <= 2:
return lambda t: -sqth*self.A(t)*(2*pi*self.d_freq_int(t,i+1)+2*pi*self.d_freq(t,i+1)/c*AU*np.sin(self.theta_S_Ec)*np.cos(omega_E*t-self.phi_S_Ec))*np.sin(2*pi*self.freq_int(t)+self.polarization_phase(t)+self.doppler_phase(t))
if i == 3:
return lambda t: -sqth*self.A(t)*(2*pi*self.d_freq_int(t,0)+2*pi*self.d_freq(t,0)/c*AU*np.sin(self.theta_S_Ec)*np.cos(omega_E*t-self.phi_S_Ec))*np.sin(2*pi*self.freq_int(t)+self.polarization_phase(t)+self.doppler_phase(t))
if i == 4:
return lambda t: self.strain(t)
print('Error in d_strain')
raise ValueError
def h_ixh_j(self,i,j,diff=1e-7):
def h_ixh_j(self,i,j,diff=1e-6):
funcA = self.h_i(i,diff)
funcB = self.h_i(j,diff)
return integrate.quad(lambda t: funcA(t)*funcB(t),0,self.T_obs,limit=int(self.T_obs*self.f_GW*20))[0]
def h_i(self,i,diff=1e-6):
B = None
C = None
if i == 2:
if i == 5:
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)
elif i == 3:
elif i == 6:
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)
elif i == 4:
elif i == 7:
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)
elif i == 5:
elif i == 8:
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)
if j == 2:
C = 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)
elif j == 3:
C = 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)
elif j == 4:
C = 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)
elif j == 5:
C = 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)
return integrate.quad(lambda t: self.d_strain(t,i,B,diff)*self.d_strain(t,j,C,diff),0,self.T_obs,limit=int(self.T_obs*self.f_GW*20),points=self.T_obs*np.linspace(0,1,int(self.T_obs*self.f_GW*2)))[0]
return self.d_strain2(i,B,diff)
def Fisher_mat(self,diff=1e-6):
mat = np.zeros((9,9),np.float64)
for key in [1,2]:
self.key = key
for i in np.arange(9): # 0: lnA, 1: f_0, 2: theta_S_Ec, 3: phi_S_Ec, 4: theta_L, 5: phi_L, 6: K, 7: P, 8: phi_0
for j in np.arange(i,9):
mat[i][j] += self.h_ixh_j(i, j,diff)
if j != i:
mat[j][i] = mat[i][j]
self.Fisher = 2.*mat/self.S_n
def Fisher_mat(self,mode=1,diff=1e-6):
if mode == 1:
mat = np.zeros((9,9),np.float64)
for key in [1,2]:
self.key = key
for i in np.arange(9): # 0: K, 1: P, 2: phi_0, 3: f_0, 4: ln(A), 5: theta_S_Ec, 6: phi_S_Ec, 7: theta_L, 8: phi_L
for j in np.arange(i,9):
mat[i][j] += self.h_ixh_j(i, j,diff)
if j != i:
mat[j][i] = mat[i][j]
else:
mat = np.zeros((3,3),np.float64)
for key in [1,2]:
self.key = key
for i in np.arange(3): # 0: K, 1: P, 2: phi_0, 3: f_0, 4: ln(A), 5: theta_S_Ec, 6: phi_S_Ec, 7: theta_L, 8: phi_L
for j in np.arange(i,3):
mat[i][j] += self.h_ixh_j(i, j,diff)
if j != i:
mat[j][i] = mat[i][j]
self.Fisher = mat #2*mat/S_n
return self.Fisher
def rel_uncertainty(self):
unc = np.sqrt(np.abs(np.diag(np.linalg.inv(self.Fisher_mat()))[6:])) #np.linalg.eigvals(self.Fisher_mat()) for diagonalized
return unc/np.array([self.K,self.P,1.])*self.signal_to_noise()
def rel_uncertainty(self,mode=1,diff=1e-6):
unc = np.sqrt(np.abs(np.diag(np.linalg.inv(self.Fisher_mat(mode,diff)))[:3])) #np.linalg.eigvals(self.Fisher_mat()) for diagonalized
return unc/np.array([self.K,self.P,1.])
def hxh(self):
res = 0.
for key in [1,2]:
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),points=self.T_obs*np.linspace(0,1,int(self.T_obs*self.f_GW*2)))[0]
res += integrate.quad(lambda t: self.strain(t)**2,0,self.T_obs,limit=int(self.T_obs*self.f_GW*20))[0]
return res
def signal_to_noise(self):
return np.sqrt(2*self.hxh()/self.S_n)
\ No newline at end of file
return np.sqrt(2*self.hxh()/self.S_n())
\ No newline at end of file
......@@ -30,37 +30,6 @@ 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)
......@@ -93,45 +62,20 @@ def sanity_plots():
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()
for i in np.arange(9):
plt.plot(hour,[C.h_i(t,i,B,1e-6) for t in hour])
plt.title(i)
plt.show()
def uncertainty_plot(n=100,mP=1):
Ps = np.logspace(-2,1,n)
def uncertainty_plot(n=100,mP=1,a=-2,b=1,mode=2):
Ps = np.logspace(a,b,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[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(mode)
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(mode)
uncs*=np.sqrt(B.S_n()/2)*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*=np.sqrt(B.S_n()/2)*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)
......@@ -188,8 +132,59 @@ def deriv_check():
plt.legend()
plt.show()
def pos_dep():
rel_unc = np.zeros((20,20,3))
for i,mu in enumerate(np.linspace(-1,1,20)):
for j,rad in enumerate(np.linspace(0,2*pi,20)):
rel_unc[j,i,:] = binary(np.arccos(mu),rad,1e3,np.arccos(-.2),4.,m1=.23,m2=.23,freq=1e-3,mP=5,P=1.5,theta_P=pi/2,phi_P=pi/2,T_obs=4).rel_uncertainty()
rel_unc = np.log10(rel_unc)
for k,label in enumerate(['K','P']):
plt.matshow(rel_unc[:,:,k])
plt.ylabel('$\mu$')
plt.xlabel('$\phi$')
plt.colorbar()
plt.title('Relative uncertainty in '+label)
plt.show()
return rel_unc
def test(n=100):
Ps = np.logspace(-2,1,n)
uncs = np.zeros((n,9),np.float64)
#uncs2 = np.zeros((n,3),np.float64)
for n,P0 in enumerate(Ps):
uncs[n] = np.diag(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).Fisher_mat())
#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/=np.diag(B.Fisher)
#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,uncs[:,3],'r-',label=r'3')
plt.legend()
plt.show()
plt.loglog(Ps,uncs[:,4],'b-',label=r'4')
plt.loglog(Ps,uncs[:,5],'g-',label=r'5')
plt.loglog(Ps,uncs[:,6],'r-',label=r'6')
plt.loglog(Ps,uncs[:,7],'b-',label=r'7')
plt.loglog(Ps,uncs[:,8],'g-',label=r'8')
plt.legend()
plt.show()
#plt.loglog(Ps,uncs2[:,0],'r--')
#plt.loglog(Ps,uncs2[:,1],'b--')
#plt.loglog(Ps,uncs2[:,2],'g--')
return uncs
#uncs = test()
#pos_dep()
#sanity_plots()
uncs, uncs2 = uncertainty_plot(15,mP=10)
uncs, uncs2 = uncertainty_plot(100,mP=10,a=-2,b=1,mode=2)
#A = np.power(uncs[:-2,:]*uncs[1:-1,:]*uncs[2:,:],1/3)[::3,:]
#one_year_degeneracy()
#deriv_check()
\ No newline at end of file
This diff is collapsed.
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