| Home | Trees | Indices | Help |
|
|---|
|
|
1 from __future__ import division
2 ## Automatically adapted for numpy.numarray Mar 05, 2007 by
3
4 import string
5 import glob
6 import re
7 import os
8 import warnings
9 import numpy as N
10 import pyfits
11
12 import spectrum
13 import units
14 import locations
15 from locations import irafconvert, _refTable
16 import planck
17 import wavetable
18 from tables import CompTable, GraphTable
19
20
21 #Flag to control verbosity
22 DEBUG = False
23
24 rootdir = locations.rootdir
25 datadir = locations.specdir
26 wavecat = locations.wavecat
27
28 #Constants to hold tables.
29 GRAPHTABLE= ''
30 GRAPHDICT = {}
31 COMPTABLE = ''
32 COMPDICT = {}
33 THERMTABLE = ''
34 THERMDICT = {}
35 HSTAREA = 45238.93416 # cm^2
36
38 global GRAPHTABLE, COMPTABLE, THERMTABLE, HSTAREA
39 # Component tables are defined here.
40
41 try:
42 GRAPHTABLE = _refTable(os.path.join('mtab','*_tmg.fits'))
43 COMPTABLE = _refTable(os.path.join('mtab','*_tmc.fits'))
44 except IOError, e:
45 GRAPHTABLE = None
46 COMPTABLE = None
47 warnings.warn("PYSYN_CDBS is undefined; No graph or component tables could be found; functionality will be SEVERELY crippled.",UserWarning)
48 try:
49 THERMTABLE = _refTable(os.path.join('mtab','*_tmt.fits'))
50 except IOError, e:
51 THERMTABLE = None
52 print "Warning: %s"%str(e)
53 print " No thermal calculations can be performed."
54
55 HSTAREA = 45238.93416 # cm^2
56
57 #Do this on import
58 _set_default_refdata()
59
60
63 """provide user access to global reference data.
64 Graph/comp/therm table names must be fully specified."""
65
66 global GRAPHTABLE, COMPTABLE, THERMTABLE, HSTAREA, GRAPHDICT, COMPDICT, THERMDICT
67
68 GRAPHDICT = {}
69 COMPDICT = {}
70 THERMDICT = {}
71
72 #Check for all None, which means reset
73 kwds=set([graphtable,comptable,thermtable,area])
74 if kwds == set([None]):
75 #then we should reset everything.
76 _set_default_refdata()
77 return
78
79 #Otherwise, check them all separately
80 if graphtable is not None:
81 GRAPHTABLE = irafconvert(graphtable)
82
83 if comptable is not None:
84 COMPTABLE = irafconvert(comptable)
85
86 if thermtable is not None:
87 THERMTABLE = irafconvert(thermtable)
88
89 #Area is a bit different:
90 if area is not None:
91 HSTAREA = area
92
93 #That's it.
94 return
95
97 """Collects & returns the current refdata as a dictionary"""
98 ans=dict(graphtable=GRAPHTABLE,
99 comptable=COMPTABLE,
100 thermtable=THERMTABLE,
101 area=HSTAREA)
102 return ans
103
105 """Prints the values settable by setref"""
106 refdata = getref()
107 for k, v in refdata.items():
108 print "%10s: %s"%(k,v)
109
110 CLEAR = 'clear'
111
113 ''' Class that handles the graph table, common to both optical and
114 thermal obsmodes.
115 '''
117 #Strip "band()" syntax if present
118 tmatch=re.search(r'band\((.*?)\)',obsmode,re.IGNORECASE)
119 if tmatch:
120 obsmode=tmatch.group(1)
121 self._obsmode = obsmode
122
123 if graphtable is None:
124 graphtable=GRAPHTABLE
125
126 self.area = HSTAREA
127
128 # For sensitivity calculations: 5.03411762e7 is hc in
129 # the appropriate units
130 self._constant = 5.03411762e7 * self.area
131 self.pardict={}
132
133 modes = obsmode.lower().split(',')
134 if '#' in obsmode:
135 self.modes=[]
136 for m in modes:
137 if '#' in m:
138 key,val=m.split('#')
139 self.pardict[key]=float(val)
140 self.modes.append("%s#"%key)
141 else:
142 self.modes.append(m)
143 else:
144 self.modes=modes
145
146 # gt = GraphTable(graphtable)
147 if graphtable in GRAPHDICT.keys():
148 gt = GRAPHDICT[graphtable]
149 else:
150 gt = GraphTable(graphtable)
151 GRAPHDICT[graphtable] = gt
152
153 self.gtname=graphtable
154
155 self.compnames,self.thcompnames = gt.GetComponentsFromGT(self.modes,1)
156
157 self.components = None #Will be filled by subclasses
158 self.pixscale = None
159
160 obm=self._obsmode.lower()
161
162 try:
163 self.binset = wavetable.wavetable[obm]
164 except KeyError,e:
165 #If zero candidates were found, that's ok.
166 pass
167 except ValueError,e:
168 #wavetable will raise a ValueError if the key was ambiguous
169 print "Warning, %s"%str(e)
170
171
174
175
178
180 files = []
181 for compname in compnames:
182 if compname not in [None, '', CLEAR]:
183 index = N.where(comptable.compnames == compname)
184 try:
185 iraffilename = comptable.filenames[index[0][0]]
186 filename = irafconvert(iraffilename)
187 files.append(filename.lstrip())
188 except IndexError:
189 raise IndexError("Can't find %s in comptable %s"%(compname,comptable.name))
190 else:
191 files.append(CLEAR)
192
193 return files
194
197
199 """ Duplicate synphot showfiles behavior"""
200 for name in self._throughput_filenames:
201 if name != 'clear':
202 print name
203
205 """ Return the binned waveset most appropriate for the obsmode,
206 as defined by the wavecat.dat file. """
207
208 if self.binset.startswith('('):
209 return self._computeBandwave(self.binset)
210 else:
211 return self._getBandwaveFomFile(self.binset)
212
214 (a,b,c,nwave) = self._computeQuadraticCoefficients(coeff)
215
216 result = N.zeros(shape=[nwave,], dtype=N.float64)
217
218 for i in range(nwave):
219 result[i] = ((a * i) + b) * i + c
220
221 return result
222
224
225 coefficients = (coeff[1:][:-1]).split(',')
226
227 c0 = float(coefficients[0])
228 c1 = float(coefficients[1])
229 c2 = (c1 - c0) / 1999.0 # arbitraily copied from synphot....
230 #In synphot.countrate/calcstep.x, it was NSPEC-1, where
231 #NSPEC was hardcoded to 2000 as the number of bins into
232 #which the wavelength set should be divided by default
233 c3 = c2
234 if len(coefficients) > 2:
235 c2 = float(coefficients[2])
236 c3 = c2
237 if len(coefficients) > 3:
238 c3 = float(coefficients[3])
239
240 nwave = int(2.0 * (c1 - c0) / (c3 + c2)) + 1
241
242 c = c0
243 b = c2
244 a = (c3 * c3 - c2 * c2) / (4.0 * (c1 - c0))
245
246 return (a,b,c,nwave)
247
249 name = irafconvert(filename)
250
251 fs = open(name, mode='r')
252 lines = fs.readlines()
253 fs.close()
254
255 tokens = []
256 for line in lines:
257 if not line.startswith('#'):
258 tokens.append(line)
259
260 return N.float_(tokens)
261
262
264
265 - def __init__(self, obsmode, method='HSTGraphTable',graphtable=None,
266 comptable=None, component_dict = {}):
267
268 if graphtable is None:
269 graphtable=GRAPHTABLE
270 if comptable is None:
271 comptable=COMPTABLE
272
273 BaseObservationMode.__init__(self, obsmode, method, graphtable)
274
275 # ct = CompTable(comptable)
276 if comptable in COMPDICT.keys():
277 ct = COMPDICT[comptable]
278 else:
279 ct = CompTable(comptable)
280 COMPDICT[comptable] = ct
281
282 self.ctname = comptable
283
284 self._throughput_filenames = self._getFileNames(ct, self.compnames)
285
286 self.components = self._getOpticalComponents(self._throughput_filenames,
287 component_dict)
288
290 components = []
291 for throughput_name in throughput_filenames:
292 if throughput_name.endswith('#]'):
293 barename,parkey=throughput_name.split('[')
294 parkey=parkey[:-2]
295 else:
296 parkey=None
297
298 if (throughput_name, self.pardict.get(parkey)) in component_dict.keys():
299 component = component_dict[(throughput_name, self.pardict.get(parkey))]
300 else:
301 component = _Component(throughput_name,
302 interpval=self.pardict.get(parkey))
303 component_dict[(throughput_name, self.pardict.get(parkey))] = component
304
305 if not component.isEmpty():
306 components.append(component)
307
308 return components
309
311 '''Calculate the sensitivity by combining the throughput curves
312 with hc/lambda to convert erg/cm^2/sec/Angstrom to counts/sec.
313 Multiplying this by the flux in erg/cm^2/sec/Angstrom will give
314 counts/sec/Angstrom'''
315 sensitivity = spectrum.TabularSpectralElement()
316
317 product = self._multiplyThroughputs()
318
319 sensitivity._wavetable = product.GetWaveSet()
320 sensitivity._throughputtable = product(sensitivity._wavetable) * \
321 sensitivity._wavetable * self._constant
322
323 return sensitivity
324
326 '''Throughput returns the TabularSpectralElement obtained by
327 multiplying the SpectralElement components together. Unitless'''
328 try:
329 throughput = spectrum.TabularSpectralElement()
330
331 product = self._multiplyThroughputs(0)
332
333 throughput._wavetable = product.GetWaveSet()
334 throughput._throughputtable = product(throughput._wavetable)
335 throughput.waveunits = product.waveunits
336 throughput.name='*'.join([str(x) for x in self.components])
337
338 ## throughput = throughput.resample(spectrum.default_waveset)
339
340 return throughput
341
342 except IndexError: # graph table is broken.
343 return None
344
345
347 product = self.components[index].throughput
348 if len(self.components) > index:
349 for component in self.components[index+1:]:
350 if component.throughput != None:
351 product = product * component.throughput
352 return product
353
354
356 try:
357 # delegate to subclass.
358 thom = _ThermalObservationMode(self._obsmode)
359 self.pixscale = thom.pixscale
360 return thom._getSpectrum()
361 except IndexError: # graph table is broken.
362 raise IndexError("Cannot calculate thermal spectrum; graphtable may be broken")
363
364
366
367 - def __init__(self, obsmode, method='HSTGraphTable',graphtable=None,
368 comptable=None, thermtable=None):
369
370 if graphtable is None:
371 graphtable = GRAPHTABLE
372 if comptable is None:
373 comptable = COMPTABLE
374 if thermtable is None:
375 thermtable = THERMTABLE
376
377
378 #The constructor of the parent class defines the self.thcompnames
379 BaseObservationMode.__init__(self, obsmode, method, graphtable)
380
381 #Check here to see if there are any.
382 if set(self.thcompnames).issubset(set(['clear',''])):
383 raise NotImplementedError("No thermal support provided for %s"%obsmode)
384
385 # ct = CompTable(comptable)
386 if comptable in COMPDICT.keys():
387 ct = COMPDICT[comptable]
388 else:
389 ct = CompTable(comptable)
390 COMPDICT[comptable] = ct
391
392 self.ctname=comptable
393
394 throughput_filenames = self._getFileNames(ct, self.compnames)
395
396 # thct = CompTable(thermtable)
397 if thermtable in THERMDICT.keys():
398 thct = THERMDICT[thermtable]
399 else:
400 thct = CompTable(thermtable)
401 THERMDICT[thermtable] = thct
402
403 self.thname = thermtable
404
405 thermal_filenames = self._getFileNames(thct, self.thcompnames)
406
407 self.components = self._getThermalComponents(throughput_filenames, \
408 thermal_filenames)
409
410 self.pixscale = self._getPixelScale()
411 self.name = obsmode+" (thermal)"
412
414 obsmode = self._obsmode.split(',')
415 obsmode = str(obsmode[0]) + ',' + str(obsmode[1])
416
417 fname=os.path.join(locations.specdir,'detectors.dat')
418 fs = open(fname,mode='r')
419 lines = fs.readlines()
420 fs.close()
421
422 regx = re.compile(r'\S+', re.IGNORECASE)
423 for line in lines:
424 try:
425 tokens = regx.findall(line)
426 if tokens[0] == obsmode:
427 break
428 except Exception, e:
429 raise ValueError("Error processing %s: %s"%(fname,str(e)))
430
431 return float(tokens[1])
432
434 components = []
435 for i in range(len(throughput_filenames)):
436 throughput_name = throughput_filenames[i]
437 thermal_name = thermal_filenames[i]
438 if throughput_name.endswith('#]'):
439 barename,parkey=throughput_name.split('[')
440 parkey=parkey[:-2]
441 else:
442 parkey=None
443
444 component = _ThermalComponent(throughput_name, thermal_name, \
445 interpval=self.pardict.get(parkey))
446 if not component.isEmpty():
447 components.append(component)
448
449 return components
450
452 ''' Overrides base class in order to deal with opaque components.
453 '''
454 index = 0
455 for component in self.components:
456 if component.throughput != None:
457 break
458 index += 1
459
460 return BaseObservationMode._multiplyThroughputs(self, index)
461
463 wave=self._getWavesetIntersection()
464 sp = spectrum.ArraySourceSpectrum(wave=wave,
465 flux=N.zeros(shape=wave.shape,dtype=N.float64),
466 waveunits='angstrom',
467 fluxunits='photlam',
468 name="%s %s"%(self.name,'ThermalSpectrum'))
469
470
471 minw = sp._wavetable[0]
472 maxw = sp._wavetable[-1]
473
474 for component in self.components:
475 # transmissive section
476 if component.throughput != None:
477 sp = sp * component.throughput
478
479 # sp = spectrum.trimSpectrum(sp, minw, maxw)
480
481 # thermal section
482 if component.emissivity != None:
483 bb = self._bb(sp.GetWaveSet(), component.emissivity.temperature)
484
485 sp_comp = component.emissivity.beamFillFactor * bb * \
486 component.emissivity
487
488 sp = sp + sp_comp
489
490 sp = spectrum.trimSpectrum(sp, minw, maxw)
491
492 return sp
493
495 minw = spectrum.default_waveset[0]
496 maxw = spectrum.default_waveset[-1]
497
498 for component in self.components[1:]:
499 if component.emissivity != None:
500 wave = component.emissivity.GetWaveSet()
501
502 minw = max(minw, wave[0])
503 maxw = min(maxw, wave[-1])
504
505 result = self._mergeEmissivityWavesets()
506
507 result = N.compress(result > minw, result)
508 result = N.compress(result < maxw, result)
509
510 # intersection with vega spectrum (why???)
511 vegasp = spectrum.TabularSourceSpectrum(locations.VegaFile)
512 vegaws = vegasp.GetWaveSet()
513 result = N.compress(result > vegaws[0], result)
514 result = N.compress(result < vegaws[-1], result)
515
516 return result
517
519 index = 1
520
521 for component in self.components:
522 emissivity = component.emissivity
523 if emissivity == None:
524 index = index + 1
525 else:
526 result = emissivity.GetWaveSet()
527 break;
528
529 for component in self.components[index:]:
530 if component.emissivity != None:
531 result = spectrum.MergeWaveSets(result, \
532 component.emissivity.GetWaveSet())
533 return result
534
536 sp = spectrum.ArraySourceSpectrum(wave=wave,
537 flux=planck.bb_photlam_arcsec(wave, temperature),
538 name='planck bb_photlam_arcsec')
539 return sp
540
541
544 self.throughput_name = throughput_name
545
546 self._empty = True
547
548 self.throughput = self._buildThroughput(throughput_name, interpval)
549 if self.throughput is not None:
550 self.waveunits = self.throughput.waveunits
551
553 return str(self.throughput)
554
556 if name != CLEAR:
557 if interpval is None:
558 self._empty = False
559 return spectrum.TabularSpectralElement(name)
560 else:
561 self._empty = False
562 return spectrum.InterpolatedSpectralElement(name, interpval)
563 else:
564 return None
565
568
569
571
573 self.throughput_name = throughput_name
574 self.thermal_name = thermal_name
575
576 self._empty = True
577
578 self.throughput = self._buildThroughput(throughput_name, interpval)
579
580 if thermal_name != CLEAR:
581 self._empty = False
582 self.emissivity = spectrum.ThermalSpectralElement(thermal_name)
583 else:
584 self.emissivity = None
585
| Home | Trees | Indices | Help |
|
|---|
| Generated by Epydoc 3.0.1 on Mon Aug 22 14:37:49 2011 | http://epydoc.sourceforge.net |