Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/solar_system_setjy.py: 70%
367 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
1#!/usr/bin/python -u
2#
3# bryan butler
4# nrao
5# spring 2012, summer 2016
6#
7# python functions to return expected flux density from solar system
8# bodies. the flux density depends on the geometry (distance, size of
9# body, subearth latitude), and on the model brightness temperature.
10# uncertainties on the flux density can also be returned, but are all
11# set to 0.0 for now, because i don't have uncertainties on the model
12# brightness temperatures.
13#
14# the model brightness temperatures for the various bodies are taken
15# from a combination of modern models and historical observations. see
16# the written description for a more complete description of the model
17# for each body.
18#
19# for all of the bodies but Mars, the model is contained in a text file
20# that has tabulated brightness temperature as a function of frequency.
21# for Mars it is also a function of time. eventually, the full-up
22# model calculations should be in the code (for those bodies that have
23# proper models) but for now, just live with the tabulated versions.
24#
25# version 2.0
26# last edited: 2016Aug15
27# Modified by TT to avoid uncessary file open: 2012Dec13
30from numpy import searchsorted
31from numpy import array
32from scipy.interpolate import interp1d
33from math import exp, pi, cos, sin, isnan, sqrt
34import os
36from casatools import table, measures, quanta, ctsys
37qa = quanta( )
38_tb = table( )
39_me = measures( )
41ctsys_resolve = ctsys.resolve
43HH = qa.constants('H')['value']
44KK = qa.constants('K')['value']
45CC = qa.constants('C')['value']
47class solar_system_setjy:
48 def __init__(self):
49 self.models={}
52 def solar_system_fd (self, source_name, MJDs, frequencies, observatory, casalog=None):
53 '''
54 find flux density for solar system bodies:
55 Venus - Butler et al. 2001
56 Mars - Butler et al. 2012
57 Jupiter - Orton et al. 2012
58 Uranus - Orton & Hofstadter 2012 (modified ESA4)
59 Neptune - Orton & Hofstadter 2012 (modified ESA3)
60 Io - Butler et al. 2012
61 Europa - Butler et al. 2012
62 Ganymede - Butler et al. 2012
63 Titan - Gurwell et al. 2012
64 Callisto - Butler et al. 2012
65 Ceres - Mueller (private communication)
66 Juno - Butler et al. 2012
67 Pallas - Mueller (private communication)
68 Vesta - Mueller (private communication)
69 Lutetia - Mueller (private communication)
71 inputs:
72 source_name = source name string. example: "Venus"
73 MJDs = list of MJD times (day + fraction). example:
74 [ 56018.232, 56018.273 ]
75 must be sorted in ascending order.
76 frequencies = list of [start, stop] frequencies for
77 which to calculate the integrated model.
78 example:
79 [ [ 224.234567e9, 224.236567e9 ],
80 [ 224.236567e9, 224.238567e9 ] ]
81 observatory = observatory name string. example: "ALMA"
83 returned is a list, first element is the return status:
84 0 -> success
85 1 -> Error: unsupported body
86 2 -> Error: unsupported frequency or time for body
87 3 -> Error: Tb model file not found
88 4 -> Error: ephemeris table not found, or time out of range
89 (note - the case where the MJD times span two ephemeris
90 files is not supported)
91 5 -> Error: unknown observatory
92 second element is a list of flux densities, one per time and
93 frequency range, frequency changes fastest.
94 third element is list of uncertainties (if known; 0 if unknown),
95 one per time and frequency range, frequency changes fastest.
96 fourth element is a list of major axis, minor axis, and PAs in
97 asec and deg, one per MJD time.
98 fifth element is a list of CASA directions, one per MJD time.
100 bjb
101 nrao
102 spring/summer/fall 2012 + spring 2016
103 '''
105 RAD2ASEC = 2.0626480624710e5
106 AU = 1.4959787066e11
107 SUPPORTED_BODIES = [ 'Venus', 'Mars', 'Jupiter', 'Uranus', 'Neptune',
108 'Io', 'Europa', 'Ganymede', 'Callisto', 'Titan',
109 'Ceres', 'Juno', 'Pallas', 'Vesta', 'Hygeia' ]
110# those which have Tb (or f.d.) that is tabulated vs. time and frequency
111 TIME_VARIABLE_BODIES = [ 'Mars', 'Ceres', 'Pallas', 'Vesta', 'Lutetia' ]
112# those which are tabulations of flux density
113 MODEL_IS_FD_BODIES = [ 'Ceres', 'Pallas', 'Vesta', 'Lutetia' ]
115 capitalized_source_name = source_name.capitalize()
116 statuses = []
117 fds = []
118 dfds = []
119 Rhats = []
120 directions = []
122 #
123 # check that body is supported
124 #
125 if (not capitalized_source_name in SUPPORTED_BODIES):
126 for MJD in MJDs:
127 estatuses = []
128 efds = []
129 edfds = []
130 for frequency in frequencies:
131 estatuses.append(1)
132 efds.append(0)
133 edfds.append(0)
134 statuses.append(estatuses)
135 fds.append(efds)
136 dfds.append(edfds)
137 Rhats.append([0.0, 0.0, 0.0])
138 directions.append(_me.direction('J2000',0.0,0.0))
139 return [ statuses, fds, dfds, Rhats, directions ]
141 #
142 # check that observatory is known
143 #
144 if not observatory in _me.obslist():
145 for MJD in MJDs:
146 estatuses = []
147 efds = []
148 edfds = []
149 for frequency in frequencies:
150 estatuses.append(5)
151 efds.append(0)
152 edfds.append(0)
153 statuses.append(estatuses)
154 fds.append(efds)
155 dfds.append(edfds)
156 Rhats.append([0.0, 0.0, 0.0])
157 directions.append(_me.direction('J2000',0.0,0.0))
158 return [ statuses, fds, dfds, Rhats, directions ]
160 #
161 # before calculating the models be sure that we have the ephemeris
162 # information. otherwise don't waste our time calculating the model.
163 # only really important for mars, but do it for them all.
164 #
165 ephemeris_path = ctsys_resolve('ephemerides/JPL-Horizons')
166 ephemeris_file_list = os.listdir(ephemeris_path)
167 ephemeris_files = []
168 for ephemeris_file in ephemeris_file_list:
169 if (ephemeris_file.split('_')[0] == capitalized_source_name and 'J2000' in ephemeris_file):
170 ephemeris_files.append(ephemeris_file)
171 ephemeris_file_OK = False
172 #
173 # ephemeris tables have the following keywords:
174 # GeoDist, GeoLat, GeoLong, MJD0, NAME, VS_CREATE, VS_DATE, VS_TYPE,
175 # VS_VERSION, dMJD, earliest, latest, meanrad, obsloc, radii, rot_per
176 # and columns:
177 # MJD, RA, DEC, Rho (geodist), RadVel, NP_ang, NP_dist, DiskLong (Ob-long),
178 # DiskLat(Ob-lat), Sl_lon, Sl_lat, r (heliodist), rdot, phang
179 # Note by TT:
180 # The column names, Obs_lon and Obs_lat have been changed to DiskLong and
181 # DiskLat respectively to be consistent with what column names assumed for
182 # ephemeris tables by casacore's MeasComet class.
184 for ephemeris_file in ephemeris_files:
185 _tb.open(os.path.join(ephemeris_path,ephemeris_file))
186 table_source_name = _tb.getkeyword('NAME').capitalize()
187 if (table_source_name != capitalized_source_name):
188 continue
189 first_time = _tb.getkeyword('earliest')['m0']['value']
190 last_time = _tb.getkeyword('latest')['m0']['value']
191 if (first_time < MJDs[0] and last_time > MJDs[-1]):
192 ephemeris_file_OK = True
193 break
194 _tb.close()
195 #
196 # if we didn't find an ephemeris file, set the statuses and return.
197 #
198 if (not ephemeris_file_OK):
199 for MJD in MJDs:
200 estatuses = []
201 efds = []
202 edfds = []
203 for frequency in frequencies:
204 estatuses.append(4)
205 efds.append(0)
206 edfds.append(0)
207 statuses.append(estatuses)
208 fds.append(efds)
209 dfds.append(edfds)
210 Rhats.append([0.0, 0.0, 0.0])
211 directions.append(_me.direction('J2000',0.0,0.0))
212 return [ statuses, fds, dfds, Rhats, directions ]
214 Req = 1000.0 * (_tb.getkeyword('radii')['value'][0] + _tb.getkeyword('radii')['value'][1]) / 2
215 Rp = 1000.0 * _tb.getkeyword('radii')['value'][2]
216 times = _tb.getcol('MJD').tolist()
217 RAs = _tb.getcol('RA').tolist()
218 DECs = _tb.getcol('DEC').tolist()
219 distances = _tb.getcol('Rho').tolist()
220 RadVels = _tb.getcol('RadVel').tolist()
221 column_names = _tb.colnames()
222 if ('DiskLat' in column_names):
223 selats = _tb.getcol('DiskLat').tolist()
224 has_selats = 1
225 else:
226 has_selats = 0
227 selat = 0.0
228 if ('NP_ang' in column_names):
229 NPangs = _tb.getcol('NP_ang').tolist()
230 has_NPangs= 1
231 else:
232 has_NPangs = 0
233 NPang = 0.0
234 _tb.close()
235 MJD_shifted_frequencies = []
236 DDs = []
237 Rmeans = []
239 for ii in range(len(MJDs)):
240 MJD = MJDs[ii]
241 DDs.append(1.4959787066e11 * self.interpolate_list (times, distances, MJD)[1])
242 if (has_selats):
243 selat = self.interpolate_list (times, selats, MJD)[1]
244 if (selat == -999.0):
245 selat = 0.0
246 # apparent polar radius
247 Rpap = sqrt (Req*Req * sin(selat)**2.0 +
248 Rp*Rp * cos(selat)**2.0)
249 Rmean = sqrt (Rpap * Req)
250 Rmeans.append(Rmean)
251 #
252 # need to check that the convention for NP angle is the
253 # same as what is needed in the component list.
254 #
255 if (has_NPangs):
256 NPang = self.interpolate_list (times, NPangs, MJD)[1]
257 if (NPang == -999.0):
258 NPang = 0.0
259 Rhats.append([2*RAD2ASEC*Req/DDs[-1], 2*RAD2ASEC*Rpap/DDs[-1], NPang])
260 RA = self.interpolate_list (times, RAs, MJD)[1]
261 RAstr=str(RA)+'deg'
262 DEC = self.interpolate_list (times, DECs, MJD)[1]
263 DECstr=str(DEC)+'deg'
264 directions.append(_me.direction('J2000',RAstr,DECstr))
265 #
266 # now get the doppler shift
267 #
268 # NOTE: this is not exactly right, because it doesn't include the
269 # distance to the body in any of these calls. the distance will matter
270 # because it will change the line-of-sight vector from the observatory
271 # to the body, which will change the doppler shift. jeff thinks using
272 # the comet measures calls might fix this, but i haven't been able to
273 # figure them out yet. i thought i had it figured out, with the
274 # call to me.framecomet(), but that doesn't give the right answer,
275 # unfortunately. i spot-checked the error introduced because of this,
276 # and it looks to be of order 1 m/s for these bodies, so i'm not going
277 # to worry about it.
278 #
279 _me.doframe(_me.observatory(observatory))
280 _me.doframe(_me.epoch('utc',str(MJD)+'d'))
281 _me.doframe(directions[-1])
282 #
283 # instead of the call to me.doframe() in the line above, i thought the
284 # following call to me.framecomet() would be right, but it doesn't give
285 # the right answer :/.
286 # me.framecomet(ephemeris_file)
287 #
288 # RadVel is currently in AU/day. we want it in km/s.
289 #
290 RadVel = self.interpolate_list (times, RadVels, MJD)[1] * AU / 86400000.0
291 rv = _me.radialvelocity('geo',str(RadVel)+'km/s')
292 shifted_frequencies = []
293 for frequency in frequencies:
294 #
295 # the measure for radial velocity could be obtained via:
296 # me.measure(rv,'topo')['m0']['value']
297 # but what we really want is a frequency shift. i could do it by
298 # hand, but i'd rather do it with casa toolkit calls. unfortunately,
299 # it's a bit convoluted in casa...
300 #
301 newfreq0 = _me.tofrequency('topo',_me.todoppler('optical',_me.measure(rv,'topo')),_me.frequency('topo',str(frequency[0])+'Hz'))['m0']['value']
302 newfreq1 = _me.tofrequency('topo',_me.todoppler('optical',_me.measure(rv,'topo')),_me.frequency('topo',str(frequency[1])+'Hz'))['m0']['value']
303 #
304 # should check units to be sure frequencies are in Hz
305 #
306 # now, we want to calculate the model shifted in the opposite direction
307 # as the doppler shift, so take that into account.
308 #
309 delta_frequency0 = frequency[0] - newfreq0
310 newfreq0 = frequency[0] + delta_frequency0
311 delta_frequency1 = frequency[1] - newfreq1
312 newfreq1 = frequency[1] + delta_frequency1
313 shifted_frequencies.append([newfreq0,newfreq1])
314 average_delta_frequency = (delta_frequency0 + delta_frequency1)/2
315 #
316 # should we print this to the log?
317 #
318 # casalog.post( 'MJD, geo & topo velocities (km/s), and shift (MHz) = %7.1f %5.1f %5.1f %6.3f' % \
319 # (MJD, RadVel, me.measure(rv,'topo')['m0']['value']/1000, average_delta_frequency/1.0e6))
320 msg='MJD, geo & topo velocities (km/s), and shift (MHz) = %7.1f %5.1f %5.1f %6.3f' % \
321 (MJD, RadVel, _me.measure(rv,'topo')['m0']['value']/1000, average_delta_frequency/1.0e6)
322 casalog.post(msg, 'INFO2')
323 MJD_shifted_frequencies.append(shifted_frequencies)
324 # _me.done()
325 for ii in range(len(MJDs)):
326 shifted_frequencies = MJD_shifted_frequencies[ii]
327 if (capitalized_source_name in TIME_VARIABLE_BODIES):
328 [tstatuses,brightnesses,dbrightnesses] = self.brightness_time_int (capitalized_source_name,[MJDs[ii]], shifted_frequencies)
329 # modified by TT: take out an extra dimension (for times), to match the rest of the operation
330 tstatuses = tstatuses[0]
331 brightnesses = brightnesses[0]
332 dbrightnesses = dbrightnesses[0]
333 else:
334 tstatuses = []
335 brightnesses = []
336 dbrightnesses = []
337 for shifted_frequency in shifted_frequencies:
338 [status,brightness,dbrightness] = self.brightness_planet_int (capitalized_source_name, shifted_frequency)
339 tstatuses.append(status)
340 brightnesses.append(brightness)
341 dbrightnesses.append(dbrightness)
342 tfds = []
343 tdfds = []
344 for jj in range (len(tstatuses)):
345#
346# if a status of 6 was returned, then it's a body with a current
347# flux density model, but a model from the past which was not
348# flux density was used, so don't multiply by apparent size.
349#
350 if (tstatuses[jj] == 0 or tstatuses[jj] == 6):
351 flux_density = brightnesses[jj] # brigtnesses contains flux density already
352 if (capitalized_source_name not in MODEL_IS_FD_BODIES or tstatuses[jj] == 6):
353 flux_density *= 1.0e26 * pi * Rmeans[ii]*Rmeans[ii]/ (DDs[ii]*DDs[ii])
354 tstatuses[jj] = 0
355 #
356 # mean apparent planet radius, in arcseconds (used if we ever
357 # calculate the primary beam reduction)
358 #
359 psize = (Rmeans[ii] / DDs[ii]) * RAD2ASEC
360 #
361 # primary beam reduction factor (should call a function, but
362 # just set to 1.0 for now...
363 #
364 pbfactor = 1.0
365 flux_density *= pbfactor
366 tfds.append(flux_density)
367 tdfds.append(0.0)
368 else:
369 tfds.append(0.0)
370 tdfds.append(0.0)
371 statuses.append(tstatuses)
372 fds.append(tfds)
373 dfds.append(tdfds)
374 return [ statuses, fds, dfds, Rhats, directions ]
377 def brightness_time_int (self, source_name, MJDs, frequencies):
378 '''
379 Planck brightness for those bodies for which the data file
380 is a function of both frequency *and* time.
381 inputs:
382 source_name = source name (first character capitalized)
383 MJDs = list of MJD times
384 frequencies = list of [start, stop] frequencies for
385 which to calculate the integrated model.
386 example: [ [ 224.234567e9, 224.236567e9 ] ]
387 '''
389 # those bodies which are tabulations of flux density
390 MODEL_IS_FD_BODIES = [ 'Ceres', 'Pallas', 'Vesta', 'Lutetia' ]
391 HAS_OLD_MODEL_BODIES = [ 'Ceres', 'Pallas', 'Vesta' ]
393 statuses = []
394 Tbs = []
395 dTbs = []
396 model_data_path = ctsys_resolve('alma/SolarSystemModels')
397 if (source_name in MODEL_IS_FD_BODIES):
398 model_data_filename = os.path.join(model_data_path,source_name + '_fd_time.dat')
399 else:
400 model_data_filename = os.path.join(model_data_path,source_name + '_Tb_time.dat')
401 try:
402 ff = open(model_data_filename)
403 except:
404 for MJD in MJDs:
405 estatuses = []
406 eTbs = []
407 edTbs = []
408 for frequency in frequencies:
409 estatuses.append(3)
410 eTbs.append(0)
411 edTbs.append(0)
412 statuses.append(estatuses)
413 Tbs.append(eTbs)
414 dTbs.append(edTbs)
415 return [ statuses, Tbs, dTbs ]
417 # first line holds frequencies, like:
418 # 30.0 80.0 115.0 150.0 200.0 230.0 260.0 300.0 330.0 360.0 425.0 650.0 800.0 950.0 1000.0
419 line = ff.readline()[:-1]
420 fields = line.split()
421 freqs = []
422 for field in fields:
423 freqs.append(1.0e9*float(field))
424 # model output lines look like:
425 #2010 01 01 00 00 55197.00000 189.2 195.8 198.9 201.2 203.7 204.9 205.9 207.1 207.8 208.5 209.8 213.0 214.6 214.8 214.5
426 modelMJDs = []
427 modelTbs = []
428 for line in ff:
429 fields = line[:-1].split()
430 modelMJDs.append(float(fields[5]))
431 lTbs = []
432 for ii in range(len(freqs)):
433 lTbs.append(float(fields[6+ii]))
434 modelTbs.append(lTbs)
435 ff.close()
436 old_model_already_called = False
437 for MJD in MJDs:
438 #
439 # first, check if it's even a supported MJD (in the file)
440 #
441 if (MJD > modelMJDs[-1]):
442 estatuses = []
443 eTbs = []
444 edTbs = []
445 for frequency in frequencies:
446 estatuses.append(2)
447 eTbs.append(0.0)
448 edTbs.append(0.0)
449 elif (MJD < modelMJDs[0] and source_name in HAS_OLD_MODEL_BODIES):
450#
451# if the time is before the first time in the model data file,
452# then use the model that is not time variable (if this body
453# has one). set the status to 6 in that case, so that we know
454# in the calling function to multiply by the apparent size
455# to get true flux density.
456#
457 if (not old_model_already_called):
458#
459# only call this once - not for every time (since the old
460# models are not a function of time)
461#
462 tstatuses = []
463 tTbs = []
464 tdTbs = []
465 for frequency in frequencies:
466 [tstatus,tTb,tdTb] = self.brightness_planet_int (source_name, frequency)
467 tstatuses.append(tstatus)
468 tTbs.append(tTb)
469 tdTbs.append(tdTb)
470 old_model_already_called = True
471 estatuses = []
472 eTbs = []
473 edTbs = []
474 ii = 0
475 for frequency in frequencies:
476 if (tstatuses[ii]):
477 estatuses.append(tstatuses[ii])
478 else:
479 estatuses.append(6)
480 eTbs.append(tTbs[ii])
481 edTbs.append(tdTbs[ii])
482 ii += 1
483 else:
484 nind = self.nearest_index (modelMJDs, MJD)
485 mTbs = []
486 mfds = []
487 for ii in range(len(freqs)):
488 lMJD = []
489 lTb = []
490 for jj in range(nind-10, nind+10):
491 lMJD.append(modelMJDs[jj])
492 lTb.append(modelTbs[jj][ii])
493 if (source_name in MODEL_IS_FD_BODIES):
494 # don't know what Tb is, so just set to 0.0
495 mTbs.append(0.0)
496 mfds.append(self.interpolate_list(lMJD, lTb, MJD)[1])
497 else:
498 mTbs.append(self.interpolate_list(lMJD, lTb, MJD)[1])
499 #
500 # note: here, when we have the planck results, get a proper
501 # estimate of the background temperature.
502 #
503 # note also that we want to do this here because the integral
504 # needs to be done on the brightness, not on the brightness
505 # *temperature*.
506 #
507 Tbg = 2.725
508 mfds.append((2.0 * HH * freqs[ii]**3.0 / CC**2.0) * \
509 ((1.0 / (exp(HH * freqs[ii] / (KK * mTbs[-1])) - 1.0)) - \
510 (1.0 / (exp(HH * freqs[ii] / (KK * Tbg)) - 1.0))))
511 estatuses = []
512 eTbs = []
513 edTbs = []
514 for frequency in frequencies:
515 if (frequency[0] < freqs[0] or frequency[1] > freqs[-1]):
516 estatuses.append(2)
517 eTbs.append(0.0)
518 edTbs.append(0.0)
519 else:
520 [estatus, eTb, edTb] = self.integrate_Tb (freqs, mfds, frequency)
521 #
522 # should we print out the Tb we found? not sure. i have a
523 # vague recollection that crystal requested it, but i'm not
524 # sure if it's really needed. we'd have to back out the
525 # planck correction (along with the background), so it wouldn't
526 # be trivial.
527 #
528 estatuses.append(estatus)
529 eTbs.append(eTb)
530 edTbs.append(edTb)
531 statuses.append(estatuses)
532 Tbs.append(eTbs)
533 dTbs.append(edTbs)
534 return [statuses, Tbs, dTbs ]
537 def brightness_planet_int (self, source_name, frequency):
538 '''
539 brightness temperature for supported planets. integrates over
540 a tabulated model. inputs:
541 source_name = source name string
542 frequency = list of [start, stop] frequencies for
543 which to calculate the integrated model.
544 example: [ 224.234567e9, 224.236567e9 ]
545 '''
547 # those bodies which are tabulations of flux density. currently, for
548 # the non-time-variable ones, none of them are flux density. but
549 # leave this in just in case somewhere down the road we have bodies
550 # that we have a flux density model for that isn't time variable
551 # (evolved stars, for instance).
552 MODEL_IS_FD_BODIES = [ ]
554 if not source_name in self.models:
555 model_data_path = ctsys_resolve('alma/SolarSystemModels')
556 if (source_name in MODEL_IS_FD_BODIES):
557 model_data_filename = os.path.join(model_data_path,source_name + '_fd.dat')
558 else:
559 model_data_filename = os.path.join(model_data_path,source_name + '_Tb.dat')
560 try:
561 ff = open(model_data_filename)
562 except:
563 return [ 3, 0.0, 0.0 ]
564 fds = []
565 Tbs = []
566 freqs = []
567 for line in ff:
568 [freq,Tb] = line[:-1].split()
569 freqs.append(1.0e9*float(freq))
570 if (source_name in MODEL_IS_FD_BODIES):
571# don't know what Tb is, so just set to 0.0
572 Tbs.append(0.0)
573 fds.append(float(Tb))
574 else:
575 Tbs.append(float(Tb))
576#
577# note: here, when we have the planck results, get a proper
578# estimate of the background temperature.
579#
580# note also that we want to do this here because the integral
581# needs to be done on the brightness, not on the brightness
582# *temperature*.
583#
584 Tbg = 2.725
585 fds.append((2.0 * HH * freqs[-1]**3.0 / CC**2.0) * \
586 ((1.0 / (exp(HH * freqs[-1] / (KK * Tbs[-1])) - 1.0)) - \
587 (1.0 / (exp(HH * freqs[-1] / (KK * Tbg)) - 1.0))))
588 ff.close()
589# TT added. store them for future calls
590 srcn=source_name
591 self.models[srcn]={}
592 self.models[srcn]['fds']=fds
593 self.models[srcn]['freqs']=freqs
594 else:
595 #recover fds, freqs, instead of reading them in from the file
596 fds=self.models[source_name]['fds']
597 freqs=self.models[source_name]['freqs']
598 if (frequency[0] < freqs[0] or frequency[1] > freqs[-1]):
599 return [ 2, 0.0, 0.0 ]
600 else:
601 #
602 # should we print out the Tb we found? not sure. i have a
603 # vague recollection that crystal requested it, but i'm not
604 # sure if it's really needed. we'd have to back out the
605 # planck correction (along with the background), so it wouldn't
606 # be trivial. and here, we'd have to return a variable and
607 # work on that.
608 #
609 return self.integrate_Tb (freqs, fds, frequency)
612 def nearest_index (self, input_list, value):
613 """
614 find the index of the list input_list that is closest to value
615 """
617 ind = searchsorted(input_list, value)
618 ind = min(len(input_list)-1, ind)
619 ind = max(1, ind)
620 if value < (input_list[ind-1] + input_list[ind]) / 2.0:
621 ind = ind - 1
622 return ind
625 def interpolate_list (self, freqs, Tbs, frequency):
626 ind = self.nearest_index (freqs, frequency)
627 low = max(0,ind-5)
628 if (low == 0):
629 high = 11
630 else:
631 high = min(len(freqs),ind+5)
632 if (high == len(freqs)):
633 low = high - 11
634 #
635 # i wanted to put in a check for tabulated values that change
636 # derivative, since that confuses the interpolator. benign cases are
637 # fine, like radial velocity, but for the model Tbs, where there are
638 # sharp spectral lines, then the fitting won't be sensible when you're
639 # right at the center of the line, because the inflection is so severe.
640 # i thought if i just only took values that had the same derivative as
641 # the location where the desired value is that would work, but it
642 # doesn't :/. i'm either not doing it right or there's something
643 # deeper.
644 #
645 # if freqs[ind] < frequency:
646 # deriv = Tbs[ind+1] - Tbs[ind]
647 # else:
648 # deriv = Tbs[ind] - Tbs[ind-1]
649 # tTbs = []
650 # tfreqs = []
651 # for ii in range(low,high):
652 # nderiv = Tbs[ii+1] - Tbs[ii]
653 # if (nderiv >= 0.0 and deriv >= 0.0) or (nderiv < 0.0 and deriv < 0.0):
654 # tTbs.append(Tbs[ii])
655 # tfreqs.append(freqs[ii])
656 # aTbs = array(tTbs)
657 # afreqs = array(tfreqs)
658 aTbs = array(Tbs[low:high])
659 afreqs = array(freqs[low:high])
660 #
661 # cubic interpolation blows up near line centers (see above comment),
662 # so check that it doesn't fail completely (put it in a try/catch), and
663 # also that it's not a NaN and within the range of the tabulated values
664 #
665 range = max(aTbs) - min(aTbs)
666 try:
667 func = interp1d (afreqs, aTbs, kind='cubic')
668 if isnan(func(frequency)) or func(frequency) < min(aTbs)-range/2 or func(frequency) > max(aTbs)+range/2:
669 func = interp1d (afreqs, aTbs, kind='linear')
670 except:
671 func = interp1d (afreqs, aTbs, kind='linear')
672 #
673 # if it still failed, even with the linear interpolation, just take the
674 # nearest tabulated point.
675 #
676 if isnan(func(frequency)) or func(frequency) < min(aTbs)-range/2 or func(frequency) > max(aTbs)+range/2:
677 brightness = Tbs[ind]
678 else:
679 brightness = float(func(frequency))
680 return [ 0, brightness, 0.0 ]
683 def integrate_Tb (self, freqs, Tbs, frequency):
684 [status,low_Tb,low_dTb] = self.interpolate_list (freqs, Tbs, frequency[0])
685 low_index = self.nearest_index (freqs, frequency[0])
686 if (frequency[0] > freqs[low_index]):
687 low_index = low_index + 1
689 [status,hi_Tb,hi_dTb] = self.interpolate_list (freqs, Tbs, frequency[1])
690 hi_index = self.nearest_index (freqs, frequency[1])
691 if (frequency[1] < freqs[hi_index]):
692 hi_index = hi_index - 1
694 if (low_index > hi_index):
695 Tb = (frequency[1] - frequency[0]) * (low_Tb + hi_Tb) / 2
696 else:
697 Tb = (freqs[low_index] - frequency[0]) * (low_Tb + Tbs[low_index]) / 2 + \
698 (frequency[1] - freqs[hi_index]) * (hi_Tb + Tbs[hi_index]) / 2
699 ii = low_index
700 while (ii < hi_index):
701 Tb += (freqs[ii+1] - freqs[ii]) * (Tbs[ii+1] + Tbs[ii]) / 2
702 ii += 1
703 Tb /= (frequency[1] - frequency[0])
704 return [ 0, Tb, 0.0 ]