pysynphot.spectrum

1 from __future__ import division 2 3 """ 4 Module: spectrum.py 5 6 Contains SourceSpectrum and SpectralElement class defnitions and 7 their subclasses. 8 9 Also contains the Vega object, which is an instance of a FileSourceSpectrum 10 that can be imported from this file and used for Vega-related calculations. 11 12 Dependencies: 13 ============= 14 pyfits, numpy 15 16 17 """ 18 19 import string 20 import re 21 import os 22 import math 23 24 import pyfits 25 import numpy as N 26 from numpy import ma as MA 27 28 import units 29 import observationmode 30 import locations 31 import planck 32 import pysynphot.exceptions as exceptions #custom pysyn exceptions 33 34 from pysynphot import __version__ 35 try : 36 from pysynphot import __svn_version__ 37 except ImportError : 38 __svn_version__ = 'unk' 39 40 # Renormalization constants from synphot: 41 PI = 3.14159265 # Mysterious math constant 42 RSUN = 6.9599E10 # Radius of sun 43 PC = 3.085678E18 # Parsec 44 RADIAN = RSUN / PC /1000. 45 RENORM = PI * RADIAN * RADIAN # Normalize to 1 solar radius @ 1 kpc 46 47 # MergeWaveSets "too close together" constant 48 MERGETHRESH = 1.e-12 49 50 #Single-precision epsilon value, taken from the synphot FAQ. 51 #This is the minimum separation in wavelength value necessary for 52 #synphot to read the entries as distinct single-precision numbers. 53 syn_epsilon=0.00032 54

55 -def _computeDefaultWaveset():

56 minwave = 500. 57 maxwave = 26000. 58 lenwave = 10000 59 60 w1 = math.log10(minwave) 61 w2 = math.log10(maxwave) 62 63 result = N.zeros(shape=[lenwave,],dtype=N.float64) 64 65 for i in range(lenwave): 66 frac = float(i) / lenwave 67 result[i] = 10 ** (w1 * (1.0 - frac) + w2 * frac) 68 69 return result

70 71 # Default waveset is computed at load time, once and for all. 72 # Note that this is not thread safe. 73 global default_waveset 74 default_waveset = _computeDefaultWaveset() 75 76

77 -def MergeWaveSets(waveset1, waveset2):

78 """Return the union of the two wavesets, unless one or 79 both of them is None.""" 80 if waveset1 is None and waveset2 is not None: 81 MergedWaveSet = waveset2 82 elif waveset2 is None and waveset1 is not None: 83 MergedWaveSet = waveset1 84 elif waveset1 is None and Waveset2 is None: 85 MergedWaveSet = None 86 else: 87 MergedWaveSet = N.union1d(waveset1, waveset2) 88 89 # The merged wave sets may sometimes contain numbers which are nearly 90 # equal but differ at levels as small as 1e-14. Having values this 91 # close together can cause problems down the line so here we test whether 92 # any such small differences are present, with a small difference 93 # defined as less than 1e-12. 94 # 95 # If small differences are present we make a copy of the union'ed array 96 # with the lower of the close together pairs removed. 97 delta = MergedWaveSet[1:] - MergedWaveSet[:-1] 98 99 if not (delta > MERGETHRESH).all(): 100 newlen = len(delta[delta > MERGETHRESH]) + 1 101 newmerged = N.zeros(newlen,dtype=MergedWaveSet.dtype) 102 newmerged[:-1] = MergedWaveSet[delta > MERGETHRESH] 103 newmerged[-1] = MergedWaveSet[-1] 104 105 MergedWaveSet = newmerged 106 107 return MergedWaveSet

108 109

110 -def trimSpectrum(sp, minw, maxw):

111 ''' Creates a new spectrum with trimmed upper and lower ranges. 112 ''' 113 wave = sp.GetWaveSet() 114 flux = sp(wave) 115 116 new_wave = N.compress(wave >= minw, wave) 117 new_flux = N.compress(wave >= minw, flux) 118 119 new_wave = N.compress(new_wave <= maxw, new_wave) 120 new_flux = N.compress(new_wave <= maxw, new_flux) 121 122 result = TabularSourceSpectrum() 123 124 result._wavetable = new_wave 125 result._fluxtable = new_flux 126 127 result.waveunits = units.Units(sp.waveunits.name) 128 result.fluxunits = units.Units(sp.fluxunits.name) 129 130 return result

131 132

133 -class Integrator(object):

134 ''' Integrator engine. 135 '''

136 - def trapezoidIntegration(self,x,y):

137 npoints = x.size 138 if npoints > 0: 139 indices = N.arange(npoints)[:-1] 140 deltas = x[indices+1] - x[indices] 141 integrand = 0.5*(y[indices+1] + y[indices])*deltas 142 sum = integrand.sum() 143 if x[-1]<x[0]: 144 sum*= -1.0 145 return sum 146 else: 147 return 0.0

148

149 - def _columnsFromASCII(self, filename):

150 """ Following synphot/TABLES, ASCII files may contain blank lines, 151 comment lines (beginning with '#'), or terminal comments. This routine 152 may be called by both Spectrum and SpectralElement objects to extract 153 the first two columns from a file.""" 154 155 wlist=[] 156 flist=[] 157 lcount=0 158 fs = open(filename,mode='r') 159 lines=fs.readlines() 160 fs.close() 161 for line in lines: 162 lcount+=1 163 cline=line.strip() 164 if ((len(cline) > 0) and (not cline.startswith('#'))): 165 try: 166 cols=cline.split() 167 if len(cols) >= 2: 168 wlist.append(float(cols[0])) 169 flist.append(float(cols[1])) 170 except Exception, e: 171 raise exceptions.BadRow("Error reading %s: %s"%(filename,str(e)),rows=lcount) 172 173 return wlist, flist

174

175 - def validate_wavetable(self):

176 "Enforce monotonic, ascending wavelengths with no zero values" 177 #First check for invalid values 178 wave=self._wavetable 179 if N.any(wave <= 0): 180 wrong=N.where(wave <= 0)[0] 181 raise exceptions.ZeroWavelength('Negative or Zero wavelength occurs in wavelength array', rows=wrong) 182 183 184 185 #Now check for monotonicity & enforce ascending 186 sorted=N.sort(wave) 187 if not N.alltrue(sorted == wave): 188 if N.alltrue(sorted[::-1] == wave): 189 #monotonic descending is allowed 190 pass 191 else: 192 wrong = N.where(sorted != wave)[0] 193 raise exceptions.UnsortedWavelength('Wavelength array is not monotonic', rows=wrong) 194 195 #Check for duplicate values 196 dw=sorted[1:]-sorted[:-1] 197 if N.any(dw==0): 198 wrong=N.where(dw==0)[0] 199 raise exceptions.DuplicateWavelength("Wavelength array contains duplicate entries",rows=wrong)

200

201 - def validate_fluxtable(self):

