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