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