202 "Enforce non-negative fluxes" 203 if ((not self.fluxunits.isMag) #neg. magnitudes are legal 204 and (self._fluxtable.min() < 0)): 205 idx=N.where(self._fluxtable < 0) 206 self._fluxtable[idx]=0.0 207 print "Warning, %d of %d bins contained negative fluxes; they have been set to zero."%(len(idx[0]),len(self._fluxtable))

208 209

210 -class SourceSpectrum(Integrator):

211 '''Base class for the Source Spectrum object. 212 ''' 213

214 - def __add__(self, other):

215 '''Source Spectra can be added. Delegate the work to the 216 CompositeSourceSpectrum class. 217 ''' 218 if not isinstance(other, SourceSpectrum): 219 raise TypeError("Can only add two SourceSpectrum objects") 220 221 return CompositeSourceSpectrum(self, other, 'add')

222

223 - def __sub__(self, other):

224 """ Source Spectra can be subtracted, which is just another way 225 of adding.""" 226 227 return self.__add__(-1.0*other)

228

229 - def __mul__(self, other):

230 '''Source Spectra can be multiplied, by constants or by 231 SpectralElement objects. 232 ''' 233 #Multiplying by numeric constants is allowed 234 if isinstance(other, (int, float) ): 235 other = UniformTransmission(other) 236 #so is by SpectralElements. Otherwise, raise an exception. 237 if not isinstance(other, SpectralElement): 238 raise TypeError("SourceSpectrum objects can only be multiplied by SpectralElement objects or constants; %s type detected"%type(other)) 239 240 241 ## Delegate the work of multiplying to CompositeSourceSpectrum 242 243 return CompositeSourceSpectrum(self, other, 'multiply')

244

245 - def __rmul__(self, other):

246 return self.__mul__(other)

247

248 - def addmag(self,magval):

249 """Adding a magnitude is like multiplying a flux. Only works for 250 numbers -- not arrays, spectrum objects, etc""" 251 if N.isscalar(magval): 252 factor = 10**(-0.4*magval) 253 return self*factor 254 else: 255 raise TypeError(".addmag() only takes a constant scalar argument")

256

257 - def getArrays(self):

258 '''Returns wavelength and flux arrays as a tuple, performing 259 units conversion. 260 ''' 261 wave = self.GetWaveSet(); 262 flux = self(wave) 263 264 flux = units.Photlam().Convert(wave, flux, self.fluxunits.name) 265 wave = units.Angstrom().Convert(wave, self.waveunits.name) 266 267 return (wave, flux)

268 269 #Define properties for consistent UI

270 - def _getWaveProp(self):

271 wave,flux=self.getArrays() 272 return wave

273

274 - def _getFluxProp(self):

275 wave,flux=self.getArrays() 276 return flux

277 278 wave=property(_getWaveProp,doc="Wavelength property") 279 flux=property(_getFluxProp,doc="Flux property") 280

281 - def validate_units(self):

282 "Ensure that waveunits are WaveUnits and fluxunits are FluxUnits" 283 if (not isinstance(self.waveunits,units.WaveUnits)): 284 raise TypeError("%s is not a valid WaveUnit"%self.waveunits) 285 if (not isinstance(self.fluxunits,units.FluxUnits)): 286 raise TypeError("%s is not a valid FluxUnit"%self.fluxunits)

287

288 - def writefits(self, filename, clobber=True, trimzero=True, 289 binned=False,precision=None,hkeys=None):

290 """Write the spectrum to a FITS file. 291 filename: name of file to write to 292 clobber=True: Will clobber existing file by default 293 trimzero=True: Will trim zero-flux elements from both ends 294 by default 295 binned=False: Will write in native waveset by default 296 precision: Will write in native precision by default; can be 297 set to "single" or "double" 298 hkeys: Optional dictionary of {keyword:(value,comment)} 299 to be added to primary FITS header 300 """ 301 302 pcodes={'d':'D','s':'E'} 303 if precision is None: 304 precision=self.flux.dtype.char 305 _precision=precision.lower()[0] 306 pcodes={'d':'D','s':'E','f':'E'} 307 308 if clobber: 309 try: 310 os.remove(filename) 311 except OSError: 312 pass 313 314 if binned: 315 wave=self.binwave 316 flux=self.binflux 317 else: 318 wave=self.wave 319 flux=self.flux 320 321 #Add a check for single/double precision clash, so 322 #that if written out in single precision, the wavelength table 323 #will still be sorted with no duplicates 324 #The value of epsilon is taken from the Synphot FAQ. 325 326 if wave.dtype == N.float64 and _precision == 's': 327 idx=N.where(abs(wave[1:]-wave[:-1]) > syn_epsilon) 328 else: 329 idx=N.where(wave) #=> idx=[:] 330 331 wave=wave[idx] 332 flux=flux[idx] 333 334 first,last=0,len(flux) 335 336 337 if trimzero: 338 #Keep one zero at each end 339 nz = flux.nonzero()[0] 340 try: 341 first=max(nz[0]-1,first) 342 last=min(nz[-1]+2,last) 343 except IndexError: 344 pass 345 346 #Construct the columns and HDUlist 347 cw = pyfits.Column(name='WAVELENGTH', 348 array=wave[first:last], 349 unit=self.waveunits.name, 350 format=pcodes[_precision]) 351 cf = pyfits.Column(name='FLUX', 352 array=flux[first:last], 353 unit=self.fluxunits.name, 354 format=pcodes[_precision]) 355 356 #Make the primary header 357 hdu = pyfits.PrimaryHDU() 358 hdulist = pyfits.HDUList([hdu]) 359 360 #User-provided keys are written to the primary header 361 #so are filename and origin 362 bkeys = dict(filename=(os.path.basename(filename),'name of file'), 363 origin=('pysynphot', 364 'Version (%s, %s)'%(__version__,__svn_version__))) 365 #User-values if present may override default values 366 if hkeys is not None: 367 bkeys.update(hkeys) 368 369 #Now update the primary header 370 for key,val in bkeys.items(): 371 hdu.header.update(key, *val) 372 373 374 #Make the extension HDU 375 cols = pyfits.ColDefs([cw, cf]) 376 hdu = pyfits.new_table(cols) 377 378 #There are some standard keywords that should be added 379 #to the extension header. 380 bkeys=dict(expr =(str(self),'pysyn expression'), 381 tdisp1 =('G15.7',), 382 tdisp2 =('G15.7',) 383 ) 384 385 386 try: 387 bkeys['grftable']=(self.bandpass.obsmode.gtname,) 388 bkeys['cmptable']=(self.bandpass.obsmode.ctname,) 389 except AttributeError: 390 pass #Not all spectra have these 391 392 for key,val in bkeys.items(): 393 hdu.header.update(key,*val) 394 395 #Add the header to the list, and write the file 396 hdulist.append(hdu) 397 hdulist.writeto(filename)

398 399

400 - def integrate(self,fluxunits='photlam'):

401 #Extract the flux in the desired units 402 u=self.fluxunits 403 self.convert(fluxunits) 404 wave,flux=self.getArrays() 405 self.convert(u) 406 #then do the integration 407 return self.trapezoidIntegration(wave,flux)

408 409

410 - def sample(self,wave):

411 """Return a flux array, in self.fluxunits, on the provided 412 wavetable""" 413 # convert input wavelengths to Angstroms since the __call__ method 414 # will be expecting that 415 angwave = self.waveunits.ToAngstrom(wave) 416 417 #First use the __call__ to get it in photlam 418 flux = self(angwave) 419 420 #Then convert to the desired units 421 ans = units.Photlam().Convert(angwave,flux,self.fluxunits.name) 422 423 return ans

424

425 - def convert(self, targetunits):

426 '''Convert to other units. This method actually just changes the 427 wavelength and flux units objects, it does not recompute the 428 internally kept wave and flux data; these are kept always in internal 429 units. Method getArrays does the actual computation. 430 ''' 431 nunits = units.Units(targetunits) 432 if nunits.isFlux: 433 self.fluxunits = nunits 434 else: 435 self.waveunits = nunits

436

437 - def redshift(self, z):

438 ''' Returns a new redshifted spectrum. 439 ''' 440 #By default, apply only the doppler shift. 441 442 waveunits=self.waveunits 443 self.convert('angstrom') 444 newwave=self.wave*(1.0+z) 445 copy = ArraySourceSpectrum(wave=newwave, 446 flux=self.flux, 447 waveunits=self.waveunits, 448 fluxunits=self.fluxunits, 449 name="%s at z=%g"%(self.name,z) 450 ) 451 452 self.convert(waveunits) 453 return copy

454

455 - def setMagnitude(self, band, value):

456 '''Makes the magnitude of the source in the band equal to value. 457 band is a SpectralElement. 458 This method is marked for deletion once the .renorm method is 459 well tested. 460 ''' 461 objectFlux = band.calcTotalFlux(self) 462 vegaFlux = band.calcVegaFlux() 463 magDiff = -2.5*math.log10(objectFlux/vegaFlux) 464 factor = 10**(-0.4*(value - magDiff)) 465 466 '''Object returned is a CompositeSourceSpectrum''' 467 468 return self * factor

469

470 - def renorm(self, RNval, RNUnits, band, force=False):

471 """Renormalize the spectrum to the specified value (in the specified 472 flux units) in the specified band. 473 Calls a function in another module to alleviate circular import 474 issues.""" 475 476 from renorm import StdRenorm 477 return StdRenorm(self,band,RNval,RNUnits,force=force)

478

479 - def effstim(self,fluxunits='photlam'):

480 print "?? %s"%fluxunits 481 raise NotImplementedError("Ticket #140: calcphot.effstim functionality")

482

483 -class CompositeSourceSpectrum(SourceSpectrum):

484 '''Composite Source Spectrum object, handles addition, multiplication 485 and keeping track of the wavelength set. 486 '''

487 - def __init__(self, source1, source2, operation):

488 self.component1 = source1 489 self.component2 = source2 490 self.operation = operation 491 self.name=str(self) 492 #Propagate warnings 493 self.warnings={} 494 self.warnings.update(source1.warnings) 495 self.warnings.update(source2.warnings) 496 # for now we keep these attributes here, in spite of the internal 497 # units model. There is code that still breaks down if these attributes 498 # are not here. 499 try: 500 self.waveunits = source1.waveunits 501 self.fluxunits = source1.fluxunits 502 except AttributeError: 503 self.waveunits = source2.waveunits 504 self.fluxunits = source2.fluxunits 505 506 self.isAnalytic = source1.isAnalytic and source2.isAnalytic

507

508 - def __str__(self):

509 opdict = {'add':'+','multiply':'*'} 510 return "%s %s %s"%(str(self.component1),opdict[self.operation],str(self.component2))

511

512 - def __call__(self, wavelength):

513 '''Add or multiply components, delegating the function calculation 514 to the individual objects. 515 ''' 516 if self.operation == 'add': 517 return self.component1(wavelength) + self.component2(wavelength) 518 if self.operation == 'multiply': 519 return self.component1(wavelength) * self.component2(wavelength)

520

521 - def complist(self):

522 ans=[] 523 for comp in (self.component1, self.component2): 524 try: 525 ans.extend(comp.complist()) 526 except AttributeError: 527 ans.append(comp) 528 return ans

529

530 - def GetWaveSet(self):

531 '''Obtain the wavelength set for the composite source by forming 532 the union of wavelengths from each component. 533 ''' 534 waveset1 = self.component1.GetWaveSet() 535 waveset2 = self.component2.GetWaveSet() 536 537 return MergeWaveSets(waveset1, waveset2)

538

539 - def tabulate(self):

540 """Evaluate the spectrum in order to return a tabular source 541 spectrum""" 542 sp=ArraySourceSpectrum(wave=self.wave, 543 flux=self.flux, 544 waveunits=self.waveunits, 545 fluxunits=self.fluxunits, 546 name='%s (tabulated)'%self.name) 547 return sp

548

549 -class TabularSourceSpectrum(SourceSpectrum):

550 '''Class for a source spectrum that is read in from a table. 551 '''

552 - def __init__(self, filename=None, fluxname=None, keepneg=False):

553 self.isAnalytic=False 554 self.warnings={} 555 if filename: 556 self._readSpectrumFile(filename, fluxname) 557 self.filename=filename 558 self.validate_units() 559 self.validate_wavetable() 560 if not keepneg: 561 self.validate_fluxtable() 562 self.ToInternal() 563 self.name=self.filename 564 self.isAnalytic=False 565 else: 566 self._wavetable = None 567 self._fluxtable = None 568 self.waveunits = None 569 self.fluxunits = None 570 self.filename = None 571 self.name=self.filename

572

573 - def _reverse_wave(self):

574 self._wavetable = self._wavetable[::-1]

575 576

577 - def __str__(self):

578 return str(self.name)

579

580 - def _readSpectrumFile(self, filename, fluxname):

581 if filename.endswith('.fits') or filename.endswith('.fit'): 582 self._readFITS(filename, fluxname) 583 else: 584 self._readASCII(filename)

585

586 - def _readFITS(self, filename, fluxname):

587 fs = pyfits.open(filename) 588 589 self._wavetable = fs[1].data.field('wavelength') 590 if fluxname == None: 591 fluxname = 'flux' 592 self._fluxtable = fs[1].data.field(fluxname) 593 594 self.waveunits = units.Units(fs[1].header['tunit1'].lower()) 595 self.fluxunits = units.Units(fs[1].header['tunit2'].lower()) 596 597 fs.close()

598

599 - def _readASCII(self, filename):

600 """ Ascii files have no headers. Following synphot, this 601 routine will assume the first column is wavelength in Angstroms, 602 and the second column is flux in Flam.""" 603 604 self.waveunits = units.Units('angstrom') 605 self.fluxunits = units.Units('flam') 606 wlist,flist = self._columnsFromASCII(filename) 607 self._wavetable=N.array(wlist,dtype=N.float64) 608 self._fluxtable=N.array(flist,dtype=N.float64)

609 610

611 - def __call__(self, wavelengths):

612 '''This is where the flux array is actually calculated given a 613 wavelength array. Returns an array of flux values calculated at 614 the wavelength values input. 615 ''' 616 if N.isscalar(wavelengths): 617 delta=0.0001 618 ww=N.array([wavelengths-delta,wavelengths,wavelengths+delta]) 619 tmp=self.resample(ww) 620 return tmp._fluxtable[1] 621 else: 622 return self.resample(wavelengths)._fluxtable

623 624

625 - def taper(self):

626 '''Taper the spectrum by adding zeros to each end. 627 ''' 628 OutSpec = TabularSourceSpectrum() 629 wcopy = N.zeros(self._wavetable.size+2,dtype=N.float64) 630 fcopy = N.zeros(self._fluxtable.size+2,dtype=N.float64) 631 wcopy[1:-1] = self._wavetable 632 fcopy[1:-1] = self._fluxtable 633 fcopy[0] = 0.0 634 fcopy[-1] = 0.0 635 636 ## The wavelengths to use for the first and last points are 637 ## calculated by using the same ratio as for the 2 interior points 638 wcopy[0] = wcopy[1]*wcopy[1]/wcopy[2] 639 wcopy[-1] = wcopy[-2]*wcopy[-2]/wcopy[-3] 640 641 OutSpec._wavetable = wcopy 642 OutSpec._fluxtable = fcopy 643 OutSpec.waveunits = units.Units(str(self.waveunits)) 644 OutSpec.fluxunits = units.Units(str(self.fluxunits)) 645 646 return OutSpec

647

648 - def resample(self, resampledWaveTab):

649 '''Interpolate flux given a wavelength array that is monotonically 650 increasing and the TabularSourceSpectrum object. 651 @param resampledWaveTab: new wavelength table IN ANGSTROMS 652 @type ressampledWaveTab: ndarray 653 ''' 654 655 ##Check whether the input wavetab is in descending order 656 if resampledWaveTab[0]<resampledWaveTab[-1]: 657 newwave=resampledWaveTab 658 newasc = True 659 else: 660 newwave=resampledWaveTab[::-1] 661 newasc = False 662 663 ## Use numpy interpolation function 664 if self._wavetable[0]<self._wavetable[-1]: 665 oldasc = True 666 ans = N.interp(newwave,self._wavetable, 667 self._fluxtable) 668 else: 669 oldasc = False 670 rev = N.interp(newwave,self._wavetable[::-1], 671 self._fluxtable[::-1]) 672 ans = rev[::-1] 673 674 675 ## If the new and old waveset don't have the same parity, 676 ## the answer has to be flipped again 677 if (newasc != oldasc): 678 ans=ans[::-1] 679 680 ## Finally, make the new object 681 # NB: these manipulations were done using the internal 682 #tables in Angstrom and photlam, so those are the units 683 #that must be fed to the constructor. 684 resampled=ArraySourceSpectrum(wave=resampledWaveTab.copy(), 685 waveunits = 'angstroms', 686 flux = ans.copy(), 687 fluxunits = 'photlam', 688 keepneg=True) 689 690 #Use the convert method to set the units desired by the user. 691 resampled.convert(self.waveunits) 692 resampled.convert(self.fluxunits) 693 694 return resampled

695

696 - def GetWaveSet(self):

697 '''For a TabularSource Spectrum, the WaveSet is just the _wavetable 698 member. Return a copy so that there is no reference to the original 699 object. 700 ''' 701 return self._wavetable.copy()

702

703 - def ToInternal(self):

704 '''Convert to the internal representation of (angstroms, photlam). 705 ''' 706 self.validate_units() 707 savewunits = self.waveunits 708 savefunits = self.fluxunits 709 angwave = self.waveunits.Convert(self.GetWaveSet(), 'angstrom') 710 phoflux = self.fluxunits.Convert(angwave, self._fluxtable, 'photlam') 711 self._wavetable = angwave.copy() 712 self._fluxtable = phoflux.copy() 713 self.waveunits = savewunits 714 self.fluxunits = savefunits

715

716 -class ArraySourceSpectrum(TabularSourceSpectrum):

717 """ spec = ArraySpectrum(numpy array containing wavelenght table, 718 numpy array containing flux table, waveunits, fluxunits, 719 name=human-readable nickname for spectrum, keepneg=True to 720 override the default behavior of setting negative flux values to zero) 721 """

722 - def __init__(self, wave=None, flux=None, 723 waveunits='angstrom', fluxunits='photlam', 724 name='UnnamedArraySpectrum', 725 keepneg=False):

726 """Create a spectrum from arrays. 727 @param wave: Wavelength array 728 @param flux: Flux array 729 @type wave,flux: Numpy array with numerical data 730 731 @param waveunits: Units of wave 732 @param fluxunits: Units of flux 733 @type waveunits: L{units.WaveUnits} or subclass 734 @type fluxunits: L{units.FluxUnits} or subclass 735 736 @param name: Description of this array 737 @type name: string 738 @param keepneg: If true, negative flux values will be retained; by default, they are forced to zero 739 @type keepneg: bool 740 """ 741 if len(wave)!=len(flux): 742 raise ValueError("wave and flux arrays must be of equal length") 743 744 self._wavetable=wave 745 self._fluxtable=flux 746 self.waveunits=units.Units(waveunits) 747 self.fluxunits=units.Units(fluxunits) 748 self.name=name 749 self.isAnalytic=False 750 self.warnings={} 751 752 self.validate_units() #must do before validate_fluxtable because it tests against unit type 753 self.validate_wavetable() #must do before ToInternal in case of descending 754 if not keepneg: 755 self.validate_fluxtable() 756 757 self.ToInternal()

758 759

760 -class FileSourceSpectrum(TabularSourceSpectrum):

761 """spec = FileSpectrum(filename (FITS or ASCII), 762 fluxname=column name containing flux (for FITS tables only), 763 keepneg=True to override thedefault behavior of setting negative 764 flux values to zero)""" 765

766 - def __init__(self, filename, fluxname=None, keepneg=False):

767 """Create a spectrum from a file. 768 @param filename: FITS or ASCII file containing the spectrum 769 @type filename: string 770 771 @param fluxname: Column name specifying the flux (FITS only) 772 @type fluxname: string 773 774 @param keepneg: If true, negative flux values will be retained; by default, they are forced to zero 775 @type keepneg: bool 776 777 """ 778 self.name = locations.irafconvert(filename) 779 self._readSpectrumFile(self.name, fluxname) 780 self.validate_units() 781 self.validate_wavetable() 782 if not keepneg: 783 self.validate_fluxtable() 784 self.ToInternal() 785 self.isAnalytic=False 786 self.warnings={}

787

788 - def _readSpectrumFile(self, filename, fluxname):

789 if filename.endswith('.fits') or filename.endswith('.fit'): 790 self._readFITS(filename, fluxname) 791 else: 792 self._readASCII(filename)

793

794 - def _readFITS(self, filename, fluxname):

795 fs = pyfits.open(filename) 796 797 self._wavetable = fs[1].data.field('wavelength') 798 if fluxname == None: 799 fluxname = 'flux' 800 self._fluxtable = fs[1].data.field(fluxname) 801 self.waveunits = units.Units(fs[1].header['tunit1'].lower()) 802 self.fluxunits = units.Units(fs[1].header['tunit2'].lower()) 803 804 #Retain the header information as a convenience for the user. 805 #If duplicate keywords exist, the value in the extension 806 #header will override that in the primary. 807 self.fheader = dict(fs[0].header) 808 self.fheader.update(dict(fs[1].header)) 809 810 fs.close()

811

812 - def _readASCII(self, filename):

813 """ Ascii files have no headers. Following synphot, this 814 routine will assume the first column is wavelength in Angstroms, 815 and the second column is flux in Flam.""" 816 817 self.waveunits = units.Units('angstrom') 818 self.fluxunits = units.Units('flam') 819 wlist,flist = self._columnsFromASCII(filename) 820 self._wavetable=N.array(wlist,dtype=N.float64) 821 self._fluxtable=N.array(flist,dtype=N.float64) 822 823 #We don't support headers from ascii files 824 self.fheader = dict()

825 826

827 -class AnalyticSpectrum(SourceSpectrum):

828 ''' Base class for analytic functions. These are spectral forms 829 which are defined, by default, on top of the default synphot 830 waveset. 831 ''' 832

833 - def __init__(self,waveunits='angstrom',fluxunits='photlam'):

834 "All AnalyticSpectra must set wave & flux units; do it here." 835 self.waveunits = units.Units(waveunits) 836 self.fluxunits = units.Units(fluxunits) 837 self.validate_units() 838 self.isAnalytic=True 839 self.warnings={}

840

841 - def GetWaveSet(self):

842 global default_waveset 843 return default_waveset.copy()

844 845

846 -class GaussianSource(AnalyticSpectrum):

847 """spec = GaussianSource(TotalFlux under Gaussian, 848 central wavelength of Gaussian, 849 FWHM of Gaussian, 850 waveunits, fluxunits) 851 852 853 """

854 - def __init__(self, flux, center, fwhm, waveunits='angstrom', 855 fluxunits='flam'):

856 AnalyticSpectrum.__init__(self,waveunits,fluxunits) 857 858 self.center = center 859 860 self.fwhm = fwhm 861 862 self.total_flux = flux 863 864 self._input_flux_units = self.fluxunits 865 self._input_wave_units = self.waveunits 866 867 self.sigma = fwhm / math.sqrt(8.0 * math.log(2.0)) 868 869 self.factor = flux / (math.sqrt(2.0 * math.pi) * self.sigma) 870 871 self.name ='Gaussian: mu=%g %s,fwhm=%g %s, total flux=%g %s' % (self.center, 872 self._input_wave_units, 873 self.fwhm, 874 self._input_wave_units, 875 self.total_flux, 876 self._input_flux_units)

877

878 - def __str__(self):

879 return self.name

880

881 - def __call__(self, wavelength):

882 # wavelength comes in as Angstom but Gaussian properties are stored 883 # in user defined units 884 wave = units.Angstrom().Convert(wavelength,self._input_wave_units.name) 885 886 # calculate flux 887 flux = self.factor * \ 888 N.exp(-0.5 * ((wave - self.center) / self.sigma)**2) 889 890 # convert flux to photlam before returning 891 return self._input_flux_units.ToPhotlam(wave,flux)

892

893 - def GetWaveSet(self):

894 '''Return a wavelength set that describes the Gaussian. 895 Overrides the base class to compute 101 values, from 896 center - 5*sigma to center + 5*sigma, in units of 897 0.1*sigma 898 ''' 899 increment = 0.1*self.sigma 900 first = self.center - 50.0*increment 901 last = self.center + 50.0*increment 902 903 return N.arange(first, last, increment)

904 905

906 -class FlatSpectrum(AnalyticSpectrum):

907 """spec = FlatSpectrum(Flux density, waveunits, fluxunits). Defines a 908 flat spectrum in units of fluxunits."""

909 - def __init__(self, fluxdensity, waveunits='angstrom', fluxunits='photlam'):

910 AnalyticSpectrum.__init__(self,waveunits,fluxunits) 911 self.wavelength = None 912 self._fluxdensity = fluxdensity 913 self._input_flux_units = self.fluxunits 914 self.name="Flat spectrum of %g %s"%(self._fluxdensity, 915 self._input_flux_units)

916 - def __str__(self):

917 return self.name

918

919 - def __call__(self, wavelength):

920 if hasattr(wavelength,'shape'): 921 flux = self._fluxdensity*N.ones(wavelength.shape,dtype=N.float64) 922 else: 923 flux = self._fluxdensity 924 925 # __call__ is supposed to return photflam so we need to do the 926 # conversion here since it doesn't make sense to store the _fluxdensity 927 # attribute in photlam 928 wave = units.Angstrom().Convert(wavelength,self.waveunits.name) 929 930 return self._input_flux_units.ToPhotlam(wave,flux)

931 932

933 - def redshift(self, z):

934 """Call the parent's method, which returns a TabularSourceSpectrum, 935 then use its results to create a new FlatSpectrum with the correct 936 value. """ 937 938 tmp=SourceSpectrum.redshift(self,z) 939 ans=FlatSpectrum(tmp.flux.max(), 940 fluxunits=tmp.fluxunits) 941 return ans

942 943 ##This change produces 5 errors and 17 failures in cos_etc_test.py 944 ## def GetWaveSet(self): 945 ## global default_waveset 946 ## return N.array([default_waveset[0],default_waveset[-1]]) 947 948

949 -class Powerlaw(AnalyticSpectrum):

950 """spec=PowerLaw(refwave, exponent, waveunits, fluxunits). 951 952 Power law spectrum of the form (lambda/refval)**exponent, 953 where refval is in Angstroms. 954 The spectrum is normalized to a flux of 1 in "fluxunits" at "refval". 955 """

956 - def __init__(self, refwave, index, waveunits='angstrom', fluxunits='photlam'):

957 AnalyticSpectrum.__init__(self,waveunits,fluxunits) 958 self.wavelength = None 959 960 self._input_flux_units = self.fluxunits 961 self._input_wave_units = self.waveunits 962 963 self._refwave = refwave 964 965 self._index = index 966 967 self.name="Power law: refwave %g %s, index %g"%(self._refwave,self._input_wave_units,self._index)

968

969 - def __str__(self):

970 return self.name

971

972 - def __call__(self, wavelength):

973 # input wavelength is assumed to be angstroms 974 # and either a scalar or a numpy array 975 976 # need to first convert input wavelength to the units the user 977 # specified when creating this object 978 wave = units.Angstrom().Convert(wavelength,self._input_wave_units.name) 979 980 flux = (wave / self._refwave) ** self._index 981 982 # convert flux to photlam before returning 983 return self._input_flux_units.ToPhotlam(wave,flux)

984 985

986 -class BlackBody(AnalyticSpectrum):

987 """ 988 spec = BlackBody(T in Kelvin) 989 990 Blackbody spectrum with specified temperature, in Kelvin. 991 The flux of the spectrum is normalized to a star of solar radius 992 at a distance of 1 kpc.L 993 """

994 - def __init__(self, temperature):

995 waveunits=units.Units('angstrom') 996 fluxunits=units.Units('photlam') 997 AnalyticSpectrum.__init__(self,waveunits,fluxunits) 998 self.wavelength = None 999 self.temperature = temperature 1000 self.name='BB(T=%d)'%self.temperature

1001

1002 - def __str__(self):

1003 return self.name

1004

1005 - def __call__(self, wavelength):

1006 return planck.bbfunc(wavelength, self.temperature) * RENORM

1007 1008

1009 -class SpectralElement(Integrator):

1010 '''Base class for a Spectral Element (e.g. Filter, Detector...). 1011 '''

1012 - def validate_units(self):

1013 "Ensure that waveunits are WaveUnits" 1014 if (not isinstance(self.waveunits,units.WaveUnits)): 1015 raise TypeError("%s is not a valid WaveUnit"%self.waveunits)

1016

1017 - def __mul__(self, other):

1018 '''Permitted to multiply a SpectralElement by another 1019 SpectralElement, or by a SourceSpectrum. In the former 1020 case we return a CompositeSpectralElement, while in the 1021 latter case a CompositeSourceSpectrum. 1022 ''' 1023 if isinstance(other, SpectralElement): 1024 return CompositeSpectralElement(self, other) 1025 1026 if isinstance(other, SourceSpectrum): 1027 return CompositeSourceSpectrum(self, other, 'multiply') 1028 1029 ## Multiplying by a constant is the same as multiplying by a 1030 ## UniformTransmission object 1031 if isinstance(other, (int, float)): 1032 return CompositeSpectralElement(self, UniformTransmission(other)) 1033 1034 else: 1035 print "SpectralElements can only be multiplied by other " + \ 1036 "SpectralElements or SourceSpectrum objects"

1037

1038 - def __rmul__(self, other):

1039 return self.__mul__(other)

1040 1041

1042 - def integrate(self,wave=None):

1043 """Integrate the throughput over the specified waveset, 1044 if None, integrate over the full waveset.""" 1045 if wave is None: 1046 wave=self.wave 1047 ans=self.trapezoidIntegration(wave,self(wave)) 1048 return ans

1049 1050 #.................................................................. 1051 # Methods to implement bandpar functionality go here 1052 #..................................................................

1053 - def avgwave(self):

1054 """Implement the equation for lambda nought as defined 1055 in Koornneef et al 1987, p 836. 1056 Should be equivalent to bandpar.avglam = bandpar.avgwv""" 1057 1058 mywaveunits = self.waveunits.name 1059 self.convert('angstroms') 1060 1061 wave=self.wave 1062 thru=self.throughput 1063 self.convert(mywaveunits) 1064 1065 num = self.trapezoidIntegration(wave, thru*wave) 1066 den = self.trapezoidIntegration(wave, thru) 1067 1068 1069 if 0.0 in (num, den): 1070 return 0.0 1071 else: 1072 return num/den

1073

1074 - def pivot(self,binned=False):

1075 """This is the calculation performed when the ETC invokes calcphot. 1076 Does this need to be calculated on binned waveset, or may 1077 it be calculated on native waveset?""" 1078 if binned: 1079 try: 1080 wave = self.binwave 1081 except AttributeError: 1082 raise AttributeError('Class ' + str(type(self)) + ' does not support binning.') 1083 else: 1084 wave = self.wave 1085 1086 countmulwave = self(wave)*wave 1087 countdivwave = self(wave)/wave 1088 1089 num = self.trapezoidIntegration(wave,countmulwave) 1090 den = self.trapezoidIntegration(wave,countdivwave) 1091 1092 if num == 0.0 or den == 0.0: 1093 return 0.0 1094 1095 return math.sqrt(num/den)

1096

1097 - def rmswidth(self, floor=0):

1098 """Defines the lambda sub rms from Koornneef et al 1987, 1099 p 836; should be definition of bandpar.bandw""" 1100 1101 mywaveunits = self.waveunits.name 1102 self.convert('angstroms') 1103 1104 wave=self.wave 1105 thru=self.throughput 1106 self.convert(mywaveunits) 1107 1108 if floor != 0: 1109 idx = N.where(thru >= floor) 1110 wave = wave[idx] 1111 thru = thru[idx] 1112 1113 integrand = (wave-self.avgwave())**2 * thru 1114 num = self.trapezoidIntegration(wave, integrand) 1115 den = self.trapezoidIntegration(wave, thru) 1116 1117 if 0.0 in (num, den): 1118 return 0.0 1119 else: 1120 ans = math.sqrt(num/den) 1121 return ans

1122

1123 - def rectwidth(self):

1124 """RECTW = INT(THRU) / MAX(THRU)""" 1125 mywaveunits = self.waveunits.name 1126 self.convert('angstroms') 1127 1128 wave=self.wave 1129 thru=self.throughput 1130 self.convert(mywaveunits) 1131 1132 num = self.trapezoidIntegration(wave, thru) 1133 den = thru.max() 1134 1135 if 0.0 in (num, den): 1136 return 0.0 1137 else: 1138 return num/den

1139

1140 - def equivwidth(self):

1141 """ EQUVW = INT(THRU) """ 1142 1143 return self.integrate()

1144

1145 - def efficiency(self):

1146 """QTLAM = dimensionless efficience 1147 = INT(THRU / LAM) 1148 """ 1149 mywaveunits = self.waveunits.name 1150 self.convert('angstroms') 1151 1152 wave=self.wave 1153 thru=self.throughput 1154 self.convert(mywaveunits) 1155 1156 ans = self.trapezoidIntegration(wave, thru/wave) 1157 return ans

1158 1159 #..................................................................

1160 - def check_sig(self, other):

1161 """Only call this if check_overlap returns 'partial'. 1162 Returns True if the LACK of overlap is INsignificant: 1163 i.e., it is ok to go ahead and do whatever we are doing.""" 1164 1165 swave=self.wave[N.where(self.throughput != 0)] 1166 s1,s2=swave.min(),swave.max() 1167 1168 owave=other.wave 1169 o1,o2=owave.min(),owave.max() 1170 1171 lorange=sorted([s1,o1]) 1172 hirange=sorted([s2,o2]) 1173 1174 #Get the full throughput 1175 total=self.integrate() 1176 1177 #Now get the other two pieces 1178 #We cannot yet do 1179 #low=self[slice(*lowrange)].integrate() 1180 wave=self.wave 1181 idxs=[N.searchsorted(wave, lorange, 'left'), 1182 N.searchsorted(wave, hirange, 'left')] 1183 1184 excluded=0.0 1185 for idx in idxs: 1186 try: 1187 excluded+=self.integrate(wave=wave[slice(*idx)]) 1188 except IndexError: 1189 pass #If the range is zero, do nothing 1190 1191 if excluded/total < 0.01: 1192 return True 1193 else: 1194 return False

1195

1196 - def check_overlap(self, other):

1197 """Check whether the wavelength range of other is defined everywhere 1198 that the wavelength range of self is defined. 1199 Returns "full", "partial", "none". 1200 Normally used for checking whether a spectrum is fully defined over 1201 the range of a bandpass. 1202 Note that the full overlap case is asymmetric: if the range of 'self' 1203 extends past the limits of 'other', this will return a partial 1204 overlap. 1205 """ 1206 1207 if other.isAnalytic: 1208 #then it's defined everywhere 1209 return 'full' 1210 1211 swave=self.wave[N.where(self.throughput != 0)] 1212 s1,s2=swave.min(),swave.max() 1213 1214 owave=other.wave 1215 o1,o2=owave.min(),owave.max() 1216 1217 if (s1>=o1 and s2<=o2): 1218 ans='full' 1219 1220 elif (s2<o1) or (o2<s1): 1221 ans='none' 1222 1223 else: 1224 ans='partial' 1225 1226 return ans

1227

1228 - def convert(self, targetunits):

1229 '''Convert to other units. This method actually just changes the 1230 wavelength unit objects, it does not recompute the 1231 internally kept wave data; these are kept always in internal 1232 units. Method getWaveSet does the actual computation.''' 1233 nunits = units.Units(targetunits) 1234 self.waveunits = nunits

1235 1236

1237 - def ToInternal(self):

1238 '''Convert wavelengths to the internal representation of angstroms.. 1239 Note: This is not yet used, but should be for safety when creating 1240 TabularSpectralElements from files. It will also be necessary for the 1241 ArraySpectralElement class that we want to create RSN. 1242 ''' 1243 self.validate_units() 1244 savewunits = self.waveunits 1245 angwave = self.waveunits.Convert(self.GetWaveSet(), 'angstrom') 1246 self._wavetable = angwave.copy() 1247 self.waveunits = savewunits

1248 1249

1250 - def __call__(self, wavelengths):

1251 '''This is where the throughput array is calculated for a given 1252 input wavelength table. 1253 @param wavelengths: an array of wavelengths in Angstroms at which the 1254 throughput should be sampled 1255 @type wavelengths: ndarray 1256 ''' 1257 if N.isscalar(wavelengths): 1258 delta=0.0001 1259 ww=N.array([wavelengths-delta,wavelengths,wavelengths+delta]) 1260 tmp=self.resample(ww) 1261 return tmp._throughputtable[1] 1262 else: 1263 return self.resample(wavelengths)._throughputtable

1264

1265 - def sample(self, wavelengths):

1266 """Provide a more normal user interface to the __call__""" 1267 return self.__call__(wavelengths)

1268 1269

1270 - def taper(self):

1271 '''Taper the spectrum by adding zeros to each end. 1272 ''' 1273 OutElement = TabularSpectralElement() 1274 1275 wcopy = N.zeros(self._wavetable.size+2,dtype=N.float64) 1276 fcopy = N.zeros(self._throughputtable.size+2,dtype=N.float64) 1277 1278 wcopy[1:-1] = self._wavetable 1279 fcopy[1:-1] = self._throughputtable 1280 1281 fcopy[0] = 0.0 1282 fcopy[-1] = 0.0 1283 1284 ## The wavelengths to use for the first and last points are 1285 ## calculated by using the same ratio as for the 2 interior points 1286 wcopy[0] = wcopy[1]*wcopy[1]/wcopy[2] 1287 wcopy[-1] = wcopy[-2]*wcopy[-2]/wcopy[-3] 1288 1289 OutElement._wavetable = wcopy 1290 OutElement._throughputtable = fcopy 1291 1292 return OutElement

1293

1294 - def writefits(self, filename, clobber=True, trimzero=True, 1295 precision=None, hkeys=None):

1296 """Write the bandpass to a FITS file. 1297 filename: name of file to write to 1298 clobber=True: Will clobber existing file by default 1299 trimzero=True: Will trim zero-flux elements from both ends 1300 by default 1301 precision: Will write in native precision by default; can be 1302 set to "single" or "double" 1303 hkeys: Optional dictionary of {keyword:(value,comment)} 1304 to be added to primary FITS header 1305 1306 """ 1307 if precision is None: 1308 precision=self.throughput.dtype.char 1309 _precision=precision.lower()[0] 1310 pcodes={'d':'D','s':'E','f':'E'} 1311 1312 if clobber: 1313 try: 1314 os.remove(filename) 1315 except OSError: 1316 pass 1317 1318 wave=self.wave 1319 thru=self.throughput 1320 1321 #Add a check for single/double precision clash, so 1322 #that if written out in single precision, the wavelength table 1323 #will still be sorted with no duplicates 1324 #The value of epsilon is taken from the Synphot FAQ. 1325 1326 if wave.dtype == N.float64 and _precision == 's': 1327 idx=N.where(abs(wave[1:]-wave[:-1]) > syn_epsilon) 1328 else: 1329 idx=N.where(wave) #=> idx=[:] 1330 1331 wave=wave[idx] 1332 thru=thru[idx] 1333 1334 1335 first,last=0,len(thru) 1336 if trimzero: 1337 #Keep one zero at each end 1338 nz = thru.nonzero()[0] 1339 try: 1340 first=max(nz[0]-1,first) 1341 last=min(nz[-1]+2,last) 1342 except IndexError: 1343 pass 1344 1345 #Construct the columns and HDUlist 1346 cw = pyfits.Column(name='WAVELENGTH', 1347 array=wave[first:last], 1348 unit=self.waveunits.name, 1349 format=pcodes[_precision]) 1350 cf = pyfits.Column(name='THROUGHPUT', 1351 array=thru[first:last], 1352 unit=' ', 1353 format=pcodes[_precision]) 1354 1355 #Make the primary header 1356 hdu = pyfits.PrimaryHDU() 1357 hdulist = pyfits.HDUList([hdu]) 1358 1359 1360 #User-provided keys are written to the primary header; 1361 #so are filename and origin 1362 bkeys = dict(filename=(os.path.basename(filename),'name of file'), 1363 origin=('pysynphot', 1364 'Version (%s, %s)'%(__version__,__svn_version__))) 1365 #User-values if present may override default values 1366 if hkeys is not None: 1367 bkeys.update(hkeys) 1368 1369 #Now update the primary header 1370 for key,val in bkeys.items(): 1371 hdu.header.update(key, *val) 1372 1373 1374 #Make the extension HDU 1375 cols = pyfits.ColDefs([cw, cf]) 1376 hdu = pyfits.new_table(cols) 1377 1378 #There are also some keys to be written to the extension 1379 #header 1380 1381 bkeys=dict(expr =(str(self),'pysyn expression'), 1382 tdisp1=('G15.7',), 1383 tdisp2=('G15.7',) 1384 ) 1385 1386 try: 1387 bkeys['grftable']=(os.path.basename(self.obsmode.gtname), 1388 'graph table used') 1389 bkeys['cmptable']=(os.path.basename(self.obsmode.ctname), 1390 'component table used') 1391 except AttributeError: 1392 pass #Not all bandpasses have these 1393 1394 for key,val in bkeys.items(): 1395 hdu.header.update(key, *val) 1396 1397 #Add the extension to the list, and write to file. 1398 hdulist.append(hdu) 1399 hdulist.writeto(filename)

1400 1401 1402

1403 - def resample(self, resampledWaveTab):

1404 '''Interpolate throughput given a wavelength array that is 1405 monotonically increasing and the TabularSpectralElement object.''' 1406 ##Check whether the input wavetab is in descending order 1407 1408 if resampledWaveTab[0]<resampledWaveTab[-1]: 1409 newwave=resampledWaveTab 1410 newasc = True 1411 else: 1412 newwave=resampledWaveTab[::-1] 1413 newasc = False 1414 1415 ## Use numpy interpolation function 1416 if self._wavetable[0]<self._wavetable[-1]: 1417 oldasc = True 1418 ans = N.interp(newwave,self._wavetable, 1419 self._throughputtable) 1420 else: 1421 oldasc = False 1422 rev = N.interp(newwave,self._wavetable[::-1], 1423 self._throughputtable[::-1]) 1424 ans = rev[::-1] 1425 1426 ## If the new and old waveset don't have the same parity, 1427 ## the answer has to be flipped again 1428 if (newasc != oldasc): 1429 ans=ans[::-1] 1430 1431 # Finally, make the new object. 1432 # NB: these manipulations were done using the internal 1433 #tables in Angstrom, so those are the units 1434 #that must be fed to the constructor. 1435 resampled=ArraySpectralElement(wave=resampledWaveTab.copy(), 1436 waveunits = 'angstroms', 1437 throughput = ans.copy()) 1438 #Use the convert method to set the units desired by the user. 1439 resampled.convert(self.waveunits) 1440 1441 return resampled

1442

1443 - def unitResponse(self):

1444 """Is this correct if waveunits != Angstrom?""" 1445 wave = self.GetWaveSet() 1446 thru = self(wave) 1447 return 1.0 / self.trapezoidIntegration(wave,thru)

1448 1449 1450

1451 - def GetWaveSet(self):

1452 "Return the waveset in the requested units." 1453 wave = units.Angstrom().Convert(self._wavetable, self.waveunits.name) 1454 return wave

1455

1456 - def GetThroughput(self):

1457 """Return the throughput for the internal wavetable.""" 1458 ## NB: Throughput never changes units no matter what the 1459 ## wavelength does. There is an implicit assumption here that 1460 ## the units of the input waveset to the __call__ are always 1461 ## Angstroms. 1462 self.convert('angstroms') 1463 return self.__call__(self.wave)

1464 1465 wave = property(GetWaveSet, doc='Waveset for bandpass') 1466 throughput = property(GetThroughput, doc='Throughput for bandpass') 1467

1468 - def fwhm(self):

1469 raise NotImplementedError("#139: Implement calcband functionality")

1470 1471

1472 -class CompositeSpectralElement(SpectralElement):

1473 '''CompositeSpectralElement Class, which knows how to calculate 1474 its throughput by delegating the calculating the calculating to 1475 its components. 1476 '''

1477 - def __init__(self, component1, component2):

1478 if (not isinstance(component1, SpectralElement) or 1479 not isinstance(component2, SpectralElement)): 1480 raise TypeError("Arguments must be SpectralElements") 1481 self.component1 = component1 1482 self.component2 = component2 1483 self.isAnalytic = component1.isAnalytic and component2.isAnalytic 1484 if component1.waveunits.name == component2.waveunits.name: 1485 self.waveunits = component1.waveunits 1486 else: 1487 msg="Components have different waveunits (%s and %s)"%(component1.waveunits,component2.waveunits) 1488 raise NotImplementedError(msg) 1489 self.throughputunits = None 1490 self.name="(%s * %s)"%(str(self.component1),str(self.component2)) 1491 self.warnings={} 1492 self.warnings.update(component1.warnings) 1493 self.warnings.update(component2.warnings)

1494

1495 - def __call__(self, wavelength):

1496 '''This is where the throughput calculation is delegated. 1497 ''' 1498 return self.component1(wavelength) * self.component2(wavelength)

1499

1500 - def __str__(self):

1501 return self.name

1502

1503 - def complist(self):

1504 ans=[] 1505 for comp in (self.component1, self.component2): 1506 try: 1507 ans.extend(comp.complist()) 1508 except AttributeError: 1509 ans.append(comp) 1510 return ans

1511

1512 - def GetWaveSet(self):

1513 '''This method returns a wavelength set appropriate for a composite 1514 object by forming the union of the wavelengths of the components. 1515 ''' 1516 wave1 = self.component1.GetWaveSet() 1517 wave2 = self.component2.GetWaveSet() 1518 1519 return MergeWaveSets(wave1, wave2)

1520 1521 wave = property(GetWaveSet,doc="wave for CompositeSpectralElement")

1522 1523

1524 -class UniformTransmission(SpectralElement):

1525 '''bandpass=UniformTransmission(dimensionless throughput) 1526 1527 @todo: Need to add a GetWaveSet method (or just return None). 1528 '''

1529 - def __init__(self, value, waveunits='angstrom'):

1530 self.waveunits = units.Units(waveunits) 1531 self.value = value 1532 self.name=str(self) 1533 self.isAnalytic=True 1534 self.warnings={} 1535 #The ._wavetable is used only by the .writefits() method at this time 1536 #It is not for general use. 1537 self._wavetable = N.array([default_waveset[0],default_waveset[-1]])

1538

1539 - def __str__(self):

1540 return "%g"%self.value

1541

1542 - def GetWaveSet(self):

1543 return None

1544 1545 ## This produced 15 test failures in cos_etc_test. 1546 ## def GetWaveSet(self): 1547 ## global default_waveset 1548 ## return N.array([default_waveset[0],default_waveset[-1]]) 1549 ## 1550 ## wave = property(GetWaveSet,doc="wave for UniformTransmission") 1551

1552 - def __call__(self, wavelength):

1553 '''__call__ returns the constant value as an array, given a 1554 wavelength array as argument. 1555 ''' 1556 return 0.0 * wavelength + self.value

1557 1558

1559 -class TabularSpectralElement(SpectralElement):

1560 """bandpass = FileBandpass(FITS or ASCII filename, thrucol= name of 1561 column containing throughput values (for FITS tables only) 1562 """

1563 - def __init__(self, fileName=None, thrucol='throughput'):

1564 '''__init__ takes a character string argument that contains the name 1565 of the file with the spectral element table. 1566 ''' 1567 self.isAnalytic=False 1568 self.warnings={} 1569 if fileName: 1570 if fileName.endswith('.fits') or fileName.endswith('.fit'): 1571 self._readFITS(fileName, thrucol) 1572 else: 1573 self._readASCII(fileName) 1574 self.name = fileName 1575 1576 else: 1577 self.name = None 1578 self._wavetable = None 1579 self._throughputtable = None 1580 self.waveunits = None 1581 self.throughputunits = None

1582

1583 - def _reverse_wave(self):

1584 self._wavetable = self._wavetable[::-1]

1585

1586 - def __str__(self):

1587 return str(self.name)

1588

1589 - def ToInternal(self):

1590 '''Convert wavelengths to the internal representation of angstroms.. 1591 ''' 1592 self.validate_units() 1593 savewunits = self.waveunits 1594 angwave = self.waveunits.Convert(self._wavetable, 'angstrom') 1595 self._wavetable = angwave.copy() 1596 self.waveunits = savewunits

1597 1598

1599 - def _readASCII(self,filename):

1600 """ Ascii files have no headers. Following synphot, this 1601 routine will assume the first column is wavelength in Angstroms, 1602 and the second column is throughput (dimensionless).""" 1603 self.waveunits = units.Units('angstrom') 1604 self.throughputunits = 'none' 1605 wlist,tlist = self._columnsFromASCII(filename) 1606 self._wavetable=N.array(wlist,dtype=N.float64) 1607 self._throughputtable=N.array(tlist,dtype=N.float64)

1608 1609

1610 - def _readFITS(self,filename,thrucol='throughput'):

1611 fs = pyfits.open(filename) 1612 1613 self._wavetable = fs[1].data.field('wavelength') 1614 self._throughputtable = fs[1].data.field(thrucol) 1615 1616 self.waveunits = units.Units(fs[1].header['tunit1'].lower()) 1617 self.throughputunits = 'none' 1618 1619 self.getHeaderKeywords(fs[1].header) 1620 1621 fs.close()

1622 1623

1624 - def getHeaderKeywords(self, header):

1625 ''' This is a placeholder for subclasses to get header keywords without 1626 having to reopen the file again. 1627 ''' 1628 pass

1629 1630

1631 -class ArraySpectralElement(TabularSpectralElement):

1632 """ spec = ArraySpectrum(numpy array containing wavelength table, 1633 numpy array containing throughput table, waveunits, 1634 name=human-readable nickname for bandpass. 1635 """

1636 - def __init__(self, wave=None, throughput=None, 1637 waveunits='angstrom', 1638 name='UnnamedArrayBandpass'):

1639 1640 """Create a spectrum from arrays. 1641 @param wave: Wavelength array 1642 @param throughput: Throughput array 1643 @type wave,throughput: Numpy array with numerical data 1644 1645 @param waveunits: Units of wave 1646 @type waveunits: L{units.WaveUnits} or subclass 1647 1648 @param name: Description of this spectral element 1649 @type name: string 1650 """ 1651 if len(wave)!=len(throughput): 1652 raise ValueError("wave and throughput arrays must be of equal length") 1653 1654 self._wavetable=wave 1655 self._throughputtable=throughput 1656 self.waveunits=units.Units(waveunits) 1657 self.name=name 1658 self.isAnalytic=False 1659 self.warnings={} 1660 1661 self.validate_units() #must do before validate_fluxtable because it tests against unit type 1662 self.validate_wavetable() #must do before ToInternal in case of descending 1663 1664 self.ToInternal()

1665 1666

1667 -class FileSpectralElement(TabularSpectralElement):

1668 """spec = FileSpectrum(filename (FITS or ASCII), 1669 throughputname=column name containing throughput (for FITS tables only), 1670 keepneg=True to override the default behavior of setting negative 1671 throughput values to zero)""" 1672

1673 - def __init__(self, filename, thrucol=None):

1674 """Create a bandpass from a file. 1675 @param filename: FITS or ASCII file containing the bandpass 1676 @type filename: string 1677 1678 @param thrucol: Column name specifying the throughput (FITS only) 1679 @type thrucol: string 1680 1681 """ 1682 self.name = locations.irafconvert(filename) 1683 self._readThroughputFile(self.name, thrucol) 1684 1685 self.validate_units() 1686 self.validate_wavetable() 1687 self.ToInternal() 1688 self.isAnalytic=False 1689 self.warnings={}

1690

1691 - def _readThroughputFile(self, filename, throughputname):

1692 if filename.endswith('.fits') or filename.endswith('.fit'): 1693 self._readFITS(filename, throughputname) 1694 else: 1695 self._readASCII(filename)

1696

1697 - def _readFITS(self, filename, throughputname):

1698 fs = pyfits.open(filename) 1699 1700 self._wavetable = fs[1].data.field('wavelength') 1701 if throughputname == None: 1702 throughputname = 'throughput' 1703 self._throughputtable = fs[1].data.field(throughputname) 1704 self.waveunits = units.Units(fs[1].header['tunit1'].lower()) 1705 1706 #Retain the header information as a convenience for the user. 1707 #If duplicate keywords exist, the value in the extension 1708 #header will override that in the primary. 1709 self.fheader = dict(fs[0].header) 1710 self.fheader.update(dict(fs[1].header)) 1711 1712 fs.close()

1713

1714 - def _readASCII(self, filename):

1715 """ Ascii files have no headers. Following synphot, this 1716 routine will assume the first column is wavelength in Angstroms, 1717 and the second column is throughput in Flam.""" 1718 1719 self.waveunits = units.Units('angstrom') 1720 wlist,flist = self._columnsFromASCII(filename) 1721 self._wavetable=N.array(wlist,dtype=N.float64) 1722 self._throughputtable=N.array(flist,dtype=N.float64) 1723 1724 #We don't support headers from asii files 1725 self.fheader = dict()

1726

1727 -class InterpolatedSpectralElement(SpectralElement):

1728 '''The InterpolatedSpectralElement class handles spectral elements 1729 that are interpolated from columns stored in FITS tables 1730 '''

1731 - def __init__(self, fileName, wavelength):

1732 ''' The file name contains a suffix with a column name specification 1733 in between square brackets, such as [fr388n#]. The wavelength 1734 parameter (poorly named -- it is not always a wavelength) is used to 1735 interpolate between two columns in the file. 1736 ''' 1737 1738 xre=re.search('\[(?P<col>.*?)\]',fileName) 1739 self.name = os.path.expandvars(fileName[0:(xre.start())]) 1740 colSpec = xre.group('col') 1741 1742 self.analytic=False 1743 self.warnings={} 1744 1745 self.interpval = wavelength 1746 1747 fs = pyfits.open(self.name) 1748 1749 # if the file has the PARAMS header keyword and if it is set to 1750 # WAVELENGTH then we want to perform a wavelength shift before 1751 # interpolation, otherwise we don't want to shift. 1752 if 'PARAMS' in fs[0].header and fs[0].header['PARAMS'].lower() == 'wavelength': 1753 doshift = True 1754 else: 1755 doshift = False 1756 1757 #The wavelength table will have to be adjusted before use 1758 wave0 = fs[1].data.field('wavelength') 1759 1760 1761 #Determine the columns that bracket the desired value 1762 # grab all columns that have '#' in them 1763 # then split off the numbers after the '#' 1764 colNames = [n for n in fs[1].data.names if n.find('#') != -1] 1765 colWaves = [float(cn.split('#')[1]) for cn in colNames] 1766 1767 if colNames == []: 1768 raise StandardError('file %s contains no interpolated columns.' % (fileName,)) 1769 1770 waves = MA.array(colWaves) 1771 greater = MA.masked_less(waves, wavelength) 1772 less = MA.masked_greater(waves, wavelength) 1773 upper = MA.minimum(greater) 1774 lower = MA.maximum(less) 1775 1776 if '--' in (str(upper),str(lower)): 1777 raise NotImplementedError("%g outside of range in %s; extrapolation not yet supported"%(wavelength,fileName)) 1778 1779 1780 #Construct the column names 1781 lcol = colNames[MA.argmax(less)] 1782 ucol = colNames[MA.argmin(greater)] 1783 1784 1785 #Extract the data from those columns 1786 lthr = fs[1].data.field(lcol) 1787 uthr = fs[1].data.field(ucol) 1788 1789 if upper != lower: 1790 if doshift: 1791 #Adjust the wavelength table to bracket the range 1792 lwave = wave0 + (lower-self.interpval) 1793 uwave = wave0 + (upper-self.interpval) 1794 1795 #Interpolate the columns at those ranges 1796 lthr = N.interp(lwave, wave0, fs[1].data.field(lcol)) 1797 uthr = N.interp(uwave, wave0, fs[1].data.field(ucol)) 1798 1799 #Then interpolate between the two columns 1800 w = (wavelength - lower) / (upper - lower) 1801 self._throughputtable = uthr * w + lthr * (1.0 - w) 1802 else: 1803 #Interpolate the matching column to the correct wave range 1804 uwave = wave0 + (upper-self.interpval) 1805 uthr = N.interp(uwave, wave0, fs[1].data.field(ucol)) 1806 self._throughputtable = uthr 1807 1808 self._wavetable = wave0 1809 self.waveunits = units.Units(fs[1].header['tunit1'].lower()) 1810 self.throughputunits = 'none' 1811 1812 fs.close()

1813

1814 - def __str__(self):

1815 return "%s#%g"%(self.name,self.interpval)

1816 1817

1818 -class ThermalSpectralElement(TabularSpectralElement):

1819 '''The ThermalSpectralElement class handles spectral elements 1820 that have associated thermal properties read from a FITS table. 1821 1822 ThermalSpectralElements differ from regular SpectralElements in 1823 that they carry thermal parameters such as temperature and beam 1824 filling factor, but otherwise they operate just as regular 1825 SpectralElements. They dont know how to apply themselves to an 1826 existing beam, in the sense that their emissivities should be 1827 handled explictly, outside the objects themselves. 1828 '''

1829 - def __init__(self, fileName):

1830 1831 TabularSpectralElement.__init__(self, fileName=fileName, thrucol='emissivity') 1832 self.warnings={}

1833

1834 - def getHeaderKeywords(self, header):

1835 ''' Overrides base class in order to get thermal keywords. 1836 ''' 1837 self.temperature = header['DEFT'] 1838 self.beamFillFactor = header['BEAMFILL']

1839 1840

1841 -class Box(SpectralElement):

1842 """bandpass = Box(central wavelength, width) - both in Angstroms"""

1843 - def __init__(self, center, width, waveunits=None):

1844 ''' Both center and width are assumed to be in Angstrom 1845 units, according to the synphot definition. 1846 ''' 1847 1848 if waveunits is None: 1849 self.waveunits=units.Units('angstrom') #per docstring: for now 1850 else: 1851 self.waveunits=units.Units(waveunits) 1852 center = self.waveunits.Convert(center, 'angstrom') 1853 width = self.waveunits.Convert(width, 'angstrom') 1854 1855 lower = center - width / 2.0 1856 upper = center + width / 2.0 1857 step = 0.05 # fixed step for now (in A) 1858 1859 self.name='Box at %g (%g wide)'%(center,width) 1860 nwaves = int(((upper - lower) / step)) + 2 1861 self._wavetable = N.zeros(shape=[nwaves,], dtype=N.float64) 1862 for i in range(nwaves): 1863 self._wavetable[i] = lower + step * i 1864 1865 self._wavetable[0] = self._wavetable[1] - step 1866 self._wavetable[-1] = self._wavetable[-2] + step 1867 1868 self._throughputtable = N.ones(shape=self._wavetable.shape, \ 1869 dtype=N.float64) 1870 self._throughputtable[0] = 0.0 1871 self._throughputtable[-1] = 0.0 1872 1873 self.isAnalytic=False 1874 self.warnings={}

1875 1876 Vega = FileSourceSpectrum(locations.VegaFile) 1877

Source Code for Module pysynphot.spectrum