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 tplarson 1.1 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((pols*pols),(ei*ei))
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 tplarson 1.1 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 tplarson 1.1 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
|
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 tplarson 1.1
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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 # 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)*(2*l+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.0*l*(2*n-1.0)/(n*(2*l-n+1))
269 b0=-(n-1.0)/n*(2*l+n)/(2*l-n+1)
270 a=a0*np.sqrt((2*n+1.0)/(2*n-1.0)*(2*l-n+1.0)/(2*l+n+1.0))
271 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))
272 help=a*x*pols[:,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 tplarson 1.1 # somewhat less accurate
280 # for j=n-1,0,-1 do begin
281 help=help-pols[:,j]*np.sum(help*pols[:,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*(2*l+n+1.0)/(2*l-n+1.0)
293 pols[:,n]=pols[:,n]*np.sqrt(c/(2*n+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 tplarson 1.4 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 tplarson 1.4 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 tplarson 1.4 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 tplarson 1.4 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 tplarson 1.4 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 tplarson 1.4 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 tplarson 1.4 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 tplarson 1.4 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 tplarson 1.4 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 tplarson 1.4
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 tplarson 1.4 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(2*np.pi*f[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 tplarson 1.4 # freqs = fftpack.fftfreq(len(x)) * sr
551
552
|
Return to
soshdata.py
CVS log
Up to [Development] / JSOC / proj / globalhs / sosh