-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathplot_kern_shravan.py
160 lines (124 loc) · 4.8 KB
/
plot_kern_shravan.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
#PLOTTING CORRECTION TO SHRAVAN'S KERNEL
import timing
kernclock = timing.stopclock() #stopclock object for timing program
tstamp = kernclock.lap
import numpy as np
import matplotlib.pyplot as plt
import functions as fn
from os import getcwd
tstamp('library loading') #printing elapsed time from beginning of runtime
n,l,m = 1,60,0
n_,l_,m_ = n,l,2
nl = fn.find_nl(n,l)
nl_ = fn.find_nl(n_,l_)
s = 22
t = m_-m
#Savitsky golay filter for smoothening
window = 45 #must be odd
order = 3
if(nl == None or nl_ == None):
print("Mode not found. Exiting."); exit()
#loading required functions
eig_dir = (getcwd() + '/eig_files')
U,V = fn.load_eig(n,l,eig_dir)
U_,V_= fn.load_eig(n_,l_,eig_dir)
r = np.loadtxt('r.dat')
rho = np.loadtxt('rho.dat')
#interpolation params
npts = 30000
r_new = np.linspace(np.amin(r),np.amax(r),npts)
tstamp('files loading')
#setting up shorthand repeatedly used in kernel evaluation
def wig_red(m1,m2,m3):
'''3j symbol with upper row fixed'''
return fn.wig(l_,s,l,m1,m2,m3)
om = fn.omega
p = (-1)**(l+l_+s) #parity of selected modes
prefac = np.sqrt((2*l_+1.) * (2*s+1.) * (2*l+1.) / (4.* np.pi)) * wig_red(-m_,t,m)
#EIGENFUCNTION DERIVATIVES
#smoothing
U,dU,d2U = fn.smooth(U,r,window,order,npts)
V,dV,d2V = fn.smooth(V,r,window,order,npts)
U_,dU_,d2U_ = fn.smooth(U_,r,window,order,npts)
V_,dV_,d2V_ = fn.smooth(V_,r,window,order,npts)
rho_sm, __, __ = fn.smooth(rho,r,window,order,npts)
#re-assigning with smoothened variables
r = r_new
rho = rho_sm
##no smoothing
#dU, dV = np.gradient(U,r), np.gradient(V,r)
#dU_, dV_ = np.gradient(U_,r), np.gradient(V_,r)
#d2U_,d2V = 0.,0.
#B-- EXPRESSION
Bmm = -r*(wig_red(0,-2,2)*om(l,0)*om(l,2)*V*dU_ + wig_red(2,-2,0)*om(l_,0)* \
om(l_,2)*V_*dU)
Bmm += wig_red(1,-2,1)*om(l_,0)*om(l,0)*(U-V)*(U_ - V_ + r*dV_)
Bmm *= ((-1)**m_)*prefac/(r**2)
#B-- EXTRA
Bmm_ = om(l_,0)*(wig_red(2,-2,0)*om(l_,2)*U*(V_ - r*dV_) + om(l,0)*V \
*(wig_red(3,-2,-1)*om(l_,2)*om(l_,3)*V_ + wig_red(1,-2,1) \
*(-U_ + V_ + om(l_,2)**2 *V_ - r*dV_)))
Bmm_ *= (-1)**(1+m_) *prefac/r**2
#B0- EXPRESSION
B0m = wig_red(1,-1,0)*om(l_,0)*(U - (om(l,0)**2)*V)*(U_ - V_ + r*dV_)
B0m += om(l,0)*(om(l_,0)*(wig_red(-1,-1,2)*om(l,2)*V*(U_ - V_ + r*dV_) \
+ 2*r*wig_red(2,-1,-1)*om(l_,2)*V_*dV) + wig_red(0,-1,1) \
*((U-V)*(2*U_ - 2*(om(l_,0)**2)*V_ - r*dU_) + r**2 * dU_*dV))
B0m *= 0.5*((-1)**m_)*prefac/r**2
#B0- EXTRA
B0m_ = om(l,0)*V*(wig_red(2,-1,-1)*om(l_,0)*om(l_,2)*(U_ - 3*V_ + r*dV_) \
+ wig_red(0,-1,1)*((2+om(l_,0)**2)*U_ - 2*r*dU_ + om(l_,0)**2 \
*(-3*V_ + r*dV_)))
B0m_ += wig_red(1,-1,0)*om(l_,0)*U*(U_ - V_ - r*(dU_ - dV_ + r*d2V_))
B0m_ *= 0.5*((-1)**m_)*prefac/r**2
#B00 OLD
B00 = -wig_red(0,0,0)*(2*U_ - 2*om(l_,0)**2 * V_ - r*dU_)*(-2*U + 2*om(l,0)**2 *V + \
r*dU)
B00 -= 2*r*(wig_red(-1,0,1) + wig_red(1,0,-1))*om(l_,0)*om(l,0) \
*(U_ - V_ + r*dV_)*dV
B00 *= 0.5*((-1)**m_)*prefac/r**2
#B00 EXTRA
B00_ = -(wig_red(-1,0,1) + wig_red(1,0,-1)) * om(l_,0)*om(l,0) * V*(-4*U_+2*(1+om(l_,0)**2)*V_+r*(dU_-2*dV_))
B00_ += wig_red(0,0,0)*U*(2*U_-2*r*dU_-2*om(l_,0)**2 *(V_-r*dV_)+r*r*d2U_)
B00_ *= 0.5*((-1)**m_)*prefac/r**2
#B+- OLD
Bpm = -r**2 * wig_red(0,0,0)*dU_*dU
Bpm += om(l_,0)*om(l,0)*(-2*(wig_red(-2,0,2)+wig_red(2,0,-2))*om(l_,2)*om(l,2)*V_*V \
+ wig_red(-1,0,1)*(U-V)*(U_ - V_ + r*dV_) + wig_red(1,0,-1) \
*(U-V)*(U_ - V_ + r*dV_))
Bpm *= 0.5*((-1)**m_)*prefac/r**2
#B0+- EXTRA
Bpm_ = (wig_red(-1,0,1) + wig_red(1,0,-1)) * om(l_,0)*om(l,0) * V * (U_-V_+r*(-dU_+dV_))
Bpm_ += wig_red(0,0,0) * r*r*U*d2U_
Bpm_ *= 0.5*((-1)**m_)*prefac/r**2
tstamp('calculations')
r_start = 0.9
start_ind = fn.nearest_index(r,r_start)
plt.subplot(2,2,1)
plt.plot(r[start_ind:],(2*rho*Bpm)[start_ind:],'k-.')
plt.plot(r[start_ind:],(2*rho*Bpm_)[start_ind:],'r--')
plt.plot(r[start_ind:],(2*rho*(Bpm+Bpm_))[start_ind:],'b-')
plt.title('$\mathcal{B}^{+-}$')
plt.subplot(2,2,2)
plt.plot(r[start_ind:],(2*rho*Bmm)[start_ind:],'k-.')
plt.plot(r[start_ind:],(2*rho*Bmm_)[start_ind:],'r--')
plt.plot(r[start_ind:],(2*rho*(Bmm+Bmm_))[start_ind:],'b-')
plt.title('$\mathcal{B}^{--}$')
plt.subplot(2,2,3)
plt.plot(r[start_ind:],(2*rho*B0m)[start_ind:],'k-.')
plt.plot(r[start_ind:],(2*rho*B0m_)[start_ind:],'r--')
plt.plot(r[start_ind:],(2*rho*(B0m+B0m_))[start_ind:],'b-')
plt.title('$\mathcal{B}^{0-}$')
plt.subplot(2,2,4)
plt.plot(r[start_ind:],(2*rho*B00)[start_ind:],'k-.')
plt.plot(r[start_ind:],(2*rho*B00_)[start_ind:],'r--')
plt.plot(r[start_ind:],(2*rho*(B00+B00_))[start_ind:],'b-')
plt.title('$\mathcal{B}^{00}$')
#plt.plot(r[start_ind:],(2*rho*Bpm)[start_ind:],'k-',label = '$\mathcal{B}^{+-}$')
##plt.plot(r[start_ind:],(rho*Bmm)[start_ind:],'r-',label = '$\mathcal{B}^{--}$')
##plt.plot(r[start_ind:],(2*rho*B0m)[start_ind:],'k-',label = '$\mathcal{B}^{0-}$')
#plt.plot(r[start_ind:],(-rho*B00)[start_ind:],'k--',label = '$\mathcal{B}^{00}$')
#plt.grid(True)
plt.legend()
plt.show()
tstamp('plotting')