pysynphot.units

1 from __future__ import division 2 """ 3 Units class hierarchy: is used to manage both wavelength and flux 4 unit conversions 5 6 Warning: vegamag unit conversions require spectrum and locations => circular 7 imports. 8 """ 9 10 import math 11 import numpy as N 12 import locations, spectrum #Circular import 13 14 import observationmode #needed for HSTAREA 15 #cannot just import the constant because it won't get updated 16 #when the setref() function is used to change it. 17 18 19 C = 2.99792458e18 # speed of light in Angstrom/sec 20 H = 6.62620E-27 # Planck's constant 21 HC = H * C 22 ABZERO = -48.60 # magnitude zero points 23 STZERO = -21.10 24 25 26 27

28 -def Units(uname):

29 """This needs to be a factory function in order to return an object 30 of the correct subclass.""" 31 if isinstance(uname,BaseUnit): 32 return uname 33 else: 34 try: 35 if issubclass(uname,BaseUnit): 36 return uname() 37 except TypeError: 38 39 try: 40 return factory(uname) 41 except KeyError: 42 if uname == str(None): 43 return None 44 else: 45 raise ValueError("Unknown units %s"%uname)

46 47 #......................................................................

48 -def ismatch(a,b):

49 """Method to allow smart comparisons between classes, instances, 50 and string representations of units and give the right answer.""" 51 match=False 52 #Try the easy case 53 if a == b: 54 return True 55 else: 56 #Try isinstance in both orders 57 try: 58 if isinstance(a,b): 59 return True 60 except TypeError: 61 try: 62 if isinstance(b,a): 63 return True 64 except TypeError: 65 #Try isinstance(a, type(b)) in both orders 66 try: 67 if isinstance(a,type(b)): 68 return True 69 except TypeError: 70 try: 71 if isinstance(b,type(a)): 72 return True 73 except TypeError: 74 #Try the string representation 75 if str(a).lower() == str(b).lower(): 76 return True 77 else: 78 return False

79 80 #...................................................................... 81 #Base classes 82

83 -class BaseUnit(object):

84 """ Base class for all units; defines UI"""

85 - def __init__(self,uname):

86 self.Dispatch=None 87 self.name=uname

88

89 - def __str__(self):

90 return self.name

91

92 - def Convert(self,wave,flux,target_units):

93 #This signature is appropriate for fluxes, not waves 94 try: 95 return self.Dispatch[target_units.lower()](wave,flux) 96 except KeyError: 97 raise TypeError("%s cannot be converted to %s"%(self.name, 98 target_units))

99

100 -class WaveUnits(BaseUnit):

101 """All WaveUnits know how to convert themselves to Angstroms"""

102 - def __init__(self):

103 self.name=None 104 self.isFlux = False 105 self.Dispatch = {'angstrom' : self.ToAngstrom}

106

107 - def Convert(self,wave,target_units):

108 """WaveUnits only need a wavelength table to do a conversion.""" 109 try: 110 return self.Dispatch[target_units.lower()](wave) 111 except KeyError: 112 raise TypeError("%s cannot be converted to %s"%(self.name, 113 target_units))

114

115 - def ToAngstrom(self,wave):

116 raise NotImplementedError("Required method ToAngstrom not yet implemented")

117

118 -class FluxUnits(BaseUnit):

119 """All FluxUnits know how to convert themselves to Photlam"""

120 - def __init__(self):

121 self.isFlux = True 122 self.isMag = False 123 self.isDensity = True #False for counts and obmag 124 self.name=None 125 self.Dispatch = {'photlam':self.ToPhotlam} 126 self.nativewave = Angstrom

127 #self.StdSpectrum = None 128 #...StdSpectrum is a placeholder. Actual values for the attributes 129 # be defined in the renorm.py module; they can't be done here 130 # because of a circular import problem. If you add a new fluxunit 131 # in this file, you must define its StdSpectrum in renorm.py. 132

133 - def Convert(self,wave,flux,target_units):

