-
Notifications
You must be signed in to change notification settings - Fork 0
/
redlaw-old.x
346 lines (285 loc) · 10.8 KB
/
redlaw-old.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
include <mach.h>
#--------------------------------------------------------------------17 Jun 96--
.help redlaw.x Jan96 nebular/lib
.ih
NAME
gal_redlaw -- Calculates the ISEF for the Galaxy by Savage & Mathis (1979)
lmc_redlaw -- Calculates the ISEF for the LMC by Howarth (1983)
smc_redlaw -- Calculates the ISEF for the SMC by Prevot, et al. (1984)
ccm_redlaw -- Calculates the ISEF for the Galaxy by Cardelli, Clayton & Mathis (1989)
.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.
.ih
ROUTINE DETAILS
.ls procedure gal_redlaw (wave, extl, npts)
.ls ARGUMENTS
.ls wave[ARB] (input: double)
wavelength of interest
.le
.ls extl[ARB] (output: double)
extinction evaluation array
.le
.ls npts (input: int)
size of evaluation/wave arrays
.le
.le
.le
.ls procedure ccm_redlaw (wave, extl, npts, rv)
.ls ARGUMENTS
.ls wave[ARB] (input: double)
wavelength of interest
.le
.ls extl[ARB] (output: double)
extinction evaluation array
.le
.ls npts (input: int)
size of evaluation/wave arrays
.le
.ls rv (input: double)
ratio of total to selective extinction
.le
.le
.le
.ls procedure lmc_redlaw (wave, extl, npts)
.ls ARGUMENTS
.ls wave[ARB] (input: double)
wavelength of interest
.le
.ls extl[ARB] (output: double)
extinction evaluation array
.le
.ls npts (input: int)
size of evaluation/wave arrays
.le
.le
.le
.ls procedure smc_redlaw (wave, extl, npts)
.ls ARGUMENTS
.ls wave[ARB] (input: double)
wavelength of interest
.le
.ls extl[ARB] (output: double)
extinction evaluation array
.le
.ls npts (input: int)
size of evaluation/wave arrays
.le
.le
.le
.endhelp
#-------------------------------------------------------------------------------
# Macro for linear interpolation/extrapolaton:
define LIN_INT ( (($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
procedure gal_redlaw (wave, extl, npts)
# Calling arguments:
double wave[ARB] #I: wavelength of interest
double extl[ARB] #O: extinction evaluation array
int npts #I: size of evaluation/wave arrays
# Local variables:
double 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
# 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] < EPSILOND || IS_INDEFD(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_INT (extab[i-1], extab[i], xtab[i-1], xtab[i], x)
# Renormalize extinction function to A(lambda)/A(V)
extl[pix] = val / 3.1D+0
}
end
#-------------------------------------------------------------------------------
# CCM_REDLAW - Computes the interstellar extinction function A(lambda)/A(V)
# of Cardelli, Clayton & Mathis 1989, ApJ 345, 245-256.
procedure ccm_redlaw (wave, extl, npts, rv)
# Calling argument:
double wave[ARB] # I: wavelength of emission line, Angstroms
double extl[ARB] # O: extinction evaluation array
int npts # I: size of evaluation/wave arrays
double rv # I: ratio of total to selective extinction
# Local variables:
double a, b # coefficients of extinction function
int pix # loop index
double x # wavelength in inverse microns
double y # function value
begin
do pix = 1, npts {
if (wave[pix] < EPSILONR || 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
# the Savage & Mathis law to 0. at 1/lambda = 0.
if (x < 0.29D+0) {
y = -2.94 * x**2
# 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
# Compute E(lambda-V)/E(B-V) <- formula for alternate normalization
# y = rv * (a - 1.) + b
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).
procedure lmc_redlaw (wave, extl, npts)
# Calling arguments:
double wave[ARB] #I: evaluation wavelength array
double extl[ARB] #O: extinction evaluation array
int npts #I: size of evaluation/wave arrays
# Local variables:
double delt # work variable
double extab[NTAB] # tabulated extinction at E(B-V)=1.
int i, pix # loop index
double val # intermediate value
double x # 1. / (wave, in microns)
double xtab[NTAB] # tabulated inverse wavelengths
define NTAB 7
# 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] < EPSILOND || IS_INDEFD(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_INT (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.0D+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).
procedure smc_redlaw (wave, extl, npts)
# Calling argument:
double wave[ARB] #I: wavelength of emission line, Angstroms
double extl[ARB] #O: extinction evaluation array
int npts #I: size of evaluation/wave arrays
# Local variables:
double extab[NTAB] # tabulated extinction function from Prevot
int i, pix # loop indexes
double val # generic
double x # wavelength in inverse microns
double xtab[NTAB] # tabulated inverse wavelengths from Prevot
define NTAB 33
# 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] < EPSILOND || IS_INDEFD(wave[pix]) )
call error (1, "SMC_REDLAW: Invalid wavelength")
# Convert wavelength in angstroms to 1/microns
x = 10000. / wave[pix]
x = min (x, 10.D+0)
# Linearly interpolate extinction law in 1/lam
for (i=2; i <= NTAB && x > xtab[i]; i=i+1)
;
val = LIN_INT (extab[i-1], extab[i], xtab[i-1], xtab[i], x)
# Renormalize extinction function to A(lambda)/A(V)
extl[pix] = val / 3.1D+0 + 1.
}
end