-
Notifications
You must be signed in to change notification settings - Fork 2
/
plot_kern_tuneer.py
143 lines (115 loc) · 4.86 KB
/
plot_kern_tuneer.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
import timing
import numpy as np
import matplotlib.pyplot as plt
import functions as fn
from os import getcwd
import get_kernels as gkerns
kernclock = timing.stopclock()
tstamp = kernclock.lap
r = np.loadtxt('r.dat')
r_start = 0.68
r_end = 1.0
start_ind = fn.nearest_index(r,r_start)
end_ind = fn.nearest_index(r,r_end)
r = r[start_ind:end_ind+1]
OM = np.loadtxt('OM.dat')
# n1,l1 = 4,3
# n2,l2 = 1,10
# n3,l3 = 0,60
n1,l1 = 15,102
n2,l2 = 15,50
n3,l3 = 15,10
omega_list = np.loadtxt('muhz.dat')
omega_nl1 = omega_list[fn.find_nl(n1, l1)]
omega_nl2 = omega_list[fn.find_nl(n2, l2)]
omega_nl3 = omega_list[fn.find_nl(n3, l3)]
m = np.array([0])
m_ = np.array([0])
s = np.array([2])
#condition about whether or not to scale by rho
multiplyrho = True
smoothen = True
# plot_fac = OM**2 * 1e12 * (4.*np.pi/3) * 1e-10 #unit = muHz G^(-2) V_sol^(-1)
plot_fac = OM**2 * 1e12 * 1e-10 #unit = muHz G^(-2)
#extracting rho in an unclean fashion
rho,__,__,__,__,__,__ = np.array(gkerns.Hkernels(n1,l1,m,n1,l1,m,s,r)\
.ret_kerns_axis_symm(smoothen,a_coeffkerns = True))
#Kernels for a-coefficients for Lorentz stress
# kern = gkerns.Hkernels(n1,l1,m,n1,l1,m,s,r)
__, Bmm1, B0m1,B001, Bpm1,_,_ = np.array(gkerns.Hkernels(n1,l1,m,n1,l1,m,s,r)\
.ret_kerns_axis_symm(smoothen,a_coeffkerns = True))*plot_fac/(-2*omega_nl1)
__, Bmm2, B0m2,B002, Bpm2,_,_ = np.array(gkerns.Hkernels(n2,l2,m,n2,l2,m,s,r).\
ret_kerns_axis_symm(smoothen,a_coeffkerns = True))*plot_fac/(-2*omega_nl2)
__, Bmm3, B0m3,B003, Bpm3,_,_ = np.array(gkerns.Hkernels(n3,l3,m,n3,l3,m,s,r).\
ret_kerns_axis_symm(smoothen,a_coeffkerns = True))*plot_fac/(-2*omega_nl3)
#############################################################################
#Purely lorentz stress kernels
# Bmm1, B0m1,B001, Bpm1,_,_ = np.array(gkerns.Hkernels(n1,l1,m,n1,l1,m,s,r)\
# .ret_kerns_axis_symm())*plot_fac
# Bmm2, B0m2,B002, Bpm2,_,_ = np.array(gkerns.Hkernels(n2,l2,m,n2,l2,m,s,r).\
# ret_kerns_axis_symm())*plot_fac
# Bmm3, B0m3,B003, Bpm3,_,_ = np.array(gkerns.Hkernels(n3,l3,m,n3,l3,m,s,r).\
# ret_kerns_axis_symm())*plot_fac
tstamp('kernel calculation time')
npts = 300 #check the npts in get_kernels
r_new = np.linspace(np.amin(r),np.amax(r),npts)
r = r_new
# rho_i = np.loadtxt('rho.dat')
# rho_i = rho_i[start_ind:end_ind+1]
# rho = np.linspace(rho_i[0],rho_i[-1],npts)
if(multiplyrho==True):
Bmm1, B0m1,B001, Bpm1 = rho*Bmm1, rho*B0m1,rho*B001, rho*Bpm1
Bmm2, B0m2,B002, Bpm2 = rho*Bmm2, rho*B0m2,rho*B002, rho*Bpm2
Bmm3, B0m3,B003, Bpm3 = rho*Bmm3, rho*B0m3,rho*B003, rho*Bpm3
#plot_fac should have units of muHz G^(-2) M_sol V_sol^(-1)
#M_sol/(4pi/3 * R_sol = 1 g/cc)
dpi = 80
plt.figure(num=None, figsize=(11, 8), dpi=dpi, facecolor='w', edgecolor='k')
plt.subplot(221)
plt.plot(r,Bmm1[0,0], label = '('+str(n1)+','+str(l1)+')')
plt.plot(r,Bmm2[0,0], label = '('+str(n2)+','+str(l2)+')')
plt.plot(r,Bmm3[0,0], label = '('+str(n3)+','+str(l3)+')')
plt.xlabel('$r/R_{\odot}$',fontsize=14)
plt.xlim(r_start,r_end)
plt.legend()
plt.grid(True)
# plt.ylabel('$\mathcal{B}^{--}_{20}$ in $(\mu Hz^2 G^{-2}V_{\odot}^{-1})$')
plt.ylabel('$\\rho\mathcal{A}^{--}_{20}$ in $(\mu Hz G^{-2}M_{\odot}V_{\odot}^{-1})$',\
fontsize = 14)
plt.subplot(222)
plt.plot(r,B0m1[0,0], label = '('+str(n1)+','+str(l1)+')')
plt.plot(r,B0m2[0,0], label = '('+str(n2)+','+str(l2)+')')
plt.plot(r,B0m3[0,0], label = '('+str(n3)+','+str(l3)+')')
plt.xlabel('$r/R_{\odot}$',fontsize=14)
plt.xlim(r_start,r_end)
plt.legend()
plt.grid(True)
# plt.ylabel('$\mathcal{B}^{0-}_{20}$ in $(\mu Hz^2 G^{-2}V_{\odot}^{-1})$')
plt.ylabel('$\\rho\mathcal{A}^{0-}_{20}$ in $(\mu Hz G^{-2}M_{\odot}V_{\odot}^{-1})$',\
fontsize = 14)
plt.subplot(223)
plt.plot(r,B001[0,0], label = '('+str(n1)+','+str(l1)+')')
plt.plot(r,B002[0,0], label = '('+str(n2)+','+str(l2)+')')
plt.plot(r,B003[0,0], label = '('+str(n3)+','+str(l3)+')')
plt.xlabel('$r/R_{\odot}$',fontsize=14)
plt.xlim(r_start,r_end)
plt.legend()
plt.grid(True)
# plt.ylabel('$\mathcal{B}^{00}_{20}$ in $(\mu Hz^2 G^{-2}V_{\odot}^{-1})$')
plt.ylabel('$\\rho\mathcal{A}^{00}_{20}$ in $(\mu Hz G^{-2}M_{\odot}V_{\odot}^{-1})$',\
fontsize = 14)
plt.subplot(224)
plt.plot(r,Bpm1[0,0], label = '('+str(n1)+','+str(l1)+')')
plt.plot(r,Bpm2[0,0], label = '('+str(n2)+','+str(l2)+')')
plt.plot(r,Bpm3[0,0], label = '('+str(n3)+','+str(l3)+')')
plt.xlabel('$r/R_{\odot}$',fontsize=14)
plt.xlim(r_start,r_end)
plt.legend()
plt.grid(True)
# plt.ylabel('$\mathcal{B}^{+-}_{20}$ in $(\mu Hz^2 G^{-2}V_{\odot}^{-1})$')
plt.ylabel('$\\rho\mathcal{A}^{+-}_{20}$ in $(\mu Hz G^{-2}M_{\odot}V_{\odot}^{-1})$',\
fontsize = 14)
plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=0.95, hspace=0.25,
wspace=0.35)
plt.savefig('./kern2.pdf')
plt.show()