134 """FluxUnits need both wavelength and flux tables to do a unit conversion.""" 135 try: 136 return self.Dispatch[target_units](wave,flux) 137 except KeyError: 138 raise TypeError("%s is not a valid flux unit"%(target_units))

139

140 - def ToPhotlam(self,wave,flux):

141 raise NotImplementedError("Required method ToPhotlam not yet implemented")

142 143

144 -class LogFluxUnits(FluxUnits):

145 """Base class for magnitudes, which often require special handling"""

146 - def __init__(self):

147 FluxUnits.__init__(self) 148 self.isMag=True 149 self.linunit=None 150 self.zeropoint=None

151 152 #............................................................. 153 # Internal wavelength units are Angstroms, so it is smarter than the others

154 -class Angstrom(WaveUnits):

155 - def __init__(self):

156 WaveUnits.__init__(self) 157 self.name = 'angstrom' 158 self.Dispatch = {'angstrom' : self.ToAngstrom, 159 'angstroms' : self.ToAngstrom, 160 'nm': self.ToNm, 161 'micron': self.ToMicron, 162 'microns': self.ToMicron, 163 '1/um': self.ToInverseMicron, 164 'inversemicron': self.ToInverseMicron, 165 'inversemicrons': self.ToInverseMicron, 166 'mm': self.ToMm, 167 'cm': self.ToCm, 168 'm': self.ToMeter, 169 'hz': self.ToHz}

170 171

172 - def ToAngstrom(self, wave):

173 return wave

174

175 - def ToNm(self, wave):

176 return wave / 10.0

177

178 - def ToMicron(self, wave):

179 return wave * 1.0e-4

180

181 - def ToInverseMicron(self, wave):

182 return 1.0e4/wave

183

184 - def ToMm(self, wave):

185 return wave * 1.0e-7

186

187 - def ToCm(self, wave):

188 return wave * 1.0e-8

189

190 - def ToMeter(self, wave):

191 return wave * 1.0e-10

192

193 - def ToHz(self, wave):

194 return C / wave

195 196 #........................................................................ 197 # Internal flux units = Photlam, so it is smarter than the others 198

199 -class Photlam(FluxUnits):

200 ''' photlam = photons cm^-2 s^-1 Ang^-1)'''

201 - def __init__(self):

202 FluxUnits.__init__(self) 203 self.name = 'photlam' 204 self.Dispatch = {'flam': self.ToFlam, 205 'fnu': self.ToFnu, 206 'photlam': self.ToPhotlam, 207 'photnu': self.ToPhotnu, 208 'jy': self.ToJy, 209 'mjy': self.TomJy, 210 'mujy': self.TomuJy, 211 'microjy': self.TomuJy, 212 'ujy': self.TomuJy, 213 'njy': self.TonJy, 214 'nanojy': self.TonJy, 215 'abmag': self.ToABMag, 216 'stmag': self.ToSTMag, 217 'obmag': self.ToOBMag, 218 'vegamag': self.ToVegaMag, 219 'counts': self.ToCounts, 220 'counts': self.ToCounts} 221 222 self.nativewave = Angstrom

223 224

225 - def unitResponse(self,band):

226 """Put a flat spectrum of 1 photlam through this band, & integrate""" 227 #sumfilt(wave,0,band) 228 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 229 total = band.trapezoidIntegration(band.wave,band.throughput) 230 return 1.0/total

231 232

233 - def ToFlam(self, wave, flux):

234 return HC * flux / wave

235

236 - def ToFnu(self, wave, flux):

237 return H * flux * wave

238

239 - def ToPhotlam(self, wave, flux):

240 return flux.copy() # No conversion, just copy the array.

241

242 - def ToPhotnu(self, wave, flux):

243 return flux * wave * wave / C

244

245 - def ToJy(self, wave, flux):

246 return 1.0e+23 * H * flux * wave

247

248 - def TomJy(self, wave, flux):

249 return 1.0e+26 * H * flux * wave

250

251 - def TomuJy(self, wave, flux):

252 return 1.0e+29 * H * flux * wave

253

254 - def TonJy(self, wave, flux):

