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
|
include <error.h>
include <ctype.h>
include <mach.h>
include <math/curfit.h>
define NEPOINTS 100 # Number of points in extinction table
# SF_EOUT -- Output a revised extinction table. This is only done if there
# is at least one residual extinction curve. No assumption is made about
# overlapping extinction curves. In the overlap the extinction corrections
# are averaged.
procedure sf_eout (wextn, extn, nextn, ecvs, necvs)
real wextn[nextn] # Standard extinction wavelengths
real extn[nextn] # Standard extinction values
int nextn # Number of standard extinction points
pointer ecvs[necvs] # Residual extinction curves (one for each beam)
int necvs # Number of residual extinction curves
int i, j, fd, open(), scan(), nscan()
real w, ext, wmin, wmax, dw, xmin, xmax, cvstatr(), cveval()
pointer sp, fname, waves, extns, navg, cv
define newfile_ 99
begin
# If there are no residual extinction values then do nothing.
for (i=1; (i<=necvs) && (ecvs[i]==NULL); i=i+1)
;
if (i > necvs)
return
# The output table consists of NEPOINTS.
call smark (sp)
call salloc (fname, SZ_FNAME, TY_CHAR)
call salloc (waves, NEPOINTS, TY_REAL)
call salloc (extns, NEPOINTS, TY_REAL)
call salloc (navg, NEPOINTS, TY_INT)
call aclrr (Memr[extns], NEPOINTS)
call aclri (Memi[navg], NEPOINTS)
# Open the extinction table. If it fails allow the user to
# enter a new name.
call clgstr ("newextinction", Memc[fname], SZ_FNAME)
for (i=fname; (Memc[i]!=EOS) && IS_WHITE(Memc[i]); i=i+1)
if (Memc[i] == EOS) {
call sfree (sp)
return
}
newfile_
iferr (fd = open (Memc[i], NEW_FILE, TEXT_FILE)) {
call printf ("Cannot create %s -- Enter new extinction filename: ")
call pargstr (Memc[fname])
call flush (STDOUT)
if (scan() != EOF) {
call gargwrd (Memc[fname], SZ_FNAME)
if (nscan() == 1)
goto newfile_
}
call sfree (sp)
call printf ("No revised extinction file created\n")
return
}
# Determine the range of the extinction table.
wmin = MAX_REAL
wmax = -MAX_REAL
do i = 1, necvs {
if (ecvs[i] == NULL)
next
wmin = min (wmin, cvstatr (ecvs[i], CVXMIN))
wmax = max (wmax, cvstatr (ecvs[i], CVXMAX))
}
dw = (wmax - wmin) / (NEPOINTS - 1)
do i = 1, NEPOINTS
Memr[waves+i-1] = wmin + (i - 1) * dw
# Average the residual extinctions and add the original extinction.
do j = 1, necvs {
if (ecvs[j] == NULL)
next
cv = ecvs[j]
xmin = cvstatr (cv, CVXMIN)
xmax = cvstatr (cv, CVXMAX)
do i = 1, NEPOINTS {
w = Memr[waves+i-1]
if ((w < xmin) || (w > xmax))
next
Memr[extns+i-1] = Memr[extns+i-1] + cveval (cv, w)
Memi[navg+i-1] = Memi[navg+i-1] + 1
}
}
do i = 1, NEPOINTS {
if (Memi[navg+i-1] > 0)
Memr[extns+i-1] = Memr[extns+i-1] / Memi[navg+i-1]
w = Memr[waves+i-1]
call intrp (1, wextn, extn, nextn, w, ext, j)
Memr[extns+i-1] = Memr[extns+i-1] + ext
}
# Output extinction table.
call fprintf (fd, "# Revised extinction table.\n")
do i = 1, NEPOINTS {
call fprintf (fd, "%7.2f %6.3f\n")
call pargr (Memr[waves+i-1])
call pargr (Memr[extns+i-1])
}
call close (fd)
call sfree (sp)
end
|
