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                freqscale=1.0/3000
 18              
 19              try:
 20                dpi
 21              except:
 22 tplarson 1.1   dpi = 300
 23              
 24              try:
 25                nframes
 26              except:
 27                nframes = 64
 28              
 29              try:
 30                animateflag
 31              except:
 32                animateflag=1
 33              
 34              try:
 35                rsurf
 36              except:
 37                rsurf=1.0
 38              
 39              try:
 40                pixels
 41              except:
 42                pixels=1000
 43 tplarson 1.1 
 44              try:
 45                bangle
 46              except:
 47                bangle=-30
 48              
 49              try:
 50                figsize
 51              except:
 52                figsize=5
 53              
 54              # not needed if using table of surface values
 55              sosh.loadmodel()
 56              
 57              #(phi,theta)=sosh.image2sphere(xpixels=pixels,ypixels=pixels,bangle=bangle)
 58              (r, theta) = sosh.image2rtheta(xpixels=pixels,ypixels=pixels)
 59              x=np.cos(theta)
 60              (nx,ny)=x.shape
 61              if (rsurf < 1.0):
 62                newind = (r <= rsurf)
 63                r.mask = np.logical_not(newindex)
 64 tplarson 1.1 
 65              #modeparms=np.loadtxt('modes.average')
 66              # use model values instead of measured values
 67              modeparms=np.loadtxt('modes.model.txt')
 68              lmod=modeparms[:,0]
 69              nmod=modeparms[:,1]
 70              
 71              arrlist={}
 72              lmaxlist={}
 73              varlist={}
 74              
 75              rlist=[]
 76              tlist=[]
 77              plist=[]
 78              freqlist=[]
 79              amplist=[]
 80              rclist=[]
 81              hclist=[]
 82              
 83              interval = 1.0/25
 84              
 85 tplarson 1.1 #Animation function to call
 86              def sumanimate(i):
 87                sumr=0.0
 88                sumt=0.0
 89                sump=0.0
 90                for k in range(nmodes):
 91                  sumr += freqlist[k]*rclist[k]*rlist[k]*np.exp(-1.0j*(2*np.pi*freqlist[k]*float(i) / np.float(nframes)))
 92                  sumt += freqlist[k]*hclist[k]*tlist[k]*np.exp(-1.0j*(2*np.pi*freqlist[k]*float(i) / np.float(nframes)))
 93                  sump += freqlist[k]*hclist[k]*plist[k]*np.exp(-1.0j*(2*np.pi*freqlist[k]*float(i) / np.float(nframes)))
 94                vh2=(np.square(sumt.real) + np.square(sump.real))
 95                v2=(np.square(sumr.real) + vh2)
 96              
 97                plotvar=sosh.plotvar
 98                if (plotvar == 'Vr'):
 99                  d=sumr.real
