globalhs/sosh/soshdata.py

import numpy as np
import os
import re
import drms
from astropy.io import fits
from urllib.request import urlretrieve
import soundfile as sf
from scipy import fftpack

import warnings 
warnings.simplefilter(action='ignore', category=FutureWarning)

#datadir = '/home/tplarson/solar/data'
datadir = '../data'

def getmodeparms(instrument=None, day=None, lmin=0, lmax=300, makemsplit=True):

  if not os.path.exists(datadir):
    os.makedirs(datadir)

  if (instrument == None):
    inst=input("Enter name of instrument: ")
  else:
    inst=instrument

  if (day == None):
    daystr=input("Enter day number: ")
  else:
    daystr=str(day)

  mfile=getmfile(inst.lower(), daystr)
  if (mfile == None):
    return None

  if (makemsplit):
    msfile=mfile+'.msplit'

    if os.path.exists(msfile):
      print("Msplit file already exists.")
    else:
      msfile=mkmsplit(mfile,lmin=lmin,lmax=lmax)

    return mfile,msfile
  else:
    return mfile


def getmfile(inst, daystr):

  if (inst == 'hmi'):
    qblank = 'hmi.V_sht_modes[{:d}d][0][300][138240]'
    avgurl = 'http://sun.stanford.edu/~tplarson/audio/HMI/hmi.average.modes'
  elif (inst == 'mdi'):
    qblank = 'mdi.vw_V_sht_modes[{:d}d][0][300][103680]'
    avgurl = 'http://sun.stanford.edu/~tplarson/audio/MDI/mdi.average.modes'
  else:
    print("Invalid instrument, try again.")
    return None

  if (daystr == 'average'):
    mblank = '%s.average.modes'
    mfile=os.path.join(datadir,mblank % inst)
    if os.path.exists(mfile):
      print("Mode fits file already exists.")
    else:
      try:
        urlretrieve(avgurl, mfile)
      except:
        print("Failure of urlretrieve() for %s." % avgurl)
        return None
  else:
    mblank = '%s.%sd.modes'
    mfile=os.path.join(datadir,mblank % (inst, daystr))
    if os.path.exists(mfile):
      print("Mode fits file already exists.")
    else:
      try:
        day=int(daystr)
      except:
        print("Invalid number, try again.")
        return None
      qstr=qblank.format(day)
      c = drms.Client()
      m = c.query(qstr, seg='m6')
      if (len(m) == 0):
        print("No data found for that day number, try again.")
        return None
      url = 'http://jsoc.stanford.edu' + m.m6[0]
      try:
        urlretrieve(url, mfile)
      except:
        print("Failure of urlretrieve() for %s." % url)
        return None

  return mfile


def mkmsplit(mfile, lmin=0, lmax=300):

  modeparms=np.loadtxt(mfile)
  lmod=modeparms[:,0]
  nmod=modeparms[:,1]

  outfmt='{:d} {:d} {:d} {:f} {:f}'
  outfile = mfile + '.msplit'
  file=open(outfile,'w')

  l=lmin
  while (l <= lmax):
    pols=apols(l,6)
    ind = (lmod == l)
    nlist = nmod[ind]
    for n in nlist:
      ind2 = (lmod == l) & (nmod == n)
      ai=np.append([0.0],modeparms[ind2,12:18]/1000)
      ei=np.append([0.0],modeparms[ind2,18:24]/1000)
      fx=np.matmul(pols,ai)
      freq=modeparms[ind2,2]+fx
      var=np.matmul((pols*pols),(ei*ei))
      sig=modeparms[ind2,7]+np.sqrt(var)
      for m in range(-l,l+1):
        outstr=outfmt.format(int(l),int(n),int(m),freq[m+l],sig[m+l])
        file.write('%s\n' % (outstr))
    l+=1
  file.close()
  return outfile


