(file) Return to sosh.py CVS log (file) (dir) Up to [Development] / JSOC / proj / globalhs / sosh

Diff for /JSOC/proj/globalhs/sosh/sosh.py between version 1.1 and 1.2

version 1.1, 2018/12/09 13:52:19 version 1.2, 2019/02/22 22:26:26
Line 3  import os
Line 3  import os
 import h5py import h5py
 import matplotlib.pyplot as plt import matplotlib.pyplot as plt
 #from matplotlib import animation #from matplotlib import animation
   from matplotlib import cm
   from matplotlib.colors import ListedColormap
 modeln=np.array([]) modeln=np.array([])
 modell=np.array([]) modell=np.array([])
 modelnu=np.array([]) modelnu=np.array([])
Line 14  R=0.0
Line 17  R=0.0
 def loadmodel(): def loadmodel():
 #  in_dir = '/home/tplarson/solar/mods' #  in_dir = '/home/tplarson/solar/mods'
   in_dir = 'mods'    in_dir = modeldir
   fname = 'mods_p3_eigfcn.h5'   fname = 'mods_p3_eigfcn.h5'
   hf = h5py.File(os.path.join(in_dir, fname))   hf = h5py.File(os.path.join(in_dir, fname))
Line 31  def loadmodel():
Line 34  def loadmodel():
 def getradial(l,n): def getradial(l,n):
 #  in_dir = '/home/tplarson/solar/mods' #  in_dir = '/home/tplarson/solar/mods'
   in_dir = 'mods'    in_dir = modeldir
   fname = 'mods_p3_eigfcn.h5'   fname = 'mods_p3_eigfcn.h5'
   hf = h5py.File(os.path.join(in_dir, fname))   hf = h5py.File(os.path.join(in_dir, fname))
Line 53  def getradial(l,n):
Line 56  def getradial(l,n):
   return (xir, xih)   return (xir, xih)
 def writemodel(lmin=0, lmax=300, rsurf=1.0):  def writesurfacemodel(lmin=0, lmax=300, rsurf=1.0):
   outfile='modes.model.txt'    outfile='model.surface.modes'
   file=open(outfile,'w')   file=open(outfile,'w')
   in_dir = 'mods'    in_dir = modeldir
   fname = 'mods_p3_eigfcn.h5'   fname = 'mods_p3_eigfcn.h5'
   hf = h5py.File(os.path.join(in_dir, fname))   hf = h5py.File(os.path.join(in_dir, fname))