255 return 1.0e+32 * H * flux * wave

256

257 - def ToABMag(self, wave, flux):

258 arg = H * flux * wave 259 return -1.085736 * N.log(arg) + ABZERO

260

261 - def ToSTMag(self, wave, flux):

262 arg = H * C* flux / wave 263 return -1.085736 * N.log(arg) + STZERO

264

265 - def ToOBMag(self, wave, flux):

266 dw = _getDeltaWave(wave) 267 arg = flux * dw * observationmode.HSTAREA 268 return -1.085736 * N.log(arg)

269

270 - def ToVegaMag(self, wave, flux):

271 272 resampled = spectrum.Vega.resample(wave) 273 normalized = flux / resampled._fluxtable 274 return -2.5 * N.log10(normalized)

275

276 - def ToCounts(self, wave, flux):

277 return flux * _getDeltaWave(wave) * observationmode.HSTAREA

278 279 280 #................................................................ 281 #Other wavelength units 282

283 -class Hz(WaveUnits):

284 - def __init__(self):

285 WaveUnits.__init__(self) 286 self.name='hz'

287

288 - def ToAngstrom(self, wave):

289 return C / wave

290

291 -class InverseMicron(WaveUnits):

292 - def __init__(self):

293 WaveUnits.__init__(self) 294 self.name = '1/um'

295

296 - def ToAngstrom(self, wave):

297 return 1.0e4/wave

298

299 -class _MetricWavelength(WaveUnits):

300 """ Encapsulates some easy unit-conversion machinery. Angstrom 301 is not subclassed from here because it needs to be especially smart in 302 other ways. (Although multiple inheritence might be possible.)"""

303 - def ToAngstrom(self,wave):

304 return wave*self.factor

305

306 -class Nm(_MetricWavelength):

307 - def __init__(self):

308 _MetricWavelength.__init__(self) 309 self.name = 'nm' 310 self.factor = 10.0

311

312 -class Micron(_MetricWavelength):

313 - def __init__(self):

314 _MetricWavelength.__init__(self) 315 self.name = 'micron' 316 self.factor = 1.0e4

317

318 -class Mm(_MetricWavelength):

319 - def __init__(self):

320 _MetricWavelength.__init__(self) 321 self.name = 'mm' 322 self.factor = 1.0e7

323

324 -class Cm(_MetricWavelength):

325 - def __init__(self):

326 _MetricWavelength.__init__(self) 327 self.name = 'cm' 328 self.factor = 1.0e8

329

330 -class Meter(_MetricWavelength):

331 - def __init__(self):

332 _MetricWavelength.__init__(self) 333 self.name = 'm' 334 self.factor = 1.0e10

335 336 337 #................................................................ 338 #Other flux units 339

340 -class Flam(FluxUnits):

341 ''' flam = erg cm^-2 s^-1 Ang^-1''' 342

343 - def __init__(self):

344 FluxUnits.__init__(self) 345 self.name='flam' 346 self.Dispatch = {'photlam':self.ToPhotlam} 347 self.nativewave = Angstrom

348

349 - def ToPhotlam(self, wave, flux):

350 return flux * wave / HC

351 352

353 - def unitResponse(self,band):

354 #sumfilt(wave,1,band) 355 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 356 wave=band.wave 357 total = band.trapezoidIntegration(wave,band.throughput*wave) 358 modtot = total / (H*C) 359 return 1.0/modtot

360 361

362 -class Photnu(FluxUnits):

363 ''' photnu = photon cm^-2 s^-1 Hz^-1'''

364 - def __init__(self):

365 FluxUnits.__init__(self) 366 self.name = 'photnu' 367 self.nativewave = Hz

368

369 - def ToPhotlam(self, wave, flux):

370 return C * flux / (wave * wave)

371 372

373 - def unitResponse(self,band):

374 #sumfilt(wave,-2,band) 375 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 376 wave=band.wave 377 total = band.trapezoidIntegration(wave,band.throughput/(wave*wave)) 378 modtot = total/C 379 return 1.0/modtot

380 381

