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 if hasattr(wave,'copy'): 174 return wave.copy() # to avoid writing over any internal wave objects 175 else: 176 return wave # probably a scalar

177

178 - def ToNm(self, wave):

179 return wave / 10.0

180

181 - def ToMicron(self, wave):

182 return wave * 1.0e-4

183

184 - def ToInverseMicron(self, wave):

185 return 1.0e4/wave

186

187 - def ToMm(self, wave):

188 return wave * 1.0e-7

189

190 - def ToCm(self, wave):

191 return wave * 1.0e-8

192

193 - def ToMeter(self, wave):

194 return wave * 1.0e-10

195

196 - def ToHz(self, wave):

197 return C / wave

198 199 #........................................................................ 200 # Internal flux units = Photlam, so it is smarter than the others 201

202 -class Photlam(FluxUnits):

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

204 - def __init__(self):

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

226 227

228 - def unitResponse(self,band):

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

234 235

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

237 return HC * flux / wave

238

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

240 return H * flux * wave

241

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

243 if hasattr(flux,'copy'): 244 return flux.copy() # No conversion, just copy the array. 245 else: 246 return flux # probably a scalar

247

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

249 return flux * wave * wave / C

250

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

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

253

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

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

256

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

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

259

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

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

262

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

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

266

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

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

270

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

272 dw = _getDeltaWave(wave) 273 arg = flux * dw * observationmode.HSTAREA 274 return -1.085736 * N.log(arg)

275

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

277 278 resampled = spectrum.Vega.resample(wave) 279 normalized = flux / resampled._fluxtable 280 return -2.5 * N.log10(normalized)

281

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

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

284 285 286 #................................................................ 287 #Other wavelength units 288

289 -class Hz(WaveUnits):

290 - def __init__(self):

291 WaveUnits.__init__(self) 292 self.name='hz'

293

294 - def ToAngstrom(self, wave):

295 return C / wave

296

297 -class InverseMicron(WaveUnits):

298 - def __init__(self):

299 WaveUnits.__init__(self) 300 self.name = '1/um'

301

302 - def ToAngstrom(self, wave):

303 return 1.0e4/wave

304

305 -class _MetricWavelength(WaveUnits):

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

309 - def ToAngstrom(self,wave):

310 return wave*self.factor

311

312 -class Nm(_MetricWavelength):

313 - def __init__(self):

314 _MetricWavelength.__init__(self) 315 self.name = 'nm' 316 self.factor = 10.0

317

318 -class Micron(_MetricWavelength):

319 - def __init__(self):

320 _MetricWavelength.__init__(self) 321 self.name = 'micron' 322 self.factor = 1.0e4

323

324 -class Mm(_MetricWavelength):

325 - def __init__(self):

326 _MetricWavelength.__init__(self) 327 self.name = 'mm' 328 self.factor = 1.0e7

329

330 -class Cm(_MetricWavelength):

331 - def __init__(self):

332 _MetricWavelength.__init__(self) 333 self.name = 'cm' 334 self.factor = 1.0e8

335

336 -class Meter(_MetricWavelength):

337 - def __init__(self):

338 _MetricWavelength.__init__(self) 339 self.name = 'm' 340 self.factor = 1.0e10

341 342 343 #................................................................ 344 #Other flux units 345

346 -class Flam(FluxUnits):

347 ''' flam = erg cm^-2 s^-1 Ang^-1''' 348

349 - def __init__(self):

350 FluxUnits.__init__(self) 351 self.name='flam' 352 self.Dispatch = {'photlam':self.ToPhotlam} 353 self.nativewave = Angstrom

354

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

356 return flux * wave / HC

357 358

359 - def unitResponse(self,band):

360 #sumfilt(wave,1,band) 361 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 362 wave=band.wave 363 total = band.trapezoidIntegration(wave,band.throughput*wave) 364 modtot = total / (H*C) 365 return 1.0/modtot

366 367

368 -class Photnu(FluxUnits):

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

370 - def __init__(self):

371 FluxUnits.__init__(self) 372 self.name = 'photnu' 373 self.nativewave = Hz

374

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

376 return C * flux / (wave * wave)

377 378

379 - def unitResponse(self,band):

380 #sumfilt(wave,-2,band) 381 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 382 wave=band.wave 383 total = band.trapezoidIntegration(wave,band.throughput/(wave*wave)) 384 modtot = total/C 385 return 1.0/modtot

386 387

388 -class Fnu(FluxUnits):

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

390 - def __init__(self):

391 FluxUnits.__init__(self) 392 self.name = 'fnu' 393 self.nativewave = Hz

