-
Notifications
You must be signed in to change notification settings - Fork 0
/
redlaw.x
374 lines (304 loc) · 13.1 KB
/
redlaw.x
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
include <error.h>
include <mach.h>
#--------------------------------------------------------------------19 Feb 97--
.help redlaw.x Mar95 nebular/fivel
.ih
NAME
gal_redlaw -- Calculates the ISEF for the Galaxy by Savage & Mathis (1979)
ccm_redlaw -- Calculates the ISEF for the Galaxy by Cardelli, Clayton & Mathis (1989)
lmc_redlaw -- Calculates the ISEF for the LMC by Howarth (1983)
smc_redlaw -- Calculates the ISEF for the SMC by Prevot, et al. (1984)
jbk_redlaw -- Returns the ISEF of Whitford (1958), as adapted by Kaler (1976)
.ih
DESCRIPTION
The following procedures return an evaluation of the interstellar
extinction function (ISEF) as published by various authors. Often
the ISEF is expressed as A(lambda)/E(B-V), and is normalized such
that the function is zero at the V passband, and 1.00 at B. These
functions are expressed in the literature as tabulated functions,
polynomial expressions, or both. The parameter is wavelength in
inverse microns, and the result is color excess in magnitudes. The
functions are defined or extended to about lambda=1000 Ang; the form
of the ISEF is not known shortward of that.
The routines here return A(lambda)/A(V) by renormaling the extinction
function to 0.0 at zero inverse wavelength, and 1.0 at V. The ratio
of total to selective absorption in the V passband, R_v, is thought to
average 3.10 for most lines of sight in the Galaxy and in the LMC; it
may be somewhat lower in the SMC. It is the responsibility of the
calling function to re-normalize (i.e., multiply) the returned value
by the choice of R_v. Note that R_v is not a parameter for any of
the ISEF's here except for CCM.
.endhelp
#-------------------------------------------------------------------------------
# Macro for linear interpolation/extrapolaton:
define LIN_INTERP ( (($2)-($1)) / (($4)-($3)) * (($5)-($3)) + ($1))
#-------------------------------------------------------------------------------
# GAL_REDLAW -- Evaluation of the average interstellar extinction function for
# the Galaxy published by Savage & Mathis 1979, ARA&A, vol. 17,
# 73-111. The published extinction function is tabulated as
# A(lambda)/E(B-V); returned value is A(lambda)/A(H-beta).
#
# Extinction:
# lambda < 0.1 um same formula as lambda = 0.1.
# 0.1 < lambda < 3.4 linear interpolation of tabulated funct.
# 3.4 < lambda extrapolate linearly in 1/lam
#
# Revision History:
# 13-Mar-95 by RAShaw Based upon synphot library routine EBMVLFUNC
procedure gal_redlaw (wave, extl, npts)
# Calling arguments:
real wave[ARB] # I: wavelength of interest
real extl[ARB] # O: extinction evaluation array
int npts # I: size of evaluation/wave arrays
# Local variables:
real extab[NTAB] # tabulated extinction at E(B-V)=1.
int i, pix # loop indexes
real val # intermediate value
double x # 1. / (wave, in microns)
double xtab[NTAB] # tabulated inverse wavelengths
define NTAB 28 #
define MIN_WAVE 1000. # Minimum wave where function is defined.
# Tabulated inverse wavelengths in microns:
data (xtab(i),i=1,NTAB) / 0.00, 0.29, 0.45, 0.80, 1.11, 1.43, 1.82,
2.27, 2.50, 2.91, 3.65, 4.00, 4.17, 4.35, 4.57, 4.76, 5.00,
5.26, 5.56, 5.88, 6.25, 6.71, 7.18, 8.00, 8.50, 9.00, 9.50,
10.00 /
# Tabulated extinction function, A(lambda)/E(B-V):
data (extab(i),i=1,NTAB) / 0.00, 0.16, 0.38, 0.87, 1.50, 2.32, 3.10,
4.10, 4.40, 4.90, 6.20, 7.29, 8.00, 8.87, 9.67, 9.33, 8.62,
8.00, 7.75, 7.87, 8.12, 8.15, 8.49, 9.65, 10.55, 11.55, 12.90,
14.40 /
begin
do pix = 1, npts {
if (wave[pix] < MIN_WAVE || IS_INDEFR(wave[pix]) )
call error (1, "GAL_REDLAW: Invalid wavelength")
# Convert wavelength in angstroms to 1/microns
x = 10000.D+0 / wave[pix]
# Linear interpolation of extinction law in 1/lam
for (i=2; i <= NTAB && x > xtab[i]; i=i+1)
;
val = LIN_INTERP (extab[i-1], extab[i], xtab[i-1], xtab[i], x)
# Renormalize extinction function to A(lambda)/A(V)
extl[pix] = val / 3.1
}
end
#-------------------------------------------------------------------------------
# CCM_REDLAW - Computes the interstellar extinction function A(lambda)/A(V)
# of Cardelli, Clayton & Mathis 1989, ApJ 345, 245-256.
#
# 18 May 93 by RAShaw: Initial implementation, based upon CCM module
# in onedspec.deredden
procedure ccm_redlaw (wave, extl, npts, rv)
# Calling argument:
real wave[ARB] # I: wavelength of emission line, Angstroms
real extl[ARB] # O: extinction evaluation array
int npts # I: size of evaluation/wave arrays
real rv # I: ratio of total to selective extinction
# Declarations:
real a, b # coefficients of extinction function
int pix # loop index
double x # wavelength in inverse microns
double y # generic value
define MIN_WAVE 1000. # Minimum wave where function is defined.
begin
do pix = 1, npts {
if (wave[pix] < MIN_WAVE || IS_INDEFR(wave[pix]) )
call error (1, "CCM_REDLAW: Invalid wavelength")
# Convert input wavelength to inverse microns
x = 10000.D+0 / wave[pix]
# For wavelengths longward of the L passband, linearly interpolate
# to 0. at 1/lambda = 0. (a=0.08, b=-0.0734 @ x=0.29)
if (x < 0.29D+0) {
a = 0.2759 * x
b = -0.2531 * x
# Compute a(x) and b(x)
} else if (x < 1.1D+0) {
y = x ** 1.61
a = 0.574 * y
b = -0.527 * y
} else if (x < 3.3D+0) {
y = x - 1.82
a = 1 + y * (0.17699 + y * (-0.50447 + y * (-0.02427 +
y * (0.72085 + y * (0.01979 + y * (-0.77530 +
y * 0.32999))))))
b = y * (1.41338 + y * (2.28305 + y * (1.07233 + y *
(-5.38434 + y * (-0.62251 + y * (5.30260 - y * 2.09002))))))
} else if (x < 5.9D+0) {
y = (x - 4.67) ** 2
a = 1.752 - 0.316 * x - 0.104 / (y + 0.341)
b = -3.090 + 1.825 * x + 1.206 / (y + 0.263)
} else if (x < 8.0D+0) {
y = (x - 4.67) ** 2
a = 1.752 - 0.316 * x - 0.104 / (y + 0.341)
b = -3.090 + 1.825 * x + 1.206 / (y + 0.263)
y = x - 5.9
a = a - 0.04473 * y**2 - 0.009779 * y**3
b = b + 0.2130 * y**2 + 0.1207 * y**3
# Truncate range of ISEF to that for 1000 Ang.
} else if (x <= 10.0D+0) {
x = min (x, 10.0D+0)
y = x - 8.D+0
a = -1.072 - 0.628 * y + 0.137 * y**2 - 0.070 * y**3
b = 13.670 + 4.257 * y - 0.420 * y**2 + 0.374 * y**3
}
# Compute A(lambda)/A(V)
y = a + b / rv
extl[pix] = y
}
end
#-------------------------------------------------------------------------------
# LMC_REDLAW - Evaluation of LMC extinction curve as published by Howarth
# 1983, MNRAS, 203, 301. Functional fits from paper are
# A(lambda)/E(B-V); returned value is A(lambda)/A(V).
#
# Revision History:
# 18-Oct-94 by RAShaw Initial implementation
# 14-Mar-95 by RAShaw Return A(lambda)/A(V) instead
procedure lmc_redlaw (wave, extl, npts)
# Calling arguments:
real wave[ARB] #I: evaluation wavelength array
real extl[ARB] #O: extinction evaluation array
int npts #I: size of evaluation/wave arrays
# Local variables:
double delt # work variable
real extab[NTAB] # tabulated extinction at E(B-V)=1.
int i, pix # loop index
real val # intermediate value
real x # 1. / (wave, in microns)
real xtab[NTAB] # tabulated inverse wavelengths
define NTAB 7
define MIN_WAVE 1000. # Minimum wave where function is defined.
# Tabulated inverse wavelengths in microns:
data (xtab(i),i=1,NTAB) / 0.00, 0.29, 0.45, 0.80, 1.11, 1.43, 1.82 /
# Tabulated extinction function, A(lambda)/E(B-V), from Savage & Mathis:
data (extab(i),i=1,NTAB) / 0.00, 0.16, 0.38, 0.87, 1.50, 2.32, 3.10 /
begin
do pix = 1, npts {
if (wave[pix] < MIN_WAVE || IS_INDEFR(wave[pix]) )
call error (1, "LMC_REDLAW: Invalid wavelength")
# Convert wavelength in angstroms to 1/microns
x = 10000. / wave[pix]
# Infrared - optical: linear interpolation of Savage & Mathis 1979
if ( x <= 1.82) {
for (i=2; i <= NTAB && x > xtab[i]; i=i+1)
;
val = LIN_INTERP (extab[i-1], extab[i], xtab[i-1], xtab[i], x)
# The following polynomial evaluations assume R = 3.1
# Violet
} else if ( x <= 2.75 ) {
delt = x - 1.82
val = 3.1 + (2.04 + 0.094 * delt) * delt
# Ultraviolet out to lambda = 1000 A
} else {
x = min (x, 10.0)
delt = x - 4.557
val = 3.1 - 0.236 + 0.462 * x + 0.105 * x * x +
0.454 / (delt**2 + 0.293)
}
# Renormalize extinction function to A(lambda)/A(V)
extl[pix] = val / 3.1
}
end
#-------------------------------------------------------------------------------
# SMC_REDLAW - Computes the interstellar extinction function of Prevot et al.
# (1984), A&A, 132, 389-392. Tabulated values from paper are
# E(lambda)/E(B-V); returned value is A(lambda)/A(v).
#
# 20-Sep-94 by RAShaw: Initial implementation
# 14-Mar-95 by RAShaw Return A(lambda)/A(V) instead
procedure smc_redlaw (wave, extl, npts)
# Calling argument:
real wave[ARB] # I: wavelength of emission line, Angstroms
real extl[ARB] # O: extinction evaluation array
int npts # I: size of evaluation/wave arrays
# Local variables:
real extab[NTAB] # tabulated extinction function from Prevot
int i, pix # loop indexes
real val # generic
real x # wavelength in inverse microns
real xtab[NTAB] # tabulated inverse wavelengths from Prevot
define NTAB 33
define MIN_WAVE 1000. # Minimum wave where function is defined.
# Tabulated inverse wavelengths in microns:
data (xtab(i),i=1,NTAB) / 0.00, 0.29, 0.45, 0.80, 1.11, 1.43, 1.82,
2.35, 2.70, 3.22, 3.34, 3.46, 3.60, 3.75, 3.92, 4.09, 4.28,
4.50, 4.73, 5.00, 5.24, 5.38, 5.52, 5.70, 5.88, 6.07, 6.27,
6.48, 6.72, 6.98, 7.23, 7.52, 7.84 /
# Tabulated extinction function, E(lambda-V)/E(B-V):
data (extab(i),i=1,NTAB) /-3.10, -2.94, -2.72, -2.23, -1.60, -0.78, 0.00,
1.00, 1.67, 2.29, 2.65, 3.00, 3.15, 3.49, 3.91, 4.24, 4.53,
5.30, 5.85, 6.38, 6.76, 6.90, 7.17, 7.71, 8.01, 8.49, 9.06,
9.28, 9.84, 10.80, 11.51, 12.52, 13.54 /
begin
do pix = 1, npts {
if (wave[pix] < MIN_WAVE || IS_INDEFR(wave[pix]) )
call error (1, "SMC_REDLAW: Invalid wavelength")
# Convert wavelength in angstroms to 1/microns
x = 10000. / wave[pix]
x = min (x, 10.0)
# Linearly interpolate extinction law in 1/lam
for (i=2; i <= NTAB && x > xtab[i]; i=i+1)
;
val = LIN_INTERP (extab[i-1], extab[i], xtab[i-1], xtab[i], x)
# Renormalize extinction function to A(lambda)/A(V)
extl[pix] = val / 3.1 + 1.
}
end
#-------------------------------------------------------------------------------
# JBK_REDLAW - Computes the interstellar extinction function of Whitford
# (1958), extended to the UV by Seaton (1977), as adapted by
# Kaler (1976).
#
# 13 May 93 by RAShaw: initial implementation
procedure jbk_redlaw (wave, extl, npts)
# Calling argument:
real wave[ARB] #I: wavelength of emission line, Angstroms
real extl[ARB] #O: extinction evaluation array
int npts #I: size of evaluation/wave arrays
# Local variables:
int i, pix # loop indexes
real x # wavelength in inverse microns
real x1, x2 # bounding reference wavelengths, in micron^-1
real extab[NTAB] # tabulated extinction function from Kaler (1976)
real refw[NTAB] # reference wavelengths
define NTAB 95 #
define MIN_WAVE 1000. # Minimum wave where function is defined.
# Tabulated wavelengths, Angstroms:
data refw /1150., 1200., 1250., 1300., 1350., 1400., 1450., 1500., 1550.,
1600., 1650., 1700., 1750., 1800., 1850., 1900., 1950., 2000., 2050.,
2100., 2150., 2200., 2250., 2300., 2350., 2400., 2450., 2500., 2550.,
2600., 2650., 2700., 2750., 2800., 2850., 2900., 2950., 3000., 3050.,
3100., 3333., 3500., 3600., 3700., 3800., 3900., 4000., 4100., 4200.,
4300., 4400., 4500., 4600., 4700., 4800., 4861.3, 5000., 5100., 5200.,
5300., 5400., 5500., 5600., 5700., 5800., 5900., 6000., 6100., 6200.,
6300., 6400., 6500., 6600., 6700., 6800., 6900., 7000., 7200., 7400.,
7600., 7800., 8000., 8200., 8400., 8600., 8800., 9000., 9500., 10000.,
11000.,12000.,14000.,16000.,20000.,1.D+6/
# Tabulated extinction function:
data extab /1.96, 1.78, 1.61, 1.49, 1.37, 1.29, 1.24, 1.20, 1.20,
1.20, 1.17, 1.13, 1.11, 1.10, 1.12, 1.17, 1.25, 1.35, 1.45,
1.53, 1.60, 1.62, 1.52, 1.40, 1.28, 1.17, 1.06, 0.98, 0.9,
0.84, 0.77, 0.72, 0.68, 0.64, 0.60, 0.57, 0.53, 0.51, 0.48,
0.46, 0.385, 0.358, 0.33, 0.306, 0.278, 0.248, 0.220, 0.195, 0.168,
0.143, 0.118, 0.095, 0.065, 0.040, 0.015, 0.000,-0.030,-0.055,-0.078,
-0.10, -0.121,-0.142,-0.164,-0.182,-0.201,-0.220,-0.238,-0.254,-0.273,
-0.291,-0.306,-0.321,-0.337,-0.351,-0.365,-0.377,-0.391,-0.416,-0.441,
-0.465,-0.490,-0.510,-0.529,-0.548,-0.566,-0.582,-0.597,-0.633,-0.663,
-0.718,-0.763,-0.840,-0.890,-0.960,-1.000/
begin
do pix = 1, npts {
if (wave[pix] < MIN_WAVE || IS_INDEFR(wave[pix]) )
call error (1, "SMC_REDLAW: Invalid wavelength")
for (i=2; i <= NTAB && refw[i] < wave[pix] ; i=i+1)
;
if ( abs (wave[pix] - refw[i]) < EPSILONR )
extl[pix] = extab[i]
else {
x = 10000. / wave[pix]
# x = min (x, 10.0)
x1 = 10000. / refw[i-1]
x2 = 10000. / refw[i]
extl[pix] = LIN_INTERP (extab[i-1], extab[i], x1, x2, x)
}
}
end