382 -class Fnu(FluxUnits):

383 ''' fnu = erg cm^-2 s^-1 Hz^-1'''

384 - def __init__(self):

385 FluxUnits.__init__(self) 386 self.name = 'fnu' 387 self.nativewave = Hz

388

389 - def ToPhotlam(self, wave, flux):

390 return flux /wave / H

391 392

393 - def unitResponse(self,band):

394 #sumfilt(wave,-1,band) 395 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 396 wave=band.wave 397 total = band.trapezoidIntegration(wave,band.throughput/wave) 398 modtot = total/H 399 return 1.0/modtot

400

401 -class Jy(FluxUnits):

402 ''' jy = 10^-23 erg cm^-2 s^-1 Hz^-1'''

403 - def __init__(self):

404 FluxUnits.__init__(self) 405 self.name = 'jy' 406 self.nativewave = Hz

407

408 - def ToPhotlam(self, wave, flux):

409 return flux / wave * (1.0e-23 / H)

410 411

412 - def unitResponse(self,band):

413 #sumfilt(wave,-1,band) 414 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 415 wave=band.wave 416 total = band.trapezoidIntegration(wave,band.throughput/wave) 417 modtot = total * (1.0e-23/H) 418 return 1.0/modtot

419

420 -class mJy(FluxUnits):

421 ''' mjy = 10^-26 erg cm^-2 s^-1 Hz^-1'''

422 - def __init__(self):

423 FluxUnits.__init__(self) 424 self.name = 'mjy' 425 self.nativewave = Hz

426

427 - def ToPhotlam(self, wave, flux):

428 return flux / wave * (1.0e-26 / H)

429 430

431 - def unitResponse(self,band):

432 #sumfilt(wave,-1,band) 433 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 434 wave=band.wave 435 total = band.trapezoidIntegration(wave,band.throughput/wave) 436 modtot = total * (1.0e-26/H) 437 return 1.0/modtot

438

439 -class muJy(FluxUnits): # New

440 ''' mujy = 10^-29 erg cm^-2 s^-1 Hz^-1'''

441 - def __init__(self):

442 FluxUnits.__init__(self) 443 self.name = 'mujy' 444 self.nativewave = Hz

445

446 - def ToPhotlam(self, wave, flux):

447 return flux / wave * (1.0e-29 / H)

448

449 - def unitResponse(self,band):

450 wave=band.wave 451 total = band.trapezoidIntegration(wave,band.throughput/wave) 452 modtot = total * (1.0e-29/H) 453 return 1.0/modtot

454

455 -class nJy(FluxUnits): # New

456 ''' njy = 10^-32 erg cm^-2 s^-1 Hz^-1'''

457 - def __init__(self):

458 FluxUnits.__init__(self) 459 self.name = 'njy' 460 self.nativewave = Hz

461

462 - def ToPhotlam(self, wave, flux):

463 return flux / wave * (1.0e-32 / H)

464

465 - def unitResponse(self,band):

466 wave=band.wave 467 total = band.trapezoidIntegration(wave,band.throughput/wave) 468 modtot = total * (1.0e-32/H) 469 return 1.0/modtot

470

471 -class ABMag(LogFluxUnits):

472 - def __init__(self):

473 LogFluxUnits.__init__(self) 474 self.name = 'abmag' 475 self.linunit = Fnu() 476 self.zeropoint = ABZERO

477 478

479 - def ToPhotlam(self, wave, flux):

480 return 1.0 / (H * wave) * 10.0**(-0.4 * (flux - ABZERO))

481

482 - def unitResponse(self,band):

483 #sumfilt(wave,-1,band) 484 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 485 wave=band.wave 486 total = band.trapezoidIntegration(wave,band.throughput/wave) 487 modtot = total/H 488 return 2.5*math.log10(modtot) + ABZERO

489

490 -class STMag(LogFluxUnits):

491 - def __init__(self):

492 LogFluxUnits.__init__(self) 493 self.name = 'stmag' 494 self.linunit = Flam() 495 self.zeropoint = STZERO

496 497