def getwavfiles(instrument=None, day=None, l=None, m=None, delfits=False):

  if not os.path.exists(datadir):
    os.makedirs(datadir)

  if (instrument == None):
    inst=input("Enter name of instrument: ")
  else:
    inst=instrument

  if (inst.lower() == 'hmi'):
    hind=1
  elif (inst.lower() == 'mdi'):
    hind=0
  else:
    print("Invalid instrument, try again.")
    return None

  if (day == None):
    daystr=input("Enter day number: ")
    try:
      day=int(daystr)
    except:
      print("Invalid number, try again.")
      return None
  if (l == None):
    lstr=input("Enter spherical harmonic degree: ")
    try:
      l=int(lstr)
    except:
      print("Invalid number, try again.")
      return None
  if (m == None):
    mstr=input("Enter azimuthal order: ")
    try:
      m=int(mstr)
      m=np.abs(m)
      if (m > l):
        print("Abs(m) must be less than l, try again.")
        return None
    except:
      print("Invalid number, try again.")
      return None

  wblank = '%s.%id_72d.l=%i_m=%i_data%s.wav'
  wfile=os.path.join(datadir,wblank % (inst.lower(),day,l,m,'r'))
  wfile2=os.path.join(datadir,wblank % (inst.lower(),day,l,m,'i'))
  if os.path.exists(wfile) and (m==0 or os.path.exists(wfile2)):
    if (m == 0):
      print("Wav file already exists.")
      return wfile
    else:
      print("Wav files already exists.")
      return wfile,wfile2

  ffile=getfitsfile(inst.lower(),day,l)
  if (ffile == None):
    return None

  wf=convertfitstowav(ffile, hind, m, wfile)
  if (wf == None):
    return None

  if delfits:
    os.remove(ffile)
  return wf


def getfitsfile(inst, day, l):

  if (inst == 'hmi'):
    qblank = 'hmi.V_sht_gf_72d[{:d}d][{:d}][{:d}][138240]'
  elif (inst == 'mdi'):
    qblank = 'mdi.vw_V_sht_gf_72d[{:d}d][{:d}][{:d}][103680]'
  else:
    print("Invalid instrument, try again.")
    return None

  fblank='%s.%id.%i.fits'
  ffile=os.path.join(datadir,fblank % (inst,day,l))
  if os.path.exists(ffile):
    print("Fits file already exists.")
  else:
    qstr=qblank.format(day,l,l)
    c = drms.Client()
    f = c.query(qstr,seg='data')
    if (len(f) == 0):
      print("No data found for that day number and degree, try again.")
      return None
    url = 'http://jsoc.stanford.edu' + f.data[0]
    try:
      urlretrieve(url, ffile)
    except:
      print("Failure of urlretrieve() for %s." % url)
      return None

  return ffile


def convertfitstowav(ffile, hind, m, wfile):

  hdul = fits.open(ffile)
  h=hdul[hind]
  d=h.data
  x=d[m,0::2]/np.abs(d[m,:]).max()
  sf.write(wfile, x, 8000, subtype='FLOAT')

  if (m > 0):
    wfile2=wfile.replace('datar','datai')
    x=d[m,1::2]/np.abs(d[m,:]).max()
    sf.write(wfile2, x, 8000, subtype='FLOAT')

  hdul.close()
  if (m==0):
    return wfile
  else:
    return wfile,wfile2


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


def splitwavfile(wavfile=None, splitfactor=2, splitboth=True):

  wblank = '%s.%id_%id.l=%i_m=%i_data%s.wav'
  outlist=[]

  if (wavfile == None):
    inst=input("Enter name of instrument: ")
    inst=inst.lower()
    if (inst != 'mdi' and inst != 'hmi'):
      print("Invalid instrument, try again.")
      return None
    daystr=input("Enter day number: ")
    try:
      firstday=int(daystr)
    except:
      print("Invalid number, try again.")
      return None
    ndstr=input("Enter number of days: ")
    try:
      ndaysin=int(ndstr)
    except:
      print("Invalid number, try again.")
      return None
    lstr=input("Enter spherical harmonic degree: ")
    try:
      l=int(lstr)
      if (l < 0):
        print("Degree l must be nonnegative, try again.")
        return None
    except:
      print("Invalid number, try again.")
      return None
    mstr=input("Enter azimuthal order: ")
    try:
      m=int(mstr)
      m=np.abs(m)
      if (m > l):
        print("Abs(m) must be less than l, try again.")
        return None
    except:
      print("Invalid number, try again.")
      return None
    ri='r'
    wavfile=os.path.join(datadir,wblank % (inst,firstday,ndaysin,l,m,ri))
  else:
    regmatch = re.search(r"(\S+)\.(\d+)d_(\d+)d\.l=(\d+)_m=(\d+)_data(\S+).wav", wavfile)
    if (regmatch == None):
      print("Invalid file name, try again.")
      return None
    else:
      inst=regmatch.group(1)
      firstday=int(regmatch.group(2))
      ndaysin=int(regmatch.group(3))
      l=int(regmatch.group(4))
      m=int(regmatch.group(5))
      ri=regmatch.group(6)

  test= ndaysin % splitfactor
  if (test != 0):
    print("Timeseries not divisible by that factor.")
    return None

