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 ivs.inout import ascii
97 from ivs.aux import loggers
98 from ivs.aux.decorators import memoized
99 from ivs.units import conversions
100 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' and photbands is not None:
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)/E(B-V), by multiplying A(lambda)/Av with Rv.
267
268 This is only defined for Rv=3.1. If it is different, this will raise an
269 AssertionError
270
271 Extra kwags are to catch unwanted keyword arguments.
272
273 UNCERTAIN NORMALISATION
274
275 @param Rv: Rv
276 @type Rv: float
277 @param curve: extinction curve
278 @type curve: string (one of 'gc' or 'ism', galactic centre or local ISM)
279 @return: wavelengths (A), A(lambda)/Av
280 @rtype: (ndarray,ndarray)
281 """
282 source = os.path.join(basename,'Chiar2006.red')
283
284
285 assert(Rv==3.1)
286 wavelengths,gc,ism = ascii.read2array(source).T
287 if curve=='gc':
288 alam_ak = gc
289 elif curve=='ism':
290 keep = ism>0
291 alam_ak = ism[keep]
292 wavelengths = wavelengths[keep]
293 else:
294 raise ValueError,'no curve %s'%(curve)
295 alam_aV = alam_ak * 0.09
296
297 return wavelengths*1e4,alam_aV
298
303 """
304 From Fitzpatrick 1999 (downloaded from ASAGIO database)
305
306 This function returns A(lambda)/A(V).
307
308 To get A(lambda)/E(B-V), multiply the return value with Rv (A(V)=Rv*E(B-V))
309
310 Extra kwags are to catch unwanted keyword arguments.
311
312 @param Rv: Rv (2.1, 3.1 or 5.0)
313 @type Rv: float
314 @return: wavelengths (A), A(lambda)/Av
315 @rtype: (ndarray,ndarray)
316 """
317 filename = 'Fitzpatrick1999_Rv_%.1f'%(Rv)
318 filename = filename.replace('.','_') + '.red'
319 myfile = os.path.join(basename,filename)
320 wave,alam_ebv = ascii.read2array(myfile).T
321 alam_av = alam_ebv/Rv
322
323 logger.info('Fitzpatrick1999 curve with Rv=%.2f'%(Rv))
324
325 return wave,alam_av
326
329 """
330 From Fitzpatrick 2004 (downloaded from FTP)
331
332 This function returns A(lambda)/A(V).
333
334 To get A(lambda)/E(B-V), multiply the return value with Rv (A(V)=Rv*E(B-V))
335
336 Extra kwags are to catch unwanted keyword arguments.
337
338 @param Rv: Rv (2.1, 3.1 or 5.0)
339 @type Rv: float
340 @return: wavelengths (A), A(lambda)/Av
341 @rtype: (ndarray,ndarray)
342 """
343 filename = 'Fitzpatrick2004_Rv_%.1f.red'%(Rv)
344 myfile = os.path.join(basename,filename)
345 wave_inv,elamv_ebv = ascii.read2array(myfile,skip_lines=15).T
346
347 logger.info('Fitzpatrick2004 curve with Rv=%.2f'%(Rv))
348
349 return 1e4/wave_inv[::-1],((elamv_ebv+Rv)/Rv)[::-1]
350
354 """
355 Small improvement on Cardelli 1989 by James E. O'Donnell (1994).
356
357 Extra kwags are to catch unwanted keyword arguments.
358
359 @keyword Rv: Rv
360 @type Rv: float
361 @keyword wave: wavelengths to compute the curve on
362 @type wave: ndarray
363 @return: wavelengths (A), A(lambda)/Av
364 @rtype: (ndarray,ndarray)
365 """
366 return cardelli1989(curve='donnell',**kwargs)
367
368
369 @memoized
370 -def cardelli1989(Rv=3.1,curve='cardelli',wave=None,**kwargs):
371 """
372 Construct extinction laws from Cardelli (1989).
373
374 Improvement in optical by James E. O'Donnell (1994)
375
376 wavelengths in Angstrom!
377
378 This function returns A(lambda)/A(V).
379
380 To get A(lambda)/E(B-V), multiply the return value with Rv (A(V)=Rv*E(B-V))
381
382 Extra kwags are to catch unwanted keyword arguments.
383
384 @param Rv: Rv
385 @type Rv: float
386 @param curve: extinction curve
387 @type curve: string (one of 'cardelli' or 'donnell')
388 @param wave: wavelengths to compute the curve on
389 @type wave: ndarray
390 @return: wavelengths (A), A(lambda)/Av
391 @rtype: (ndarray,ndarray)
392 """
393 if wave is None:
394 wave = np.r_[100.:100000.:10]
395
396 all_x = 1./(wave/1.0e4)
397 alam_aV = np.zeros_like(all_x)
398
399
400 infrared = all_x<1.1
401 x = all_x[infrared]
402 ax = +0.574*x**1.61
403 bx = -0.527*x**1.61
404 alam_aV[infrared] = ax + bx/Rv
405
406
407 optical = (1.1<=all_x) & (all_x<3.3)
408 x = all_x[optical]
409 y = x-1.82
410 if curve=='cardelli':
411 ax = 1 + 0.17699*y - 0.50447*y**2 - 0.02427*y**3 + 0.72085*y**4 \
412 + 0.01979*y**5 - 0.77530*y**6 + 0.32999*y**7
413 bx = 1.41338*y + 2.28305*y**2 + 1.07233*y**3 - 5.38434*y**4 \
414 - 0.62251*y**5 + 5.30260*y**6 - 2.09002*y**7
415 elif curve=='donnell':
416 ax = 1 + 0.104*y - 0.609*y**2 + 0.701*y**3 + 1.137*y**4 \
417 - 1.718*y**5 - 0.827*y**6 + 1.647*y**7 - 0.505*y**8
418 bx = 1.952*y + 2.908*y**2 - 3.989*y**3 - 7.985*y**4 \
419 + 11.102*y**5 + 5.491*y**6 -10.805*y**7 + 3.347*y**8
420 else:
421 raise ValueError,'curve %s not found'%(curve)
422 alam_aV[optical] = ax + bx/Rv
423
424
425 ultraviolet = (3.3<=all_x) & (all_x<8.0)
426 x = all_x[ultraviolet]
427 Fax = -0.04473*(x-5.9)**2 - 0.009779*(x-5.9)**3
428 Fbx = +0.21300*(x-5.9)**2 + 0.120700*(x-5.9)**3
429 Fax[x<5.9] = 0
430 Fbx[x<5.9] = 0
431 ax = +1.752 - 0.316*x - 0.104 / ((x-4.67)**2 + 0.341) + Fax
432 bx = -3.090 + 1.825*x + 1.206 / ((x-4.62)**2 + 0.263) + Fbx
433 alam_aV[ultraviolet] = ax + bx/Rv
434
435
436 fuv = 8.0<=all_x
437 x = all_x[fuv]
438 ax = -1.073 - 0.628*(x-8) + 0.137*(x-8)**2 - 0.070*(x-8)**3
439 bx = 13.670 + 4.257*(x-8) - 0.420*(x-8)**2 + 0.374*(x-8)**3
440 alam_aV[fuv] = ax + bx/Rv
441
442 logger.info('%s curve with Rv=%.2f'%(curve.title(),Rv))
443
444 return wave,alam_aV
445
446
447 @memoized
448 -def seaton1979(Rv=3.1,wave=None,**kwargs):
449 """
450 Extinction curve from Seaton, 1979.
451
452 This function returns A(lambda)/A(V).
453
454 To get A(lambda)/E(B-V), multiply the return value with Rv (A(V)=Rv*E(B-V))
455
456 Extra kwags are to catch unwanted keyword arguments.
457
458 @param Rv: Rv
459 @type Rv: float
460 @param wave: wavelengths to compute the curve on
461 @type wave: ndarray
462 @return: wavelengths (A), A(lambda)/Av
463 @rtype: (ndarray,ndarray)
464 """
465 if wave is None: wave = np.r_[1000.:10000.:10]
466
467 all_x = 1e4/(wave)
468 alam_aV = np.zeros_like(all_x)
469
470
471 x_ = np.r_[1.0:2.8:0.1]
472 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])
473 fir = all_x<=2.7
474 alam_aV[fir] = np.interp(all_x[fir][::-1],x_,X_,left=0)[::-1]
475
476
477 infrared = (2.70<=all_x) & (all_x<3.65)
478 x = all_x[infrared]
479 alam_aV[infrared] = 1.56 + 1.048*x + 1.01 / ( (x-4.60)**2 + 0.280)
480
481
482 optical = (3.65<=all_x) & (all_x<7.14)
483 x = all_x[optical]
484 alam_aV[optical] = 2.29 + 0.848*x + 1.01 / ( (x-4.60)**2 + 0.280)
485
486
487 ultraviolet = (7.14<=all_x) & (all_x<=10)
488 x = all_x[ultraviolet]
489 alam_aV[ultraviolet] = 16.17 - 3.20*x + 0.2975*x**2
490
491 logger.info('Seaton curve with Rv=%.2f'%(Rv))
492
493 return wave,alam_aV/Rv
494
495
496
497 if __name__=="__main__":
498 import doctest
499 import pylab as pl
500 doctest.testmod()
501 pl.show()
502