1
2 """
3 Definitions of interstellar reddening curves
4
5 Section 1. General interface
6 ============================
7
8 Use the general interface to get different curves:
9
10 >>> wave = np.r_[1e3:1e5:10]
11 >>> for name in ['chiar2006','fitzpatrick1999','fitzpatrick2004','cardelli1989','seaton1979']:
12 ... wave_,mag_ = get_law(name,wave=wave)
13 ... p = pl.plot(1e4/wave_,mag_,label=name)
14 >>> p = pl.xlim(0,10)
15 >>> p = pl.ylim(0,12)
16
17 Use the general interface to get the same curves but with different Rv:
18
19 >>> for Rv in [2.0,3.1,5.1]:
20 ... wave_,mag_ = get_law('cardelli1989',wave=wave,Rv=Rv)
21 ... p = pl.plot(1e4/wave_,mag_,'--',lw=2,label='cardelli1989 Rv=%.1f'%(Rv))
22 >>> p = pl.xlim(0,10)
23 >>> p = pl.ylim(0,12)
24 >>> p = pl.xlabel('1/$\lambda$ [1/$\mu$m]')
25 >>> p = pl.ylabel(r'Extinction E(B-$\lambda$) [mag]')
26 >>> p = pl.legend(prop=dict(size='small'),loc='lower right')
27
28 ]include figure]]ivs_sed_reddening_curves.png]
29
30 Section 2. Individual curve definitions
31 =======================================
32
33 Get the curves seperately:
34
35 >>> wave1,mag1 = cardelli1989()
36 >>> wave2,mag2 = chiar2006()
37 >>> wave3,mag3 = seaton1979()
38 >>> wave4,mag4 = fitzpatrick1999()
39 >>> wave5,mag5 = fitzpatrick2004()
40
41
42 Section 3. Normalisations
43 =========================
44
45 >>> wave = np.logspace(3,6,1000)
46 >>> photbands = ['JOHNSON.V','JOHNSON.K']
47
48 Retrieve two interstellar reddening laws, normalise them to Av and see what
49 the ratio between Ak and Av is. Since the L{chiar2006} law is not defined
50 in the optical, this procedure actually doesn't make sense for that law. In the
51 case of L{cardelli1989}, compute the ratio C{Ak/Av}. Note that you cannot
52 ask for the L{chiar2006} to be normalised in the visual, since the curve is
53 not defined their! If you really want to do that anyway, you need to derive
54 a Ak/Av factor from another curve (e.g. the L{cardelli1989}).
55
56 >>> p = pl.figure()
57 >>> for name,norm in zip(['chiar2006','cardelli1989'],['Ak','Av']):
58 ... wave_,mag_ = get_law(name,wave=wave,norm=norm)
59 ... photwave,photflux = get_law(name,wave=wave,norm=norm,photbands=photbands)
60 ... p = pl.plot(wave_/1e4,mag_,label=name)
61 ... p = pl.plot(photwave/1e4,photflux,'s')
62 ... if name=='cardelli1989': print 'Ak/Av = %.3g'%(photflux[1]/photflux[0])
63 Ak/Av = 0.114
64 >>> p = pl.gca().set_xscale('log')
65 >>> p = pl.gca().set_yscale('log')
66 >>> p = pl.xlabel('Wavelength [micron]')
67 >>> p = pl.ylabel('Extinction A($\lambda$)/Av [mag]')
68
69 ]include figure]]ivs_sed_reddening_02.png]
70
71 Compute the Cardelli law normalised to Ak and Av.
72
73 >>> p = pl.figure()
74 >>> wave_av1,mag_av1 = get_law('cardelli1989',wave=wave,norm='Av')
75 >>> wave_av2,mag_av2 = get_law('cardelli1989',wave=wave,norm='Av',photbands=['JOHNSON.V'])
76 >>> p = pl.plot(wave_av1,mag_av1,'k-',lw=2,label='Av')
77 >>> p = pl.plot(wave_av2,mag_av2,'ks')
78
79 >>> wave_ak1,mag_ak1 = get_law('cardelli1989',wave=wave,norm='Ak')
80 >>> wave_ak2,mag_ak2 = get_law('cardelli1989',wave=wave,norm='Ak',photbands=['JOHNSON.K'])
81 >>> p = pl.plot(wave_ak1,mag_ak1,'r-',lw=2,label='Ak')
82 >>> p = pl.plot(wave_ak2,mag_ak2,'rs')
83
84 >>> p = pl.gca().set_xscale('log')
85 >>> p = pl.gca().set_yscale('log')
86 >>> p = pl.xlabel('Wavelength [micron]')
87 >>> p = pl.ylabel('Extinction A($\lambda$)/A($\mu$) [mag]')
88
89 ]include figure]]ivs_sed_reddening_03.png]
90 """
91
92 import os
93 import numpy as np
94 import logging
95
96 from cc.ivs.io import ascii
97 from cc.ivs.aux import loggers
98 from cc.ivs.aux.decorators import memoized
99 from cc.ivs.units import conversions
100 from cc.ivs.sed import filters
101 import model
102
103 logger = logging.getLogger("SED.RED")
104 logger.addHandler(loggers.NullHandler())
105
106 basename = os.path.join(os.path.dirname(__file__),'redlaws')
107
108
109
110 -def get_law(name,norm='E(B-V)',wave_units='AA',photbands=None,**kwargs):
111 """
112 Retrieve an interstellar reddening law.
113
114 Parameter C{name} must be the function name of one of the laws defined in
115 this module.
116
117 By default, the law will be interpolated on a grid from 100 angstrom to
118 10 micron in steps of 10 angstrom. This can be adjusted with the parameter
119 C{wave} (array), which B{must} be in angstrom. You can change the units
120 ouf the returned wavelength array via C{wave_units}.
121
122 By default, the curve is normalised with respect to E(B-V) (you get
123 A(l)/E(B-V)). You can set the C{norm} keyword to Av if you want A(l)/Av.
124 Remember that
125
126 A(V) = Rv * E(B-V)
127
128 The parameter C{Rv} is by default 3.1, other reasonable values lie between
129 2.0 and 5.1
130
131 Extra accepted keywords depend on the type of reddening law used.
132
133 Example usage:
134
135 >>> wave = np.r_[1e3:1e5:10]
136 >>> wave,mag = get_law('cardelli1989',wave=wave,Rv=3.1)
137
138 @param name: name of the interstellar law
139 @type name: str, one of the functions defined here
140 @param norm: type of normalisation of the curve
141 @type norm: str (one of E(B-V), Av)
142 @param wave_units: wavelength units
143 @type wave_units: str (interpretable for units.conversions.convert)
144 @param photbands: list of photometric passbands
145 @type photbands: list of strings
146 @keyword wave: wavelength array to interpolate the law on
147 @type wave: ndarray
148 @return: wavelength, reddening magnitude
149 @rtype: (ndarray,ndarray)
150 """
151
152 wave_ = kwargs.pop('wave',None)
153 Rv = kwargs.setdefault('Rv',3.1)
154
155
156 wave,mag = globals()[name.lower()](**kwargs)
157 wave_orig,mag_orig = wave.copy(),mag.copy()
158
159
160 if wave_ is not None:
161 if wave_units != 'AA':
162 wave_ = conversions.convert(wave_units,'AA',wave_)
163 mag = np.interp(wave_,wave,mag,right=0)
164 wave = wave_
165
166
167 if norm.lower()=='e(b-v)':
168 mag *= Rv
169 else:
170
171
172 if norm.lower()=='ak':
173 norm = 'JOHNSON.K'
174 elif norm.lower()=='av':
175 norm = 'JOHNSON.V'
176 norm_reddening = model.synthetic_flux(wave_orig,mag_orig,[norm])[0]
177 logger.info('Normalisation via %s: Av/%s = %.6g'%(norm,norm,1./norm_reddening))
178 mag /= norm_reddening
179
180
181 if photbands is not None:
182 mag = model.synthetic_flux(wave,mag,photbands)
183 wave = filters.get_info(photbands)['eff_wave']
184
185
186
187 if wave_units != 'AA':
188 wave = conversions.convert('AA',wave_units,wave)
189
190 return wave,mag
191
192
193 -def redden(flux,wave=None,photbands=None,ebv=0.,rtype='flux',law='cardelli1989',**kwargs):
194 """
195 Redden flux or magnitudes
196
197 The reddening parameters C{ebv} means E(B-V).
198
199 If it is negative, we B{deredden}.
200
201 If you give the keyword C{wave}, it is assumed that you want to (de)redden
202 a B{model}, i.e. a spectral energy distribution.
203
204 If you give the keyword C{photbands}, it is assumed that you want to (de)redden
205 B{photometry}, i.e. integrated fluxes.
206
207 @param flux: fluxes to (de)redden (magnitudes if C{rtype='mag'})
208 @type flux: ndarray (floats)
209 @param wave: wavelengths matching the fluxes (or give C{photbands})
210 @type wave: ndarray (floats)
211 @param photbands: photometry bands matching the fluxes (or give C{wave})
212 @type photbands: ndarray of str
213 @param ebv: reddening parameter E(B-V)
214 @type ebv: float
215 @param rtype: type of dereddening (magnituds or fluxes)
216 @type rtype: str ('flux' or 'mag')
217 @return: (de)reddened flux/magnitude
218 @rtype: ndarray (floats)
219 """
220 if photbands is not None:
221 wave = filters.get_info(photbands)['eff_wave']
222
223 old_settings = np.seterr(all='ignore')
224 wave, reddeningMagnitude = get_law(law,wave=wave,**kwargs)
225
226 if rtype=='flux':
227
228 flux_reddened = flux / 10**(reddeningMagnitude*ebv/2.5)
229 np.seterr(**old_settings)
230 return flux_reddened
231 elif rtype=='mag':
232
233 magnitude = flux
234 magnitude_reddened = magnitude + reddeningMagnitude*ebv
235 np.seterr(**old_settings)
236 return magnitude_reddened
237
238 -def deredden(flux,wave=None,photbands=None,ebv=0.,rtype='flux',**kwargs):
239 """
240 Deredden flux or magnitudes.
241
242 @param flux: fluxes to (de)redden (NOT magnitudes)
243 @type flux: ndarray (floats)
244 @param wave: wavelengths matching the fluxes (or give C{photbands})
245 @type wave: ndarray (floats)
246 @param photbands: photometry bands matching the fluxes (or give C{wave})
247 @type photbands: ndarray of str
248 @param ebv: reddening parameter E(B-V)
249 @type ebv: float
250 @param rtype: type of dereddening (magnituds or fluxes)
251 @type rtype: str ('flux' or 'mag')
252 @return: (de)reddened flux
253 @rtype: ndarray (floats)
254 """
255 return redden(flux,wave=wave,photbands=photbands,ebv=-ebv,rtype=rtype,**kwargs)
256
257
258
259
260
261 @memoized
262 -def chiar2006(Rv=3.1,curve='ism',**kwargs):
263 """
264 Extinction curve at infrared wavelengths from Chiar and Tielens (2006)
265
266 We return A(lambda)/Av, by multiplying A(lambda)/Ak Ak/Av=0.09 (see Chiar
267 and Tielens (2006).
268
269 This is only defined for Rv=3.1. If it is different, this will raise an
270 AssertionError
271
272 Extra kwags are to catch unwanted keyword arguments.
273
274 UNCERTAIN NORMALISATION
275
276 @param Rv: Rv
277 @type Rv: float
278 @param curve: extinction curve
279 @type curve: string (one of 'gc' or 'ism', galactic centre or local ISM)
280 @return: wavelengths (A), A(lambda)/Av
281 @rtype: (ndarray,ndarray)
282 """
283 source = os.path.join(basename,'Chiar2006.red')
284
285
286 assert(Rv==3.1)
287 wavelengths,gc,ism = ascii.read2array(source).T
288 if curve=='gc':
289 alam_ak = gc
290 elif curve=='ism':
291 keep = ism>0
292 alam_ak = ism[keep]
293 wavelengths = wavelengths[keep]
294 else:
295 raise ValueError,'no curve %s'%(curve)
296 alam_aV = alam_ak * 0.09
297
298 return wavelengths*1e4,alam_aV
299
303 """
304 Combined and extrapolated extinction curve Fitzpatrick 2004 and
305 from Chiar and Tielens (2006).
306
307 We return A(lambda)/E(B-V), by multiplying A(lambda)/Av with Rv.
308
309 This is only defined for Rv=3.1. If it is different, this will raise an
310 AssertionError
311
312 Extra kwags are to catch unwanted keyword arguments.
313
314 @param Rv: Rv
315 @type Rv: float
316 @param curve: extinction curve
317 @type curve: string (one of 'gc' or 'ism', galactic centre or local ISM)
318 @return: wavelengths (A), A(lambda)/Av
319 @rtype: (ndarray,ndarray)
320 """
321
322 if curve.lower() not in ['ism','gc']:
323 raise ValueError,'No Fitzpatrick2004/Chiar2006 curve available for %s.'\
324 %(curve)
325
326 fn = 'fitzpatrick2004_chiar2006%s_extrapol.dat'%curve
327 source = os.path.join(basename,fn)
328
329
330 assert(Rv==3.1)
331 wavelengths,alam_ak = ascii.read2array(source).T
332
333
334 wavelengths *= 1e4
335
336
337 norm_reddening = model.synthetic_flux(wavelengths,alam_ak,\
338 ['JOHNSON.V','JOHNSON.K'])
339 ak_to_av = norm_reddening[1]/norm_reddening[0]
340 alam_aV = alam_ak * ak_to_av
341
342 return wavelengths,alam_aV
343
347 """
348 From Fitzpatrick 1999 (downloaded from ASAGIO database)
349
350 This function returns A(lambda)/A(V).
351
352 To get A(lambda)/E(B-V), multiply the return value with Rv (A(V)=Rv*E(B-V))
353
354 Extra kwags are to catch unwanted keyword arguments.
355
356 @param Rv: Rv (2.1, 3.1 or 5.0)
357 @type Rv: float
358 @return: wavelengths (A), A(lambda)/Av
359 @rtype: (ndarray,ndarray)
360 """
361 filename = 'Fitzpatrick1999_Rv_%.1f'%(Rv)
362 filename = filename.replace('.','_') + '.red'
363 myfile = os.path.join(basename,filename)
364 wave,alam_ebv = ascii.read2array(myfile).T
365 alam_av = alam_ebv/Rv
366
367 logger.info('Fitzpatrick1999 curve with Rv=%.2f'%(Rv))
368
369 return wave,alam_av
370
373 """
374 From Fitzpatrick 2004 (downloaded from FTP)
375
376 This function returns A(lambda)/A(V).
377
378 To get A(lambda)/E(B-V), multiply the return value with Rv (A(V)=Rv*E(B-V))
379
380 Extra kwags are to catch unwanted keyword arguments.
381
382 @param Rv: Rv (2.1, 3.1 or 5.0)
383 @type Rv: float
384 @return: wavelengths (A), A(lambda)/Av
385 @rtype: (ndarray,ndarray)
386 """
387 filename = 'Fitzpatrick2004_Rv_%.1f.red'%(Rv)
388 myfile = os.path.join(basename,filename)
389 wave_inv,elamv_ebv = ascii.read2array(myfile,skip_lines=15).T
390
391 logger.info('Fitzpatrick2004 curve with Rv=%.2f'%(Rv))
392
393 return 1e4/wave_inv[::-1],((elamv_ebv+Rv)/Rv)[::-1]
394
398 """
399 Small improvement on Cardelli 1989 by James E. O'Donnell (1994).
400
401 Extra kwags are to catch unwanted keyword arguments.
402
403 @keyword Rv: Rv
404 @type Rv: float
405 @keyword wave: wavelengths to compute the curve on
406 @type wave: ndarray
407 @return: wavelengths (A), A(lambda)/Av
408 @rtype: (ndarray,ndarray)
409 """
410 return cardelli1989(curve='donnell',**kwargs)
411
412
413 @memoized
414 -def cardelli1989(Rv=3.1,curve='cardelli',wave=None,**kwargs):
415 """
416 Construct extinction laws from Cardelli (1989).
417
418 Improvement in optical by James E. O'Donnell (1994)
419
420 wavelengths in Angstrom!
421
422 This function returns A(lambda)/A(V).
423
424 To get A(lambda)/E(B-V), multiply the return value with Rv (A(V)=Rv*E(B-V))
425
426 Extra kwags are to catch unwanted keyword arguments.
427
428 @param Rv: Rv
429 @type Rv: float
430 @param curve: extinction curve
431 @type curve: string (one of 'cardelli' or 'donnell')
432 @param wave: wavelengths to compute the curve on
433 @type wave: ndarray
434 @return: wavelengths (A), A(lambda)/Av
435 @rtype: (ndarray,ndarray)
436 """
437 if wave is None:
438 wave = np.r_[100.:100000.:10]
439
440 all_x = 1./(wave/1.0e4)
441 alam_aV = np.zeros_like(all_x)
442
443
444 infrared = all_x<1.1
445 x = all_x[infrared]
446 ax = +0.574*x**1.61
447 bx = -0.527*x**1.61
448 alam_aV[infrared] = ax + bx/Rv
449
450
451 optical = (1.1<=all_x) & (all_x<3.3)
452 x = all_x[optical]
453 y = x-1.82
454 if curve=='cardelli':
455 ax = 1 + 0.17699*y - 0.50447*y**2 - 0.02427*y**3 + 0.72085*y**4 \
456 + 0.01979*y**5 - 0.77530*y**6 + 0.32999*y**7
457 bx = 1.41338*y + 2.28305*y**2 + 1.07233*y**3 - 5.38434*y**4 \
458 - 0.62251*y**5 + 5.30260*y**6 - 2.09002*y**7
459 elif curve=='donnell':
460 ax = 1 + 0.104*y - 0.609*y**2 + 0.701*y**3 + 1.137*y**4 \
461 - 1.718*y**5 - 0.827*y**6 + 1.647*y**7 - 0.505*y**8
462 bx = 1.952*y + 2.908*y**2 - 3.989*y**3 - 7.985*y**4 \
463 + 11.102*y**5 + 5.491*y**6 -10.805*y**7 + 3.347*y**8
464 else:
465 raise ValueError,'curve %s not found'%(curve)
466 alam_aV[optical] = ax + bx/Rv
467
468
469 ultraviolet = (3.3<=all_x) & (all_x<8.0)
470 x = all_x[ultraviolet]
471 Fax = -0.04473*(x-5.9)**2 - 0.009779*(x-5.9)**3
472 Fbx = +0.21300*(x-5.9)**2 + 0.120700*(x-5.9)**3
473 Fax[x<5.9] = 0
474 Fbx[x<5.9] = 0
475 ax = +1.752 - 0.316*x - 0.104 / ((x-4.67)**2 + 0.341) + Fax
476 bx = -3.090 + 1.825*x + 1.206 / ((x-4.62)**2 + 0.263) + Fbx
477 alam_aV[ultraviolet] = ax + bx/Rv
478
479
480 fuv = 8.0<=all_x
481 x = all_x[fuv]
482 ax = -1.073 - 0.628*(x-8) + 0.137*(x-8)**2 - 0.070*(x-8)**3
483 bx = 13.670 + 4.257*(x-8) - 0.420*(x-8)**2 + 0.374*(x-8)**3
484 alam_aV[fuv] = ax + bx/Rv
485
486 logger.info('%s curve with Rv=%.2f'%(curve.title(),Rv))
487
488 return wave,alam_aV
489
490
491 @memoized
492 -def seaton1979(Rv=3.1,wave=None,**kwargs):
493 """
494 Extinction curve from Seaton, 1979.
495
496 This function returns A(lambda)/A(V).
497
498 To get A(lambda)/E(B-V), multiply the return value with Rv (A(V)=Rv*E(B-V))
499
500 Extra kwags are to catch unwanted keyword arguments.
501
502 @param Rv: Rv
503 @type Rv: float
504 @param wave: wavelengths to compute the curve on
505 @type wave: ndarray
506 @return: wavelengths (A), A(lambda)/Av
507 @rtype: (ndarray,ndarray)
508 """
509 if wave is None: wave = np.r_[1000.:10000.:10]
510
511 all_x = 1e4/(wave)
512 alam_aV = np.zeros_like(all_x)
513
514
515 x_ = np.r_[1.0:2.8:0.1]
516 X_ = np.array([1.36,1.44,1.84,2.04,2.24,2.44,2.66,2.88,3.14,3.36,3.56,3.77,3.96,4.15,4.26,4.40,4.52,4.64])
517 fir = all_x<=2.7
518 alam_aV[fir] = np.interp(all_x[fir][::-1],x_,X_,left=0)[::-1]
519
520
521 infrared = (2.70<=all_x) & (all_x<3.65)
522 x = all_x[infrared]
523 alam_aV[infrared] = 1.56 + 1.048*x + 1.01 / ( (x-4.60)**2 + 0.280)
524
525
526 optical = (3.65<=all_x) & (all_x<7.14)
527 x = all_x[optical]
528 alam_aV[optical] = 2.29 + 0.848*x + 1.01 / ( (x-4.60)**2 + 0.280)
529
530
531 ultraviolet = (7.14<=all_x) & (all_x<=10)
532 x = all_x[ultraviolet]
533 alam_aV[ultraviolet] = 16.17 - 3.20*x + 0.2975*x**2
534
535 logger.info('Seaton curve with Rv=%.2f'%(Rv))
536
537 return wave,alam_aV/Rv
538
539
540
541 if __name__=="__main__":
542 import doctest
543 import pylab as pl
544 doctest.testmod()
545 pl.show()
546