#  wfile = os.path.join(datadir,wavfile)
  try:
    x, sr = sf.read(wavfile)
  except:
    print("Invalid file, try again.")
    return None

  if (inst == 'mdi'):
    nperday = 1440
  elif (inst == 'hmi'):
    nperday = 1920

  ndaysout = ndaysin / splitfactor
  ndt=int(ndaysout*nperday)
  day = firstday
  index=0
  while (day < firstday+ndaysin):
    xout=x[index:index+ndt]
    wfileout = os.path.join(datadir,wblank % (inst,day,ndaysout,l,m,ri))
    sf.write(wfileout, xout, sr, subtype='FLOAT')
    day += ndaysout
    index += ndt
    outlist.append(wfileout)

  if (splitboth and m > 0):
    if (ri == 'i'):
      ir = 'r'
    elif (ri == 'r'):
      ir = 'i'
    wfile=wblank % (inst,firstday,ndaysin,l,m,ir)
    wavfile = os.path.join(datadir,wfile)
    try:
      x, sr = sf.read(wavfile)
    except:
      print("Invalid file, try again.")
      return None
    day=firstday
    index=0
    while (day < firstday+ndaysin):
      xout=x[index:index+ndt]
      wfileout = os.path.join(datadir,wblank % (inst,day,ndaysout,l,m,ir))
      sf.write(wfileout, xout, sr, subtype='FLOAT')
      day += ndaysout
      index += ndt
      outlist.append(wfileout)

  return outlist


def scanspectrum(wavfile=None,downshift=4,firstbin=None,lastbin=None,nbins=200,ndt=10000,binstep=None,ramp=50):

  if (binstep == None):
    binstep = int(nbins/2)

  wblank = '%s.%id_%id.l=%i_m=%i_data%s.wav'
  if (wavfile == None):
    inst=input("Enter name of instrument: ")
    inst=inst.lower()
    if (inst != 'mdi' and inst != 'hmi'):
      print("Invalid instrument, try again.")
      return None
    daystr=input("Enter day number: ")
    try:
      firstday=int(daystr)
    except:
      print("Invalid number, try again.")
      return None
    ndstr=input("Enter number of days: ")
    try:
      ndaysin=int(ndstr)
    except:
      print("Invalid number, try again.")
      return None
    lstr=input("Enter spherical harmonic degree: ")
    try:
      l=int(lstr)
      if (l < 0):
        print("Degree l must be nonnegative, try again.")
        return None
    except:
      print("Invalid number, try again.")
      return None
    mstr=input("Enter azimuthal order: ")
    try:
      m=int(mstr)
      m=np.abs(m)
      if (m > l):
        print("Abs(m) must be less than l, try again.")
        return None
    except:
      print("Invalid number, try again.")
      return None
    ri='r'
    wavfile=os.path.join(datadir,wblank % (inst,firstday,ndaysin,l,m,ri))
    print("Wav file is %s" % wavfile)
  else:
    regmatch = re.search(r"(\S+)\.(\d+)d_(\d+)d\.l=(\d+)_m=(\d+)_data(\S+).wav", wavfile)
    if (regmatch == None):
      print("Invalid file name, try again.")
      return None
    else:
      inst=regmatch.group(1)
      firstday=int(regmatch.group(2))
      ndaysin=int(regmatch.group(3))
      l=int(regmatch.group(4))
      m=int(regmatch.group(5))
      ri=regmatch.group(6)

  try:
    x, sr = sf.read(wavfile)
  except:
    print("Invalid file, try again.")
    return None

  z=np.zeros((len(x)),dtype=np.complex_)
  if (ri == 'r'):
    z.real=x
  else:
    z.imag=x

  if (m > 0):
    if (ri == 'i'):
      ir = 'r'
    elif (ri == 'r'):
      ir = 'i'
    wfile=wblank % (inst,firstday,ndaysin,l,m,ir)
    wavfile = os.path.join(datadir,wfile)
    try:
      x, sr = sf.read(wavfile)
    except:
      print("Invalid file, try again.")
      return None
  else:
    x=0.0*x

  if (ri == 'r'):
    z.imag=x
  else:
    z.real=x

  window=np.ones(ndt)
  if (ramp != 0):
    rampsamples=int(sr*ramp/1000)
    if (rampsamples > ndt/2):
      print("Ramp longer than ndt/2, try again.")
      return None
    window[0:rampsamples]=np.linspace(0.0,1.0,rampsamples,endpoint=False)
    window[ndt-1:ndt-rampsamples-1:-1]=np.linspace(0.0,1.0,rampsamples,endpoint=False)

  nx=len(x)
  nx2=int(nx/2)
  if (firstbin == None):
    firstbin = int(nx2/6)
  if (lastbin == None):
    lastbin = int(0.9*nx2)
  ftz = fftpack.fft(z)
  freq = fftpack.fftfreq(nx) #*sr
  if (m >= 0):
    mag = np.abs(ftz[0:nx2])
    phase = np.arctan2(ftz[0:nx2].imag,ftz[0:nx2].real)
    f = freq[0:nx2]/downshift
  else:
    mag = np.abs(ftz[nx:nx2:-1])
    phase = np.arctan2(-1*ftz[nx:nx2:-1].imag,ftz[nx:nx2:-1].real)
    f = -1*freq[nx:nx2:-1]/downshift

  nstep=int((lastbin-firstbin)/binstep)
  tlen=nstep*ndt
  t=np.linspace(0,ndt-1,ndt)
  xout=np.zeros((tlen))
  bindex=firstbin
  tindex=0
  t0in=0
  while (bindex + nbins <= lastbin):
    for i in range(nbins):
      xout[tindex:tindex+ndt]+=mag[bindex+i]*np.sin(2*np.pi*f[bindex+i]*(t+tindex)+phase[bindex+i])
    xout[tindex:tindex+ndt]*=window
    bindex+=binstep
    tindex+=ndt

  mx=np.abs(xout).max()
  xout /= mx
  sf.write('output.wav', xout, sr, subtype='FLOAT')

  return xout