498 - def ToPhotlam(self, wave, flux):

499 return wave / H / C * 10.0**(-0.4 * (flux - STZERO))

500

501 - def unitResponse(self,band):

502 #sumfilt(wave,1,band) 503 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 504 wave=band.wave 505 total = band.trapezoidIntegration(wave,band.throughput*wave) 506 modtot = total/(H*C) 507 return 2.5*math.log10(modtot) + STZERO

508

509 -class OBMag(LogFluxUnits):

510 - def __init__(self):

511 LogFluxUnits.__init__(self) 512 self.name = 'obmag' 513 self.linunit = Counts() 514 self.zeropoint = 0.0 515 self.isDensity = False

516

517 - def ToPhotlam(self, wave, flux):

518 dw = _getDeltaWave(wave) 519 return 10.0**(-0.4 * flux) / (dw * observationmode.HSTAREA)

520

521 - def unitResponse(self,band):

522 #sum = asumr(band,nwave) 523 total = band.throughput.sum() 524 return 2.5*math.log10(total)

525

526 -class VegaMag(LogFluxUnits):

527 - def __init__(self):

528 LogFluxUnits.__init__(self) 529 self.name = 'vegamag' 530 self.vegaspec = spectrum.Vega

531

532 - def ToPhotlam(self, wave, flux):

533 resampled = self.vegaspec.resample(wave) 534 return resampled.flux * 10.0**(-0.4 * flux)

535

536 - def unitResponse(self,band):

537 sp=band*self.vegaspec 538 total=sp.integrate() 539 return 2.5*math.log10(total)

540

541 -class Counts(FluxUnits):

542 - def __init__(self):

543 FluxUnits.__init__(self) 544 self.name = 'counts' 545 self.isDensity = False

546

547 - def ToPhotlam(self, wave, flux):

548 return flux / (_getDeltaWave(wave) * observationmode.HSTAREA)

549 550

551 - def unitResponse(self,band):

552 #sum = asumr(band,nwave) 553 total = band.throughput.sum() 554 return 1.0/total

555 556 557 ################ Factory for Units subclasses. ##################### 558 559

560 -def factory(uname, *args, **kwargs):

561 unitsClasses = {'flam' : Flam, 562 'fnu' : Fnu, 563 'photlam' : Photlam, 564 'photnu' : Photnu, 565 'jy' : Jy, 566 'mjy' : mJy, 567 'mujy' : muJy, 568 'microjy' : muJy, 569 'ujy' : muJy, 570 'njy' : nJy, 571 'nanojy' : nJy, 572 'abmag' : ABMag, 573 'stmag' : STMag, 574 'obmag' : OBMag, 575 'vegamag' : VegaMag, 576 'counts' : Counts, 577 'count' : Counts, 578 'angstrom' : Angstrom, 579 'angstroms' : Angstrom, 580 'nm' : Nm, 581 'micron' : Micron, 582 'microns' : Micron, 583 'um' : Micron, 584 'inversemicron': InverseMicron, 585 'inversemicrons': InverseMicron, 586 '1/um' : InverseMicron, 587 'mm' : Mm, 588 'cm' : Cm, 589 'm' : Meter, 590 'meter' : Meter, 591 'hz' : Hz} 592 593 key=uname.lower() 594 ans= unitsClasses[key]() 595 return ans

596

597 -def _getDeltaWave(wave):

598 """ Compute delta wavelngth for an array of wavelengths. 599 If we had a WaveTable class, this function would be a method 600 on that class: possible refactoring candidate. """ 601 602 last = wave.shape[0]-1 603 604 hold1 = N.empty(shape=wave.shape, dtype=N.float64) 605 hold2 = N.empty(shape=wave.shape, dtype=N.float64) 606 607 hold1[1::] = wave[0:last] 608 hold2[0:last] = wave[1::] 609 610 delta = (hold2 - hold1) / 2.0 611 612 delta[0] = wave[1] - wave[0] 613 delta[last] = wave[last] - wave[last-1] 614 615 return delta

616

Source Code for Module pysynphot.units