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