#def plotspectrum(wavfile):
#  x, sr = sf.read(wavfile)
#  ftx = fftpack.fft(x)
#  freqs = fftpack.fftfreq(len(x)) * sr


Revision:	1.4
Committed:	Mon Apr 8 16:52:19 2019 UTC (4 years, 5 months ago) by tplarson
Content type:	text/x-python
Branch:	MAIN
CVS Tags:	Ver_9-5, Ver_9-4, Ver_LATEST, Ver_9-41, HEAD
Changes since 1.3:	+258 -1 lines
Log Message:	new versions
#	User	Rev	Content
1	tplarson	1.1	import numpy as np
2			import os
3	tplarson	1.4	import re
4	tplarson	1.1	import drms
5			from astropy.io import fits
6			from urllib.request import urlretrieve
7			import soundfile as sf
8	tplarson	1.4	from scipy import fftpack
9	tplarson	1.1
10			import warnings
11			warnings.simplefilter(action='ignore', category=FutureWarning)
12
13			#datadir = '/home/tplarson/solar/data'
14			datadir = '../data'
15
16	tplarson	1.2	def getmodeparms(instrument=None, day=None, lmin=0, lmax=300, makemsplit=True):
17	tplarson	1.1
18			if not os.path.exists(datadir):
19			os.makedirs(datadir)
20
21			if (instrument == None):
22			inst=input("Enter name of instrument: ")
23			else:
24			inst=instrument
25
26			if (day == None):
27			daystr=input("Enter day number: ")
28			else:
29			daystr=str(day)
30
31			mfile=getmfile(inst.lower(), daystr)
32			if (mfile == None):
33			return None
34
35	tplarson	1.2	if (makemsplit):
36			msfile=mfile+'.msplit'
37	tplarson	1.1
38	tplarson	1.2	if os.path.exists(msfile):
39			print("Msplit file already exists.")
40			else:
41			msfile=mkmsplit(mfile,lmin=lmin,lmax=lmax)
42
43			return mfile,msfile
44	tplarson	1.1	else:
45	tplarson	1.2	return mfile
46	tplarson	1.1
47
48			def getmfile(inst, daystr):
49
50			if (inst == 'hmi'):
51			qblank = 'hmi.V_sht_modes[{:d}d][0][300][138240]'
52			avgurl = 'http://sun.stanford.edu/~tplarson/audio/HMI/hmi.average.modes'
53			elif (inst == 'mdi'):
54			qblank = 'mdi.vw_V_sht_modes[{:d}d][0][300][103680]'
55			avgurl = 'http://sun.stanford.edu/~tplarson/audio/MDI/mdi.average.modes'
56			else:
57			print("Invalid instrument, try again.")
58			return None
59
60	tplarson	1.3	if (daystr == 'average'):
61			mblank = '%s.average.modes'
62			mfile=os.path.join(datadir,mblank % inst)
63			if os.path.exists(mfile):
64			print("Mode fits file already exists.")
65			else:
66			try:
67			urlretrieve(avgurl, mfile)
68			except:
69			print("Failure of urlretrieve() for %s." % avgurl)
70			return None
71	tplarson	1.1	else:
72	tplarson	1.3	mblank = '%s.%sd.modes'
73			mfile=os.path.join(datadir,mblank % (inst, daystr))
74			if os.path.exists(mfile):
75			print("Mode fits file already exists.")
76			else:
77			try:
78			day=int(daystr)
79			except:
80	tplarson	1.1	print("Invalid number, try again.")
81			return None
82	tplarson	1.3	qstr=qblank.format(day)
83			c = drms.Client()
84			m = c.query(qstr, seg='m6')
85			if (len(m) == 0):
86			print("No data found for that day number, try again.")
87			return None
88			url = 'http://jsoc.stanford.edu' + m.m6[0]
89			try:
90			urlretrieve(url, mfile)
91			except:
92			print("Failure of urlretrieve() for %s." % url)
93			return None
94	tplarson	1.1
95			return mfile
96
97
98			def mkmsplit(mfile, lmin=0, lmax=300):
99
100			modeparms=np.loadtxt(mfile)
101			lmod=modeparms[:,0]
102			nmod=modeparms[:,1]
103
104			outfmt='{:d} {:d} {:d} {:f} {:f}'
105			outfile = mfile + '.msplit'
106			file=open(outfile,'w')
107
108			l=lmin
109			while (l <= lmax):
110			pols=apols(l,6)
111			ind = (lmod == l)
112			nlist = nmod[ind]
113			for n in nlist:
114			ind2 = (lmod == l) & (nmod == n)
115			ai=np.append([0.0],modeparms[ind2,12:18]/1000)
116			ei=np.append([0.0],modeparms[ind2,18:24]/1000)
117			fx=np.matmul(pols,ai)
118			freq=modeparms[ind2,2]+fx
119			var=np.matmul((polspols),(eiei))
120			sig=modeparms[ind2,7]+np.sqrt(var)
121			for m in range(-l,l+1):
122			outstr=outfmt.format(int(l),int(n),int(m),freq[m+l],sig[m+l])
123			file.write('%s\n' % (outstr))
124			l+=1
125			file.close()
126			return outfile
127
128
129			def getwavfiles(instrument=None, day=None, l=None, m=None, delfits=False):
130
131			if not os.path.exists(datadir):
132			os.makedirs(datadir)
133
134			if (instrument == None):
135			inst=input("Enter name of instrument: ")
136			else:
137			inst=instrument
138
139			if (inst.lower() == 'hmi'):
140			hind=1
141			elif (inst.lower() == 'mdi'):
142			hind=0
143			else:
144			print("Invalid instrument, try again.")
145			return None
146
147			if (day == None):
148			daystr=input("Enter day number: ")
149			try:
150			day=int(daystr)
151			except:
152			print("Invalid number, try again.")
153			return None
154			if (l == None):
155			lstr=input("Enter spherical harmonic degree: ")
156			try:
157			l=int(lstr)
158			except:
159			print("Invalid number, try again.")
160			return None
161			if (m == None):
162			mstr=input("Enter azimuthal order: ")
163			try:
164			m=int(mstr)
165			m=np.abs(m)
166			if (m > l):
167			print("Abs(m) must be less than l, try again.")
168			return None
169			except:
170			print("Invalid number, try again.")
171			return None
172
173	tplarson	1.4	wblank = '%s.%id_72d.l=%i_m=%i_data%s.wav'
174	tplarson	1.1	wfile=os.path.join(datadir,wblank % (inst.lower(),day,l,m,'r'))
175			wfile2=os.path.join(datadir,wblank % (inst.lower(),day,l,m,'i'))
176			if os.path.exists(wfile) and (m==0 or os.path.exists(wfile2)):
177			if (m == 0):
178			print("Wav file already exists.")
179			return wfile
180			else:
181			print("Wav files already exists.")
182			return wfile,wfile2
183
184			ffile=getfitsfile(inst.lower(),day,l)
185			if (ffile == None):
186			return None
187
188			wf=convertfitstowav(ffile, hind, m, wfile)
189			if (wf == None):
190			return None
191
192			if delfits:
193			os.remove(ffile)
194			return wf
195
196
197			def getfitsfile(inst, day, l):
198
199			if (inst == 'hmi'):
200			qblank = 'hmi.V_sht_gf_72d[{:d}d][{:d}][{:d}][138240]'
201			elif (inst == 'mdi'):
202			qblank = 'mdi.vw_V_sht_gf_72d[{:d}d][{:d}][{:d}][103680]'
203			else:
204			print("Invalid instrument, try again.")
205			return None
206
207			fblank='%s.%id.%i.fits'
208			ffile=os.path.join(datadir,fblank % (inst,day,l))
209			if os.path.exists(ffile):
210			print("Fits file already exists.")
211			else:
212			qstr=qblank.format(day,l,l)
213			c = drms.Client()
214			f = c.query(qstr,seg='data')
215			if (len(f) == 0):
216			print("No data found for that day number and degree, try again.")
217			return None
218			url = 'http://jsoc.stanford.edu' + f.data[0]
219			try:
220			urlretrieve(url, ffile)
221			except:
222			print("Failure of urlretrieve() for %s." % url)
223			return None
224
225			return ffile
226
227
228			def convertfitstowav(ffile, hind, m, wfile):
229
230			hdul = fits.open(ffile)
231			h=hdul[hind]
232			d=h.data
233			x=d[m,0::2]/np.abs(d[m,:]).max()
234			sf.write(wfile, x, 8000, subtype='FLOAT')
235
236			if (m > 0):
237			wfile2=wfile.replace('datar','datai')
238			x=d[m,1::2]/np.abs(d[m,:]).max()
239			sf.write(wfile2, x, 8000, subtype='FLOAT')
240
241			hdul.close()
242			if (m==0):
243			return wfile
244			else:
245			return wfile,wfile2
246
247
248			def apols(lin, amaxin):
249
250			# adapted from IDL procedure apols.pro written by Jesper Schou
251
252			l=np.long(lin)
253			amax=np.minimum(amaxin,2*l)
254
255			pols=np.zeros((2*l+1,amaxin+1))
256			# It is ESSENTIAL that x is set up such that x(-m)=-x(m) to full machine
257			# accuracy or that the re-orthogonalization is done with respect to all
258			# previous polynomials (second option below).
259			# x=(dindgen(2*l+1)-l)/l
260			x=np.linspace(-1,1,2*l+1)
261
262			pols[:,0]=1/np.sqrt(2*l+1.0)
263			if (amax >= 1):
264			pols[:,1]=x/np.sqrt((l+1.0)(2l+1.0)/(3.0*l))
265			# for n=2l,amax do begin
266			# Set up polynomials using exact recurrence relation.
267			for n in range(2,amax+1):
268			a0=2.0l(2n-1.0)/(n(2*l-n+1))
269			b0=-(n-1.0)/n(2l+n)/(2*l-n+1)
270			a=a0np.sqrt((2n+1.0)/(2n-1.0)(2l-n+1.0)/(2l+n+1.0))
271			b=b0np.sqrt((2n+1.0)/(2n-3.0)(2l-n+1.0)(2l-n+2.0)/(2l+n+1.0)/(2*l+n))
272			help=axpols[:,n-1]+b*pols[:,n-2]
273			# Turns out that roundoff kills the algorithm, so we have to re-orthogonalize.
274			# Two choices here. First one is twice as fast and more accurate, but
275			# depends on the rounding being done in the standard IEEE way.
276			# for j=n-2,0,-2 do begin
277			for j in range(n-2,-1,-2):
278			# This choice is robust to the roundoff, but twice as slow and generally
279			# somewhat less accurate
280			# for j=n-1,0,-1 do begin
281			help=help-pols[:,j]np.sum(helppols[:,j])
282			# end
283			# Reset norm to 1.
284			pols[:,n]=help/np.sqrt(np.sum(np.square(help)))
285
286			# Reset polynomials to have P_l(l)=l by setting overall norm.
287			# Note that this results in more accurate overall normalization, but
288			# that P_l(l) may be very far from l. This is the key to the algorithm.
289			c=l*2(2*l+1.0)
290			# for n=0l,amax do begin
291			for n in range(0,amax+1):
292			c=c(2l+n+1.0)/(2*l-n+1.0)
293			pols[:,n]=pols[:,n]np.sqrt(c/(2n+1))
294
295			return pols
296
297
298	tplarson	1.4	def splitwavfile(wavfile=None, splitfactor=2, splitboth=True):
299
300			wblank = '%s.%id_%id.l=%i_m=%i_data%s.wav'
301			outlist=[]
302
303			if (wavfile == None):
304			inst=input("Enter name of instrument: ")
305			inst=inst.lower()
306			if (inst != 'mdi' and inst != 'hmi'):
307			print("Invalid instrument, try again.")
308			return None
309			daystr=input("Enter day number: ")
310			try:
311			firstday=int(daystr)
312			except:
313			print("Invalid number, try again.")
314			return None
315			ndstr=input("Enter number of days: ")
316			try:
317			ndaysin=int(ndstr)
318			except:
319			print("Invalid number, try again.")
320			return None
321			lstr=input("Enter spherical harmonic degree: ")
322			try:
323			l=int(lstr)
324			if (l < 0):
325			print("Degree l must be nonnegative, try again.")
326			return None
327			except:
328			print("Invalid number, try again.")
329			return None
330			mstr=input("Enter azimuthal order: ")
331			try:
332			m=int(mstr)
333			m=np.abs(m)
334			if (m > l):
335			print("Abs(m) must be less than l, try again.")
336			return None
337			except:
338			print("Invalid number, try again.")
339			return None
340			ri='r'
341			wavfile=os.path.join(datadir,wblank % (inst,firstday,ndaysin,l,m,ri))
342			else:
343			regmatch = re.search(r"(\S+)\.(\d+)d_(\d+)d\.l=(\d+)_m=(\d+)_data(\S+).wav", wavfile)
344			if (regmatch == None):
345			print("Invalid file name, try again.")
346			return None
347			else:
348			inst=regmatch.group(1)
349			firstday=int(regmatch.group(2))
350			ndaysin=int(regmatch.group(3))
351			l=int(regmatch.group(4))
352			m=int(regmatch.group(5))
353			ri=regmatch.group(6)
354
355			test= ndaysin % splitfactor
356			if (test != 0):
357			print("Timeseries not divisible by that factor.")
358			return None
359
360			# wfile = os.path.join(datadir,wavfile)
361			try:
362			x, sr = sf.read(wavfile)
363			except:
364			print("Invalid file, try again.")
365			return None
366
367			if (inst == 'mdi'):
368			nperday = 1440
369			elif (inst == 'hmi'):
370			nperday = 1920
371
372			ndaysout = ndaysin / splitfactor
373			ndt=int(ndaysout*nperday)
374			day = firstday
375			index=0
376			while (day < firstday+ndaysin):
377			xout=x[index:index+ndt]
378			wfileout = os.path.join(datadir,wblank % (inst,day,ndaysout,l,m,ri))
379			sf.write(wfileout, xout, sr, subtype='FLOAT')
380			day += ndaysout
381			index += ndt
382			outlist.append(wfileout)
383
384			if (splitboth and m > 0):
385			if (ri == 'i'):
386			ir = 'r'
387			elif (ri == 'r'):
388			ir = 'i'
389			wfile=wblank % (inst,firstday,ndaysin,l,m,ir)
390			wavfile = os.path.join(datadir,wfile)
391			try:
392			x, sr = sf.read(wavfile)
393			except:
394			print("Invalid file, try again.")
395			return None
396			day=firstday
397			index=0
398			while (day < firstday+ndaysin):
399			xout=x[index:index+ndt]
400			wfileout = os.path.join(datadir,wblank % (inst,day,ndaysout,l,m,ir))
401			sf.write(wfileout, xout, sr, subtype='FLOAT')
402			day += ndaysout
403			index += ndt
404			outlist.append(wfileout)
405
406			return outlist
407
408
409			def scanspectrum(wavfile=None,downshift=4,firstbin=None,lastbin=None,nbins=200,ndt=10000,binstep=None,ramp=50):
410
411			if (binstep == None):
412			binstep = int(nbins/2)
413
414			wblank = '%s.%id_%id.l=%i_m=%i_data%s.wav'
415			if (wavfile == None):
416			inst=input("Enter name of instrument: ")
417			inst=inst.lower()
418			if (inst != 'mdi' and inst != 'hmi'):
419			print("Invalid instrument, try again.")
420			return None
421			daystr=input("Enter day number: ")
422			try:
423			firstday=int(daystr)
424			except:
425			print("Invalid number, try again.")
426			return None
427			ndstr=input("Enter number of days: ")
428			try:
429			ndaysin=int(ndstr)
430			except:
431			print("Invalid number, try again.")
432			return None
433			lstr=input("Enter spherical harmonic degree: ")
434			try:
435			l=int(lstr)
436			if (l < 0):
437			print("Degree l must be nonnegative, try again.")
438			return None
439			except:
440			print("Invalid number, try again.")
441			return None
442			mstr=input("Enter azimuthal order: ")
443			try:
444			m=int(mstr)
445			m=np.abs(m)
446			if (m > l):
447			print("Abs(m) must be less than l, try again.")
448			return None
449			except:
450			print("Invalid number, try again.")
451			return None
452			ri='r'
453			wavfile=os.path.join(datadir,wblank % (inst,firstday,ndaysin,l,m,ri))
454			print("Wav file is %s" % wavfile)
455			else:
456			regmatch = re.search(r"(\S+)\.(\d+)d_(\d+)d\.l=(\d+)_m=(\d+)_data(\S+).wav", wavfile)
457			if (regmatch == None):
458			print("Invalid file name, try again.")
459			return None
460			else:
461			inst=regmatch.group(1)
462			firstday=int(regmatch.group(2))
463			ndaysin=int(regmatch.group(3))
464			l=int(regmatch.group(4))
465			m=int(regmatch.group(5))
466			ri=regmatch.group(6)
467
468			try:
469			x, sr = sf.read(wavfile)
470			except:
471			print("Invalid file, try again.")
472			return None
473
474			z=np.zeros((len(x)),dtype=np.complex_)
475			if (ri == 'r'):
476			z.real=x
477			else:
478			z.imag=x
479
480			if (m > 0):
481			if (ri == 'i'):
482			ir = 'r'
483			elif (ri == 'r'):
484			ir = 'i'
485			wfile=wblank % (inst,firstday,ndaysin,l,m,ir)
486			wavfile = os.path.join(datadir,wfile)
487			try:
488			x, sr = sf.read(wavfile)
489			except:
490			print("Invalid file, try again.")
491			return None
492			else:
493			x=0.0*x
494
495			if (ri == 'r'):
496			z.imag=x
497			else:
498			z.real=x
499
500			window=np.ones(ndt)
501			if (ramp != 0):
502			rampsamples=int(sr*ramp/1000)
503			if (rampsamples > ndt/2):
504			print("Ramp longer than ndt/2, try again.")
505			return None
506			window[0:rampsamples]=np.linspace(0.0,1.0,rampsamples,endpoint=False)
507			window[ndt-1:ndt-rampsamples-1:-1]=np.linspace(0.0,1.0,rampsamples,endpoint=False)
508
509			nx=len(x)
510			nx2=int(nx/2)
511			if (firstbin == None):
512			firstbin = int(nx2/6)
513			if (lastbin == None):
514			lastbin = int(0.9*nx2)
515			ftz = fftpack.fft(z)
516			freq = fftpack.fftfreq(nx) #*sr
517			if (m >= 0):
518			mag = np.abs(ftz[0:nx2])
519			phase = np.arctan2(ftz[0:nx2].imag,ftz[0:nx2].real)
520			f = freq[0:nx2]/downshift
521			else:
522			mag = np.abs(ftz[nx:nx2:-1])
523			phase = np.arctan2(-1*ftz[nx:nx2:-1].imag,ftz[nx:nx2:-1].real)
524			f = -1*freq[nx:nx2:-1]/downshift
525
526			nstep=int((lastbin-firstbin)/binstep)
527			tlen=nstep*ndt
528			t=np.linspace(0,ndt-1,ndt)
529			xout=np.zeros((tlen))
530			bindex=firstbin
531			tindex=0
532			t0in=0
533			while (bindex + nbins <= lastbin):
534			for i in range(nbins):
535			xout[tindex:tindex+ndt]+=mag[bindex+i]np.sin(2np.pif[bindex+i](t+tindex)+phase[bindex+i])
536			xout[tindex:tindex+ndt]*=window
537			bindex+=binstep
538			tindex+=ndt
539
540			mx=np.abs(xout).max()
541			xout /= mx
542			sf.write('output.wav', xout, sr, subtype='FLOAT')
543
544			return xout
545
546
547			#def plotspectrum(wavfile):
548			# x, sr = sf.read(wavfile)
549			# ftx = fftpack.fft(x)
550			# freqs = fftpack.fftfreq(len(x)) * sr
551
552