JSOC/proj/globalhs/sosh/sosh.py - diff

Return to sosh.py CVS log

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

import h5py

import matplotlib.pyplot as plt

#from matplotlib import animation

from matplotlib import cm

from matplotlib.colors import ListedColormap

modeldir='../mods'

modeln=np.array([])

modell=np.array([])

modelnu=np.array([])

Line 14 R=0.0

Line 17 R=0.0

def loadmodel():

# in_dir = '/home/tplarson/solar/mods'

in_dir = 'mods'

in_dir = modeldir

fname = 'mods_p3_eigfcn.h5'

hf = h5py.File(os.path.join(in_dir, fname))

Line 31 def loadmodel():

Line 34 def loadmodel():

def getradial(l,n):

# in_dir = '/home/tplarson/solar/mods'

in_dir = 'mods'

in_dir = modeldir

fname = 'mods_p3_eigfcn.h5'

hf = h5py.File(os.path.join(in_dir, fname))

Line 53 def getradial(l,n):

Line 56 def getradial(l,n):

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')

in_dir = 'mods'

in_dir = modeldir

fname = 'mods_p3_eigfcn.h5'

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']]

outfmt='{:d} {:d} {:f} {:f} {:f} {:f} {:f} {:f} {:f}'

l=0

l=lmin

while (l <= lmax):

ind = (modell == l)

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)

rr = np.sqrt(xx*xx+yy*yy)

index = (rr <= 1.0)

lat = np.arctan(yy/np.abs(xx))

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))

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

x1 = np.maximum(np.minimum(x,1-eps),eps-1)

# x1 = x

x2 = 1.0/(x1*x1-1.0)

x3 = x1*x2

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])

def apols(lin, amaxin):

# adapted from IDL procedure apols.pro written by Jesper Schou

l=np.long(lin)

amax=np.minimum(amaxin,2*l)

pols=np.zeros((2*l+1,amaxin+1))

# 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

x=np.linspace(-1,1,2*l+1)

pols[:,0]=1/np.sqrt(2*l+1.0)

if (amax >= 1):

pols[:,1]=x/np.sqrt((l+1.0)*(2*l+1.0)/(3.0*l))

# for n=2l,amax do begin

# Set up polynomials using exact recurrence relation.

for n in range(2,amax+1):

a0=2.0*l*(2*n-1.0)/(n*(2*l-n+1))

b0=-(n-1.0)/n*(2*l+n)/(2*l-n+1)

a=a0*np.sqrt((2*n+1.0)/(2*n-1.0)*(2*l-n+1.0)/(2*l+n+1.0))

b=b0*np.sqrt((2*n+1.0)/(2*n-3.0)*(2*l-n+1.0)*(2*l-n+2.0)/(2*l+n+1.0)/(2*l+n))

help=a*x*pols[:,n-1]+b*pols[:,n-2]

# 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

help=help-pols[:,j]*np.sum(help*pols[:,j])

# end

# Reset norm to 1.

pols[:,n]=help/np.sqrt(np.sum(np.square(help)))

# 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.

c=l**2*(2*l+1.0)

# for n=0l,amax do begin

for n in range(0,amax+1):

c=c*(2*l+n+1.0)/(2*l-n+1.0)

pols[:,n]=pols[:,n]*np.sqrt(c/(2*n+1))

return pols

lsave=0

msave=0

nsave=1

Line 483 def nothing():

Line 438 def nothing():

nframes=64

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

calcimage=nothing

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()

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']:

mn,mx=var.min(),var.max()

maxabs=np.abs([mn,mx]).max()

# 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):

maxval=0.0

for i in range(nframes):

d=calcimage(i)

val=np.abs(d).max()

if (val > maxval):

maxval=val

if plotvar in ['Vr','Vt','Vp']:

# im=ax.imshow(var, cmap=colormap, vmin=-maxval, vmax=maxval)

im.set_clim(vmin=-maxval, vmax=maxval)

else:

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()

fig.canvas.draw()

return (fig,im)

def savefigure(animateflag=1):

def savefigure(ianimate=1, label=''):

l=lsave

m=msave

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)

print("File saved.")

else:

pngdir = pngdirfmt.format(l, m, n, plotvar)

pngdir = pngdirfmt.format(plotvar, label, l, m, n)

if not os.path.exists(pngdir):

os.makedirs(pngdir)

print("Writing files...", end="")

for i in range(nframes):

print(i, end=" ", flush=True)

animate(i)

fpath = os.path.join(pngdir, '{:03d}.png'.format(i))

print(i, end=" ", flush=True)

plt.savefig(fpath, dpi=dpi)

print("done.")

Legend:

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

Karen Tian