1 tplarson 1.1 import matplotlib.pyplot as plt
2 from matplotlib import animation
3 import numpy as np
4 import sosh
5 import sys
6 import os
7
8 nargs=len(sys.argv)
9 #print(sys.argv)
10 if (nargs > 1):
11 for i in range(1, nargs):
12 exec(sys.argv[i])
13
14 try:
15 freqscale
16 except:
|
17 tplarson 1.2 freqscale=1.0/3
|
18 tplarson 1.1
19 try:
20 dpi
21 except:
22 dpi = 300
23
24 try:
25 nframes
26 except:
27 nframes = 64
28
29 try:
|
30 tplarson 1.2 colorshift
|
31 tplarson 1.1 except:
|
32 tplarson 1.4 colorshift=1
|
33 tplarson 1.2
34 try:
35 animate
36 except:
37 animate=1
38
39 try:
40 show
41 except:
42 show=1
|
43 tplarson 1.1
44 try:
45 rsurf
46 except:
47 rsurf=1.0
48
49 try:
50 pixels
51 except:
52 pixels=1000
53
54 try:
|
55 tplarson 1.5 distobs
56 except:
57 distobs=220.0
58
59 try:
60 xoffset
61 except:
62 xoffset=0.0
63
64 try:
65 yoffset
66 except:
67 yoffset=0.0
68
69 try:
70 pangle
71 except:
72 pangle=0.0
73
74 try:
|
75 tplarson 1.1 figsize
76 except:
77 figsize=5
78
|
79 tplarson 1.2 sosh.nframes=nframes
80 sosh.dpi=dpi
81 sosh.icolshift=colorshift
82
|
83 tplarson 1.1 sosh.loadmodel()
84
|
85 tplarson 1.5 (r, theta, phi) = sosh.image2rtheta(xpixels=pixels,ypixels=pixels,distobs=distobs,pangle=pangle,xoffset=xoffset,yoffset=yoffset)
|
86 tplarson 1.1 x=np.cos(theta)
87 (nx,ny)=x.shape
88 if (rsurf < 1.0):
89 newind = (r <= rsurf)
|
90 tplarson 1.5 r.mask = np.logical_not(newind)
91 x.mask = np.logical_not(newind)
|
92 tplarson 1.1
|
93 tplarson 1.2 #modeparms=np.loadtxt('mdi.average.modes')
|
94 tplarson 1.1 # use model values instead of measured values
|
95 tplarson 1.2 #modeparms=np.loadtxt('model.surface.modes')
96 #lmod=modeparms[:,0]
97 #nmod=modeparms[:,1]
|
98 tplarson 1.1
99 arrlist={}
100 lmaxlist={}
101 varlist={}
102
103 rlist=[]
104 tlist=[]
105 plist=[]
106 freqlist=[]
107 amplist=[]
108 rclist=[]
109 hclist=[]
110
111 interval = 1.0/25
|
112 tplarson 1.2 maxsave=0.0
|
113 tplarson 1.1
114 #Animation function to call
|
115 tplarson 1.2 def calcsumt(i):
116
|
117 tplarson 1.1 sumr=0.0
118 sumt=0.0
119 sump=0.0
|
120 tplarson 1.2 plotvar=sosh.plotvar
121
122 if plotvar in ['Vr','Vmag','Vsq']:
123 for k in range(nmodes):
124 sumr += freqlist[k]*rclist[k]*rlist[k]*np.exp(-1.0j*2*np.pi*freqlist[k]*freqscale*float(i) / np.float(nframes))
125 if plotvar in ['Vt', 'Vh', 'Vmag','Vsq']:
126 for k in range(nmodes):
127 sumt += freqlist[k]*hclist[k]*tlist[k]*np.exp(-1.0j*2*np.pi*freqlist[k]*freqscale*float(i) / np.float(nframes))
128 if plotvar in ['Vt', 'Vh', 'Vmag','Vsq']:
129 for k in range(nmodes):
130 sump += freqlist[k]*hclist[k]*plist[k]*np.exp(-1.0j*2*np.pi*freqlist[k]*freqscale*float(i) / np.float(nframes))
131
132 if plotvar in ['Vh', 'Vmag', 'Vsq']:
133 vh2=np.square(sumt.real) + np.square(sump.real)
134 if plotvar in ['Vmag', 'Vsq']:
135 v2=np.square(sumr.real) + vh2
|
136 tplarson 1.1
137 if (plotvar == 'Vr'):
138 d=sumr.real
139 if (plotvar == 'Vt'):
140 d=sumt.real
141 if (plotvar == 'Vp'):
142 d=sump.real
143 if (plotvar == 'Vh'):
144 d=np.sqrt(vh2)
145 if (plotvar == 'Vmag'):
146 d=np.sqrt(v2)
147 if (plotvar == 'Vsq'):
148 d=v2
149
|
150 tplarson 1.2 return d
151
152 def sumanimate(i):
153
154 d=calcsumt(i)
155 global maxsave
156 mn,mx=d.min(),d.max()
157 maxabs=np.abs([mn,mx]).max()
158 if (maxabs > maxsave):
159 maxsave=maxabs
160 # print("new max=%f"%maxabs, end=" ", flush=True)
161
|
162 tplarson 1.1 im.set_data(d)
163 fig.canvas.draw()
|
164 tplarson 1.2 return maxabs
|
165 tplarson 1.1
166 sosh.animate=sumanimate
|
167 tplarson 1.2 sosh.calcimage=calcsumt
|
168 tplarson 1.1 sosh.isave=0
169
|
170 tplarson 1.2
|
171 tplarson 1.1 print("You may enter 'q' to quit or 'c' to change parameters.")
172 instr=sosh.catchc("Enter number of modes: ",2)
173 if (instr == 'q'):
174 sys.exit()
175 nmodes=int(instr)
|
176 tplarson 1.4 for i in range(abs(nmodes)):
|
177 tplarson 1.1 print("MODE #%i" % (i+1))
178
|
179 tplarson 1.4 (l,m,n)=sosh.querylmn(nmodes)
|
180 tplarson 1.1 if l == -1:
181 sys.exit()
|
182 tplarson 1.2 ntest = sosh.modeln[l==sosh.modell]
|
183 tplarson 1.1 while (n not in ntest):
184 print("That n was not modelled for that l. Modelled values are in the range %i to %i. Try again." \
185 % (ntest.min(),ntest.max()))
|
186 tplarson 1.4 (l,m,n)=sosh.querylmn(nmodes)
|
187 tplarson 1.2 ntest = sosh.modeln[l==sosh.modell]
|
188 tplarson 1.1
189 signedm=m
190 m=np.abs(m)
191 if m not in lmaxlist.keys() or l > lmaxlist[m]:
192 lmaxlist[m]=l
193 plm=np.ma.array(np.zeros((nx,ny,l-m+1)))
194 dplm=np.ma.array(np.zeros((nx,ny,l-m+1)))
195 (y,dy)=sosh.setplm(l,m,x,plm,dplm)
196 arrlist[m]=(plm,dplm)
197 print("Spherical harmonic calculated and saved.")
198 else:
199 (plm,dplm)=arrlist[m]
200 (y,dy)=(plm[...,l-m],dplm[...,l-m])
201 print("Spherical harmonic retrieved.")
202
|
203 tplarson 1.3 if (i == 0):
204 if (m > 0):
|
205 tplarson 1.5 phi+=(np.pi/4)/m
206 # phi[:,int(nx/2):nx]=(np.pi/4)/m
207 # phi[:,0:int(nx/2)]=(np.pi/4)/m + np.pi
208 # else:
209 # phi=0.0*theta
|
210 tplarson 1.1 ylm=y*np.exp(1.0j*signedm*phi)
211 dylmt=-np.sin(theta)*dy*np.exp(1.0j*signedm*phi)
|
212 tplarson 1.2 dylmp=1.0j*signedm*y*np.exp(1.0j*signedm*phi)
|
213 tplarson 1.1
214 rlist.append(ylm)
215 tlist.append(dylmt)
216 plist.append(dylmp/np.sin(theta))
217
218 xir, xih = sosh.getradial(l,n)
219 fr=np.interp(r,sosh.rmesh,xir*np.sqrt(sosh.rho))
220 fh=np.interp(r,sosh.rmesh,xih*np.sqrt(sosh.rho))
221
|
222 tplarson 1.2 idx = ((sosh.modeln == n) & (sosh.modell == l))
|
223 tplarson 1.3 freq = sosh.modelnu[idx]*1000.0
|
224 tplarson 1.2 # ind = ((l==lmod) & (n==nmod))
225 # freq = float(modeparms[ind,2])/1000.0
|
226 tplarson 1.1 rclist.append(fr)
227 hclist.append(fh)
228 freqlist.append(freq)
|
229 tplarson 1.2 print("Scaled frequency = %f." % (freq*freqscale))
|
230 tplarson 1.1
231 # if using measured mode parameters, this would
232 # allow you to compute frequency as a function of m
233 # from the fitted a-coefficients
234 # ai=np.append([0.0],modeparms[ind,12:18]/1000)
235 # pols=sosh.apols(l,6)
236 # fx=np.matmul(pols,ai)
237 # freq=modeparms[ind,2]+fx
238 # freqlist.append(freq[l+signedm]*freqscale)
239 # amplist.append(modeparms[ind,3])
240
|
241 tplarson 1.4 nmodes=abs(nmodes)
|
242 tplarson 1.1 sumr=0.0
243 sumt=0.0
244 sump=0.0
245 for i in range(nmodes):
246 sumr += rclist[i]*rlist[i]*freqlist[i]
247 sumt += hclist[i]*tlist[i]*freqlist[i]
248 sump += hclist[i]*plist[i]*freqlist[i]
249
250 vh2=(np.square(sumt.real) + np.square(sump.real))
251 v2=(np.square(sumr.real) + vh2)
252
253 varlist['Vr']=sumr.real
254 varlist['Vt']=sumt.real
255 varlist['Vp']=sump.real
256 varlist['Vh']=np.sqrt(vh2)
257 varlist['Vmag']=np.sqrt(v2)
258 varlist['Vsq']=v2
259
260 while True:
261
262 var=varlist[sosh.plotvar]
|
263 tplarson 1.2 maxsave=0.0
|
264 tplarson 1.1
265 fig,im = sosh.drawfigure(var, fsize=figsize)
266 if (sosh.isave != 0):
|
267 tplarson 1.4 sosh.savefigure(ianimate=animate, label='_interiorsum')
|
268 tplarson 1.2 if (show != 0):
269 if (animate != 0):
270 anim = animation.FuncAnimation(fig, sumanimate, frames=nframes, interval=interval)
271 plt.show()
|
272 tplarson 1.1
273 print()
274 c=sosh.queryplotparms()
275 if (c == 'q'):
276 break
|
277 tplarson 1.4 fs=sosh.catchc("Enter new freqscale: ", freqscale)
278 if (fs == 'q'):
279 break
280 else:
281 freqscale=float(fs)
282 nf=sosh.catchc("Enter new nframes: ", nframes)
283 if (nf == 'q'):
284 break
285 else:
286 nframes=int(nf)
287 sosh.nframes=nframes
|
288 tplarson 1.1
289
290
|
Return to
addradial.py
CVS log
Up to [Development] / JSOC / proj / globalhs / sosh