394

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

396 return flux /wave / H

397 398

399 - def unitResponse(self,band):

400 #sumfilt(wave,-1,band) 401 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 402 wave=band.wave 403 total = band.trapezoidIntegration(wave,band.throughput/wave) 404 modtot = total/H 405 return 1.0/modtot

406

407 -class Jy(FluxUnits):

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

409 - def __init__(self):

410 FluxUnits.__init__(self) 411 self.name = 'jy' 412 self.nativewave = Hz

413

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

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

416 417

418 - def unitResponse(self,band):

419 #sumfilt(wave,-1,band) 420 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 421 wave=band.wave 422 total = band.trapezoidIntegration(wave,band.throughput/wave) 423 modtot = total * (1.0e-23/H) 424 return 1.0/modtot

425

426 -class mJy(FluxUnits):

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

428 - def __init__(self):

429 FluxUnits.__init__(self) 430 self.name = 'mjy' 431 self.nativewave = Hz

432

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

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

435 436

437 - def unitResponse(self,band):

438 #sumfilt(wave,-1,band) 439 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 440 wave=band.wave 441 total = band.trapezoidIntegration(wave,band.throughput/wave) 442 modtot = total * (1.0e-26/H) 443 return 1.0/modtot

444

445 -class muJy(FluxUnits): # New

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

447 - def __init__(self):

448 FluxUnits.__init__(self) 449 self.name = 'mujy' 450 self.nativewave = Hz

451

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

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

454

455 - def unitResponse(self,band):

456 wave=band.wave 457 total = band.trapezoidIntegration(wave,band.throughput/wave) 458 modtot = total * (1.0e-29/H) 459 return 1.0/modtot

460

461 -class nJy(FluxUnits): # New

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

463 - def __init__(self):

464 FluxUnits.__init__(self) 465 self.name = 'njy' 466 self.nativewave = Hz

467

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

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

470

471 - def unitResponse(self,band):

472 wave=band.wave 473 total = band.trapezoidIntegration(wave,band.throughput/wave) 474 modtot = total * (1.0e-32/H) 475 return 1.0/modtot

476

477 -class ABMag(LogFluxUnits):

478 - def __init__(self):

479 LogFluxUnits.__init__(self) 480 self.name = 'abmag' 481 self.linunit = Fnu() 482 self.zeropoint = ABZERO

483 484

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

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

487

488 - def unitResponse(self,band):

489 #sumfilt(wave,-1,band) 490 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 491 wave=band.wave 492 total = band.trapezoidIntegration(wave,band.throughput/wave) 493 modtot = total/H 494 return 2.5*math.log10(modtot) + ABZERO

495

496 -class STMag(LogFluxUnits):

497 - def __init__(self):

498 LogFluxUnits.__init__(self) 499 self.name = 'stmag' 500 self.linunit = Flam() 501 self.zeropoint = STZERO

502 503

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

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

506

507 - def unitResponse(self,band):

508 #sumfilt(wave,1,band) 509 # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ] 510 wave=band.wave 511 total = band.trapezoidIntegration(wave,band.throughput*wave) 512 modtot = total/(H*C) 513 return 2.5*math.log10(modtot) + STZERO

514

515 -class OBMag(LogFluxUnits):

516 - def __init__(self):

517 LogFluxUnits.__init__(self) 518 self.name = 'obmag' 519 self.linunit = Counts() 520 self.zeropoint = 0.0 521 self.isDensity = False

522

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

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

526

527 - def unitResponse(self,band):

528 #sum = asumr(band,nwave) 529 total = band.throughput.sum() 530 return 2.5*math.log10(total)

531

532 -class VegaMag(LogFluxUnits):

533 - def __init__(self):

534 LogFluxUnits.__init__(self) 535 self.name = 'vegamag' 536 self.vegaspec = spectrum.Vega

537

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

539 resampled = self.vegaspec.resample(wave) 540 return resampled.flux * 10.0**(-0.4 * flux)

541

542 - def unitResponse(self,band):

543 sp=band*self.vegaspec 544 total=sp.integrate() 545 return 2.5*math.log10(total)

546

547 -class Counts(FluxUnits):

548 - def __init__(self):

549 FluxUnits.__init__(self) 550 self.name = 'counts' 551 self.isDensity = False

552

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

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

555 556

557 - def unitResponse(self,band):

558 #sum = asumr(band,nwave) 559 total = band.throughput.sum() 560 return 1.0/total

561 562 563 ################ Factory for Units subclasses. ##################### 564 565

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

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

602

603 -def _getDeltaWave(wave):

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

622

Source Code for Module pysynphot.units