100                if (plotvar == 'Vt'):
101                  d=sumt.real
102                if (plotvar == 'Vp'):
103                  d=sump.real
104                if (plotvar == 'Vh'):
105                  d=np.sqrt(vh2)
106 tplarson 1.1   if (plotvar == 'Vmag'):
107                  d=np.sqrt(v2)
108                if (plotvar == 'Vsq'):
109                  d=v2
110              
111                im.set_data(d)
112                fig.canvas.draw()
113                return
114              
115              sosh.animate=sumanimate
116              sosh.isave=0
117              
118              print("You may enter 'q' to quit or 'c' to change parameters.")
119              instr=sosh.catchc("Enter number of modes: ",2)
120              if (instr == 'q'):
121                sys.exit()
122              nmodes=int(instr)
123              for i in range(nmodes):
124                print("MODE #%i" % (i+1))
125              
126                (l,m,n)=sosh.querylmn()
127 tplarson 1.1   if l == -1:
128                  sys.exit()
129                ntest = nmod[l==lmod]
130                while (n not in ntest):
131                  print("That n was not modelled for that l. Modelled values are in the range %i to %i. Try again." \
132                        % (ntest.min(),ntest.max()))
133                  (l,m,n)=sosh.querylmn()
134                  ntest = nmod[l==lmod]
135              
136                signedm=m
137                m=np.abs(m)
138                if m not in lmaxlist.keys() or l > lmaxlist[m]:
139                  lmaxlist[m]=l
140                  plm=np.ma.array(np.zeros((nx,ny,l-m+1)))
141                  dplm=np.ma.array(np.zeros((nx,ny,l-m+1)))
142                  (y,dy)=sosh.setplm(l,m,x,plm,dplm)
143                  arrlist[m]=(plm,dplm)
144                  print("Spherical harmonic calculated and saved.")
145                else:
146                  (plm,dplm)=arrlist[m]
147                  (y,dy)=(plm[...,l-m],dplm[...,l-m])
148 tplarson 1.1     print("Spherical harmonic retrieved.")
149              
150                if (m > 0):
151                  phi=(np.pi/4)/m
152                else:
153                  phi=0.0
154                ylm=y*np.exp(1.0j*signedm*phi)
155                dylmt=-np.sin(theta)*dy*np.exp(1.0j*signedm*phi)
156                dylmp=1.0j*m*y*np.exp(1.0j*signedm*phi)
157              
158                rlist.append(ylm)
159                tlist.append(dylmt)
160                plist.append(dylmp/np.sin(theta))
161              
162                xir, xih = sosh.getradial(l,n)
163                fr=np.interp(r,sosh.rmesh,xir*np.sqrt(sosh.rho))
164                fh=np.interp(r,sosh.rmesh,xih*np.sqrt(sosh.rho))
165              
166                ind = ((l==lmod) & (n==nmod))
167                freq = float(modeparms[ind,2])*freqscale
168              #  rc = float(modeparms[ind,5])/1e8
169 tplarson 1.1 #  hc = float(modeparms[ind,6])/1e8
170                rclist.append(fr)
171                hclist.append(fh)
172                freqlist.append(freq)
173              #  print("Using vertical to horizontal ratio of %f, theoretical is %f." % (rc/hc, float(modeparms[ind,4])))
174                print("Scaled frequency = %f." % (freq))
175              
176              # if using measured mode parameters, this would 
177              # allow you to compute frequency as a function of m
178              # from the fitted a-coefficients
179              #  ai=np.append([0.0],modeparms[ind,12:18]/1000)
180              #  pols=sosh.apols(l,6)
181              #  fx=np.matmul(pols,ai)
182              #  freq=modeparms[ind,2]+fx
183              #  freqlist.append(freq[l+signedm]*freqscale)
184              #  amplist.append(modeparms[ind,3])
185              
186              sumr=0.0
187              sumt=0.0
188              sump=0.0
189              for i in range(nmodes):
190 tplarson 1.1   sumr += rclist[i]*rlist[i]*freqlist[i]
191                sumt += hclist[i]*tlist[i]*freqlist[i]
192                sump += hclist[i]*plist[i]*freqlist[i]
193              
194              vh2=(np.square(sumt.real) + np.square(sump.real))
195              v2=(np.square(sumr.real) + vh2)
196              
197              varlist['Vr']=sumr.real
198              varlist['Vt']=sumt.real
199              varlist['Vp']=sump.real
200              varlist['Vh']=np.sqrt(vh2)
201              varlist['Vmag']=np.sqrt(v2)
202              varlist['Vsq']=v2
203              
204              while True:
205              
206                var=varlist[sosh.plotvar]
207              
208                fig,im = sosh.drawfigure(var, fsize=figsize)
209                if (sosh.isave != 0):
210                  sosh.savefigure(animateflag=animateflag)
211 tplarson 1.1   if (animateflag != 0):
212                  anim = animation.FuncAnimation(fig, sumanimate, frames=nframes, interval=interval)
213                plt.show()
214              
215                print()
216                c=sosh.queryplotparms()
217                if (c == 'q'):
218                  break
219              
220              
221