Line 69  def writemodel(lmin=0, lmax=300, rsurf=1
Line 72  def writemodel(lmin=0, lmax=300, rsurf=1
   modeln, modell, modelnu, modelsig2, modelE = [hf['modes'][k].value for k in ['n', 'l', 'nu', 'sigma2', 'E']]   modeln, modell, modelnu, modelsig2, modelE = [hf['modes'][k].value for k in ['n', 'l', 'nu', 'sigma2', 'E']]
   outfmt='{:d} {:d} {:f} {:f} {:f} {:f} {:f} {:f} {:f}'   outfmt='{:d} {:d} {:f} {:f} {:f} {:f} {:f} {:f} {:f}'
   l=0    l=lmin
   while (l <= lmax):   while (l <= lmax):
     ind = (modell == l)     ind = (modell == l)
     nlist = modeln[ind]     nlist = modeln[ind]
Line 197  def image2rtheta(xpixels=1000, ypixels=1
Line 200  def image2rtheta(xpixels=1000, ypixels=1
   xx, yy = np.meshgrid(xx, yy)   xx, yy = np.meshgrid(xx, yy)
   rr = np.sqrt(xx*xx+yy*yy)   rr = np.sqrt(xx*xx+yy*yy)
   index = (rr <= 1.0)   index = (rr <= 1.0)
   lat = np.arctan(yy/np.abs(xx))  
   r = np.ma.array(rr, mask=np.logical_not(index))   r = np.ma.array(rr, mask=np.logical_not(index))
   #  lat = np.arctan(yy/np.abs(xx))
     lat = np.arctan2(yy,np.abs(xx))
   theta = np.ma.array(np.pi/2 - lat, mask=np.logical_not(index))   theta = np.ma.array(np.pi/2 - lat, mask=np.logical_not(index))
   return (r, theta)   return (r, theta)
Line 214  def setplm(l, m, x, plm, dplm):
Line 218  def setplm(l, m, x, plm, dplm):
   eps = 1.0e-12   eps = 1.0e-12
   x1 = np.maximum(np.minimum(x,1-eps),eps-1)   x1 = np.maximum(np.minimum(x,1-eps),eps-1)
   #  x1 = x
   x2 = 1.0/(x1*x1-1.0)   x2 = 1.0/(x1*x1-1.0)
   x3 = x1*x2   x3 = x1*x2
   c = np.sqrt((2*m+1)/2.0)   c = np.sqrt((2*m+1)/2.0)
Line 240  def setplm(l, m, x, plm, dplm):
Line 245  def setplm(l, m, x, plm, dplm):
   return (plm[...,l-m],dplm[...,l-m])   return (plm[...,l-m],dplm[...,l-m])
 def apols(lin, amaxin):  
 # adapted from IDL procedure apols.pro written by Jesper Schou  
 # It is ESSENTIAL that x is set up such that x(-m)=-x(m) to full machine  
 # accuracy or that the re-orthogonalization is done with respect to all  
 # previous polynomials (second option below).  
 # x=(dindgen(2*l+1)-l)/l  
   if (amax >= 1):  
 # for n=2l,amax do begin  
 # Set up polynomials using exact recurrence relation.  
   for n in range(2,amax+1):  
 # Turns out that roundoff kills the algorithm, so we have to re-orthogonalize.  
 # Two choices here. First one is twice as fast and more accurate, but  
 # depends on the rounding being done in the standard IEEE way.  
 #  for j=n-2,0,-2 do begin  
     for j in range(n-2,-1,-2):  
 # This choice is robust to the roundoff, but twice as slow and generally  
 # somewhat less accurate  
 # for j=n-1,0,-1 do begin  
 #  end  
 # Reset norm to 1.  
 # Reset polynomials to have P_l(l)=l by setting overall norm.  
 # Note that this results in more accurate overall normalization, but  
 # that P_l(l) may be very far from l. This is the key to the algorithm.  
 # for n=0l,amax do begin  
   for n in range(0,amax+1):  
   return pols  
 lsave=0 lsave=0
 msave=0 msave=0
 nsave=1 nsave=1
Line 483  def nothing():
Line 438  def nothing():
 nframes=64 nframes=64
 dpi=300 dpi=300
 pngdirfmt = './png/l{:d}_m{:d}_n{:d}_{:s}'  icolshift=0
   pngdirfmt = './png_out/{:s}{:s}.l{:d}m{:d}n{:d}'
 animate=nothing animate=nothing
 def drawfigure(var, fsize=5): def drawfigure(var, fsize=5):
   fig, ax = plt.subplots(num=1, figsize=(fsize, fsize))  #  below was written for single images. same functionality is now
 #  ax.set(xlim=(-lim, lim), ylim=(-lim, lim))  #  achieved by setting vmin and vmax.  probably exactly the same
   #  but new code is clearer.
   #  for animations.
   #  if (icolshift != 0) and plotvar in ['Vr','Vt','Vp']:
   #    mn,mx=var.min(),var.max()
   #    absarr=np.abs([mn,mx])
   #    frac=(mx-mn)/(2*absarr.max())
   #    table=cm.get_cmap(colormap,256/frac)
   #    if (absarr[0]>absarr[1]):
   #      newcolors=table(np.linspace(0,frac,256))
   #    else:
   #      newcolors=table(np.linspace(1-frac,1,256))
   #    newmap=ListedColormap(newcolors, name='soshtmp')
   #    cm.register_cmap(cmap=newmap)
   #    usemap='soshtmp'
   #  else:
   #    usemap=colormap
     fig, ax = plt.subplots(num=1, figsize=(fsize, fsize))
   ax.clear()   ax.clear()
   im=ax.imshow(var,cmap=colormap)  #  im=ax.imshow(var,cmap=usemap)
     im=ax.imshow(var, origin='lower', cmap=colormap)
     if (icolshift == 1) and plotvar in ['Vr','Vt','Vp']:
   #    im=ax.imshow(var, cmap=colormap, vmin=-maxabs, vmax=maxabs)
       im.set_clim(vmin=-maxabs, vmax=maxabs)
       print("Scaling to maxval = %f"%maxabs)
     elif (icolshift == 2):
       for i in range(nframes):
         if (val > maxval):
       if plotvar in ['Vr','Vt','Vp']:
   #      im=ax.imshow(var, cmap=colormap, vmin=-maxval, vmax=maxval)
         im.set_clim(vmin=-maxval, vmax=maxval)
         im.set_clim(vmin=0.0, vmax=maxval)
       print("Scaling to maxval = %f"%maxval)
   #  else:
   #    im=ax.imshow(var, cmap=colormap)
   ax.set_axis_off()   ax.set_axis_off()
   fig.canvas.draw()   fig.canvas.draw()
   return (fig,im)   return (fig,im)
 def savefigure(animateflag=1):  def savefigure(ianimate=1, label=''):
   l=lsave   l=lsave
   m=msave   m=msave
   n=nsave   n=nsave
   if (animateflag == 0):    if (ianimate == 0):
     savestr="l%im%in%i.png" % (l,m,n)      savestr='png_out/'+plotvar+label+'.l%im%in%i.png' % (l,m,n)
     plt.savefig(savestr, dpi=dpi)     plt.savefig(savestr, dpi=dpi)
     print("File saved.")     print("File saved.")
   else:   else:
     pngdir = pngdirfmt.format(l, m, n, plotvar)      pngdir = pngdirfmt.format(plotvar, label, l, m, n)
     if not os.path.exists(pngdir):     if not os.path.exists(pngdir):
       os.makedirs(pngdir)       os.makedirs(pngdir)
     print("Writing files...", end="")     print("Writing files...", end="")
     for i in range(nframes):     for i in range(nframes):
         print(i, end=" ", flush=True)
       animate(i)       animate(i)
       fpath = os.path.join(pngdir, '{:03d}.png'.format(i))       fpath = os.path.join(pngdir, '{:03d}.png'.format(i))
       print(i, end=" ", flush=True)  
       plt.savefig(fpath, dpi=dpi)       plt.savefig(fpath, dpi=dpi)
     print("done.")     print("done.")

Removed from v.1.1  
changed lines
  Added in v.1.2

Karen Tian
Powered by
ViewCVS 0.9.4