1 tplarson 1.1 import numpy as np
  2              import os
  3              import drms
  4              from astropy.io import fits
  5              from urllib.request import urlretrieve
  6              import soundfile as sf
  7              
  8              import warnings 
  9              warnings.simplefilter(action='ignore', category=FutureWarning)
 10              
 11              #datadir = '/home/tplarson/solar/data'
 12              datadir = '../data'
 13              
 14              def getmodeparms(instrument=None, day=None, lmin=0, lmax=300):
 15              
 16                if not os.path.exists(datadir):
 17                  os.makedirs(datadir)
 18              
 19                if (instrument == None):
 20                  inst=input("Enter name of instrument: ")
 21                else:
 22 tplarson 1.1     inst=instrument
 23              
 24                if (day == None):
 25                  daystr=input("Enter day number: ")
 26                else:
 27                  daystr=str(day)
 28              
 29                mfile=getmfile(inst.lower(), daystr)
 30                if (mfile == None):
 31                  return None
 32              
 33                msfile=mfile+'.msplit'
 34              
 35                if os.path.exists(msfile):
 36                  print("Msplit file already exists.")
 37                else:
 38                  msfile=mkmsplit(mfile,lmin=lmin,lmax=lmax)
 39                  
 40                return mfile,msfile
 41              
 42              
 43 tplarson 1.1 def getmfile(inst, daystr):
 44              
 45                if (inst == 'hmi'):
 46                  qblank = 'hmi.V_sht_modes[{:d}d][0][300][138240]'
 47                  mblank = 'hmi.%sd.modes'
 48                  avgurl = 'http://sun.stanford.edu/~tplarson/audio/HMI/hmi.average.modes'
 49                elif (inst == 'mdi'):
 50                  qblank = 'mdi.vw_V_sht_modes[{:d}d][0][300][103680]'
 51                  mblank = 'mdi.%sd.modes'
 52                  avgurl = 'http://sun.stanford.edu/~tplarson/audio/MDI/mdi.average.modes'
 53                else:
 54                  print("Invalid instrument, try again.")
 55                  return None
 56              
 57                mblank = '%s.%sd.modes'
 58                mfile=os.path.join(datadir,mblank % (inst, daystr))
 59                if os.path.exists(mfile):
 60                  print("Mode fits file already exists.")
 61                elif (daystr == 'average'):
 62                  try:
 63                    urlretrieve(avgurl, mfile)
 64 tplarson 1.1     except:
 65                    print("Failure of urlretrieve() for %s." % avgurl)
 66                    return None
 67                else:
 68                  try:
 69                    day=int(daystr)
 70                  except:
 71                      print("Invalid number, try again.")
 72                      return None
 73                  qstr=qblank.format(day)
 74                  c = drms.Client()
 75                  m = c.query(qstr, seg='m6')
 76                  if (len(m) == 0):
 77                    print("No data found for that day number, try again.")
 78                    return None
 79                  url = 'http://jsoc.stanford.edu' + m.m6[0]
 80                  try:
 81                    urlretrieve(url, mfile)
 82                  except:
 83                    print("Failure of urlretrieve() for %s." % url)
 84                    return None
 85 tplarson 1.1 
 86                return mfile
 87              
 88              
 89              def mkmsplit(mfile, lmin=0, lmax=300):
 90              
 91                modeparms=np.loadtxt(mfile)
 92                lmod=modeparms[:,0]
 93                nmod=modeparms[:,1]
 94              
 95                outfmt='{:d} {:d} {:d} {:f} {:f}'
 96                outfile = mfile + '.msplit'
 97                file=open(outfile,'w')
 98              
 99                l=lmin
