On Nov 2, 2020 LDAP authentication will be disabled - see http://status.ncsa.illinois.edu for details.

Commit 7db365dc authored by Brockton Brendal's avatar Brockton Brendal

removed eq_16 from power

parent 49893bbf
......@@ -51,7 +51,6 @@ import warnings
import scipy.optimize
import scipy.interpolate
import scipy.integrate
import matplotlib.pyplot as plt
#-----Function Definitions-----#
......@@ -375,229 +374,7 @@ def get_angular_momentum(python_strain):
if __name__ == "__main__":
#--------------------------------------------------------------------------
#Inserted equation 16 here
#--------------------------------------------------------------------------
def eq_16(psi4_path, modes_and_radius, path_for_strain, quasilocalmeasures_path):
#ar = np.loadtxt("C:\\Users\\Brock\\Documents\\UIUC\\Gravity Group\\POWER_project\\Extrapolated_Strain\\J0040_N40\\J0040_N40_strain_at_100.0_l2_m2.dat")
#ar = loadHDF5Series("simulations/J0040_N40/output-????/J0040_N40/mp_psi4.h5", "l2_m2_r100.0") ####Try this
ar = loadHDF5Series(psi4_path, modes_and_radius)
t = ar[:,0]
psi = ar[:,1]
impsi = ar[:,2]
plt.plot(t,psi)
plt.plot(t,impsi)
plt.show()
#%%
plt.show()
#%%
s_in = scipy.integrate.cumtrapz(psi,t)
ims_in = scipy.integrate.cumtrapz(impsi,t)
s_in = s_in - s_in[-1]
ims_in = ims_in - ims_in[-1]
plt.plot(t[1:],s_in)
plt.plot(t[1:],ims_in)
#%%
#%%
d_in = scipy.integrate.cumtrapz(s_in,t[1:])
d_in = d_in - d_in[-1]
imd_in = scipy.integrate.cumtrapz(ims_in,t[1:])
imd_in = imd_in - imd_in[-1]
plt.plot(t[2:],d_in)
plt.plot(t[2:],imd_in)
#%%
#A_val = np.loadtxt("C:\\Users\\Brock\\Documents\\UIUC\\Gravity Group\\POWER_project\\Gravitational_Waveform_Extractor\\POWER\\simulations\\J0040_N40\\output-0018\\J0040_N40\\quasilocalmeasures-qlm_scalars..asc")
A_val = loadHDF5Series(quasilocalmeasures_path)
r = float(167)
l = float(3)
m = float(2)
M = A_val[:,58][-1]
a = (A_val[:,37]/A_val[:,58])[-1]
#ORIGINAL: ar_a = np.loadtxt("C:\\Users\\Brock\\Documents\\UIUC\\Gravity Group\\POWER_project\\Extrapolated_Strain\\J0040_N40\\J0040_N40_strain_at_167.0_l3_m2.dat")
ar_a = np.loadtxt(path_for_strain)
t_a = ar_a[:,0]
psi_a = ar_a[:,1]
impsi_a = ar_a[:,2]
t_b = t_a
psi_b = np.zeros(len(psi_a))
impsi_b = np.zeros(len(impsi_a))
print (a,M)
#%%
A = 1-(2*M/r)
a_1 = r
a_2 = ((l-1)*(l+2))/(2*r)
a_3 = ((l-1)*(l+2)*(l**2 + l -4))/(8*r*r)
B = ((0+a*2j)/((l+1)**2))*((((l+3)*(l-1)*(l+m+1)*(l-m+1))/((2*l+1)*(2*l+3)))**(1/2))
b_1 = r
b_2 = l*(l+3)
C = ((0+a*2j)/((l)**2))*((((l+2)*(l-2)*(l+m)*(l-m))/((2*l-1)*(2*l+1)))**(1/2))
c_1 = r
c_2 = (l-2)*(l+1)
#%%
ans = A*(a_1*psi[2:] - a_2*s_in[1:] + a_3*d_in) + B*(b_1*np.gradient(psi_a, t_a)[2:] - b_2*psi_a[2:]) - C*(c_1*np.gradient(psi_b, t_b)[2:] - c_2*psi_b[2:])
imans = A*(a_1*impsi[2:] - a_2*ims_in[1:] + a_3*imd_in) + B*(b_1*np.gradient(impsi_a, t_a)[2:] - b_2*impsi_a[2:]) - C*(c_1*np.gradient(impsi_b, t_b)[2:] - c_2*impsi_b[2:])
#%%
f1 = scipy.integrate.cumtrapz(ans,t[2:])
f1 = f1-f1[-1]
imf1 = scipy.integrate.cumtrapz(imans,t[2:])
imf1 = imf1-imf1[-1]
f2 = scipy.integrate.cumtrapz(f1,t[3:])
f2 = f2-f2[-1]
imf2 = scipy.integrate.cumtrapz(imf1,t[3:])
imf2 = imf2-imf2[-1]
plt.plot(t[4:],f2)
f3_cmp = f2 + imf2*1j
imf3 = f3_cmp.imag
f3 = f3_cmp.real
#%%
#ORIGINAL: arr3 = np.loadtxt("C:\\Users\\Brock\\Documents\\UIUC\\Gravity Group\\POWER_project\\Extrapolated_Strain\\J0040_N40\\J0040_N40_strain_at_167.0_l2_m2.dat")
arr3 = np.loadtxt(path_for_strain)
rans = arr3[:,1]
rt = arr3[:,0]
rt.shape
#%%
plt.plot(t[4:],f2)
plt.plot(rt,rans)
#%%
rans_mod = scipy.interpolate.spline(np.copy(rt),np.copy(rans),np.copy(t))
#%%
plt.plot(t[4:],f2)
plt.plot(t,rans_mod)
#%%
#ORIGINAL: pwr = np.loadtxt("C:\\Users\\Brock\\Documents\\UIUC\\Gravity Group\\POWER_project\\Extrapolated_Strain\\J0040_N40\\J0040_N40_strain_at_167.0_l2_m2.dat")
pwr = np.loadtxt(path_for_strain)
tp = pwr[:,0]
rep = pwr[:,1]
imp = pwr[:,2]
plt.plot(tp,rep,label="Re{h(t)}:POWER")
plt.plot(t[4:],f2,label="Re{h(t)}:Nakano")
plt.legend()
plt.xlabel("time [100 solar masses]")
plt.ylabel("Re{h(t)} = strain [dimensionless]")
plt.show()
tp1 = tp - tp[0]
#%%
Ap = (rep**2 + imp**2)**(1/2.0)
Ae = (f3**2 + imf3**2)**(1/2.0)
Ap_m = scipy.interpolate.spline(np.copy(tp1),np.copy(Ap),np.copy(t))
Ae_m = scipy.interpolate.spline(np.copy(t[4:]),np.copy(Ae),np.copy(t))
plt.plot(t,Ap_m)
plt.plot(t,Ae_m)
#%%
print (f3, imf3)
#%%
argp = np.arctan2(imp, rep)
arge = np.arctan2(imf3, f3)
argp_m = scipy.interpolate.spline(np.copy(tp1),np.copy(argp),np.copy(t))
arge_m = scipy.interpolate.spline(np.copy(t[4:]),np.copy(arge),np.copy(t))
plt.plot(t,argp_m, label="Phase(POWER)")
plt.plot(t,arge_m, label="Phase(Nakano)")
plt.legend()
plt.xlabel("time [100 solar masses]")
plt.ylabel("phase [rad]")
plt.show()
#%%
m = max(Ap_m)
max_ap = 0.0
for i in range(len(t)):
if Ap_m[i]==m:
max_ap = i
break
tmax_ap = t[max_ap]
m1 = max(Ae_m)
max_ae = 0.0
for i in range(len(t)):
if Ae_m[i]==m1:
max_ae = i
break
tmax_ae = t[max_ae]
print (tmax_ap, tmax_ae)
#%%
tp2 = t-tmax_ap
te2 = t-tmax_ae
plt.plot(tp2,Ap_m, label="Amplitude(POWER)")
plt.plot(te2,Ae_m, label="Amplitude(Nakano)")
plt.legend()
plt.xlabel("time [100 solar masses]")
plt.ylabel("amplitude = strain [dimentionless]")
plt.show()
#%%
phitp = np.unwrap(argp_m)
phite = np.unwrap(arge_m)
plt.plot(t,phitp)
plt.plot(t,phite)
#%%
phihp = phitp - phitp[max_ap]
phihe = phite - phite[max_ae]
plt.plot(t,phihp, label="Phase(POWER)")
plt.plot(t,phihe, label="Phase(Nakano)")
plt.legend()
plt.xlabel("time [100 solar masses]")
plt.ylabel("phase [rad]")
plt.show()
#%%
diffx = phihp-phihe
plt.plot(t,diffx-0.25)
plt.ylim(-1.,1.)
plt.xlabel("time [100 solar masses]")
plt.ylabel("phase [rad]")
plt.show()
#--------------------------------------------------------------------------
#End of equation 16
#--------------------------------------------------------------------------
#Initialize simulation data
if(len(sys.argv) < 2):
print("Pass in the number n of the n innermost detector radii to be used in the extrapolation (optional, default=all) and the simulation folders (e.g., ./power.py 6 ./simulations/J0040_N40 /path/to/my_simulation_folder).")
......
Markdown is supported
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