100                while (l <= lmax):
101                  pols=apols(l,6)
102                  ind = (lmod == l)
103                  nlist = nmod[ind]
104                  for n in nlist:
105                    ind2 = (lmod == l) & (nmod == n)
106 tplarson 1.1       ai=np.append([0.0],modeparms[ind2,12:18]/1000)
107                    ei=np.append([0.0],modeparms[ind2,18:24]/1000)
108                    fx=np.matmul(pols,ai)
109                    freq=modeparms[ind2,2]+fx
110                    var=np.matmul((pols*pols),(ei*ei))
111                    sig=modeparms[ind2,7]+np.sqrt(var)
112                    for m in range(-l,l+1):
113                      outstr=outfmt.format(int(l),int(n),int(m),freq[m+l],sig[m+l])
114                      file.write('%s\n' % (outstr))
115                  l+=1
116                file.close()
117                return outfile
118              
119              
120              def getwavfiles(instrument=None, day=None, l=None, m=None, delfits=False):
121              
122                if not os.path.exists(datadir):
123                  os.makedirs(datadir)
124              
125                if (instrument == None):
126                  inst=input("Enter name of instrument: ")
127 tplarson 1.1   else:
128                  inst=instrument
129              
130                if (inst.lower() == 'hmi'):
131                  hind=1
132                elif (inst.lower() == 'mdi'):
133                  hind=0
134                else:
135                  print("Invalid instrument, try again.")
136                  return None
137              
138                if (day == None):
139                  daystr=input("Enter day number: ")
140                  try:
141                    day=int(daystr)
142                  except:
143                    print("Invalid number, try again.")
144                    return None
145                if (l == None):
146                  lstr=input("Enter spherical harmonic degree: ")
147                  try:
148 tplarson 1.1       l=int(lstr)
149                  except:
150                    print("Invalid number, try again.")
151                    return None
152                if (m == None):
153                  mstr=input("Enter azimuthal order: ")
154                  try:
155                    m=int(mstr)
156                    m=np.abs(m)
157                    if (m > l):
158                      print("Abs(m) must be less than l, try again.")
159                      return None
160                  except:
161                    print("Invalid number, try again.")
162                    return None
163              
164                wblank = '%s.%id_l=%i_m=%i_data%s.wav'
165                wfile=os.path.join(datadir,wblank % (inst.lower(),day,l,m,'r'))
166                wfile2=os.path.join(datadir,wblank % (inst.lower(),day,l,m,'i'))
167                if os.path.exists(wfile) and (m==0 or os.path.exists(wfile2)):
168                  if (m == 0):
169 tplarson 1.1       print("Wav file already exists.")
170                    return wfile
171                  else:
172                    print("Wav files already exists.")
173                    return wfile,wfile2
174              
175                ffile=getfitsfile(inst.lower(),day,l)
176                if (ffile == None):
177                  return None
178              
179                wf=convertfitstowav(ffile, hind, m, wfile)
180                if (wf == None):
181                  return None
182              
183                if delfits:
184                  os.remove(ffile)
185                return wf
186              
187              
188              def getfitsfile(inst, day, l):
189              
190 tplarson 1.1   if (inst == 'hmi'):
191                  qblank = 'hmi.V_sht_gf_72d[{:d}d][{:d}][{:d}][138240]'
192                elif (inst == 'mdi'):
193                  qblank = 'mdi.vw_V_sht_gf_72d[{:d}d][{:d}][{:d}][103680]'
194                else:
195                  print("Invalid instrument, try again.")
196                  return None
197              
198                fblank='%s.%id.%i.fits'
199                ffile=os.path.join(datadir,fblank % (inst,day,l))
200                if os.path.exists(ffile):
201                  print("Fits file already exists.")
202                else:
203                  qstr=qblank.format(day,l,l)
204                  c = drms.Client()
205                  f = c.query(qstr,seg='data')
206                  if (len(f) == 0):
207                    print("No data found for that day number and degree, try again.")
208                    return None
209                  url = 'http://jsoc.stanford.edu' + f.data[0]
210                  try:
211 tplarson 1.1       urlretrieve(url, ffile)
212                  except:
213                    print("Failure of urlretrieve() for %s." % url)
214                    return None
215              
216                return ffile
217              
218              
219              def convertfitstowav(ffile, hind, m, wfile):
220              
221                hdul = fits.open(ffile)
222                h=hdul[hind]
223                d=h.data
224                x=d[m,0::2]/np.abs(d[m,:]).max()
225                sf.write(wfile, x, 8000, subtype='FLOAT')
226              
227                if (m > 0):
228                  wfile2=wfile.replace('datar','datai')
229                  x=d[m,1::2]/np.abs(d[m,:]).max()
230                  sf.write(wfile2, x, 8000, subtype='FLOAT')
231              
232 tplarson 1.1   hdul.close()
233                if (m==0):
234                  return wfile
235                else:
236                  return wfile,wfile2
237              
238              
239              def apols(lin, amaxin):
240              
241              # adapted from IDL procedure apols.pro written by Jesper Schou
242              
243                l=np.long(lin)
244                amax=np.minimum(amaxin,2*l)
245              
246                pols=np.zeros((2*l+1,amaxin+1))
247              # It is ESSENTIAL that x is set up such that x(-m)=-x(m) to full machine
248              # accuracy or that the re-orthogonalization is done with respect to all
249              # previous polynomials (second option below).
250              # x=(dindgen(2*l+1)-l)/l
251                x=np.linspace(-1,1,2*l+1)
252              
253 tplarson 1.1   pols[:,0]=1/np.sqrt(2*l+1.0)
254                if (amax >= 1):
255                  pols[:,1]=x/np.sqrt((l+1.0)*(2*l+1.0)/(3.0*l))
256              # for n=2l,amax do begin
257              # Set up polynomials using exact recurrence relation.
258                for n in range(2,amax+1):
259                  a0=2.0*l*(2*n-1.0)/(n*(2*l-n+1))
260                  b0=-(n-1.0)/n*(2*l+n)/(2*l-n+1)
261                  a=a0*np.sqrt((2*n+1.0)/(2*n-1.0)*(2*l-n+1.0)/(2*l+n+1.0))
262                  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))
263                  help=a*x*pols[:,n-1]+b*pols[:,n-2]
264              # Turns out that roundoff kills the algorithm, so we have to re-orthogonalize.
265              # Two choices here. First one is twice as fast and more accurate, but
266              # depends on the rounding being done in the standard IEEE way.
267              #  for j=n-2,0,-2 do begin 
268                  for j in range(n-2,-1,-2):    
269              # This choice is robust to the roundoff, but twice as slow and generally
270              # somewhat less accurate
271              # for j=n-1,0,-1 do begin
272                    help=help-pols[:,j]*np.sum(help*pols[:,j])
273              #  end
274 tplarson 1.1 # Reset norm to 1.
275                  pols[:,n]=help/np.sqrt(np.sum(np.square(help)))
276              
277              # Reset polynomials to have P_l(l)=l by setting overall norm.
278              # Note that this results in more accurate overall normalization, but
279              # that P_l(l) may be very far from l. This is the key to the algorithm.
280                c=l**2*(2*l+1.0)
281              # for n=0l,amax do begin
282                for n in range(0,amax+1):
283                  c=c*(2*l+n+1.0)/(2*l-n+1.0)
284                  pols[:,n]=pols[:,n]*np.sqrt(c/(2*n+1))
285              
286                return pols
287              
288