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
|
include <math/curfit.h>
# T_HDSHIFT -- Task hdshift in the dtoi package. This task is provided
# to support Kormendy's method of combining related characteristic curves.
# A zero point shift in log exp unique to each set of spots is calculated and
# subtracted. A single curve is fit to the combined, shifted data in
# a separate task (hdfit).
procedure t_hdshift ()
pointer sp, de, fun, save, db, cv, ps_coeff, den, exp, cvf
int fd, rec, nsave, ncoeff, nvalues, i, nfile, junk
real a0, ref_a0
pointer ddb_map()
int clpopni(), ddb_locate(), ddb_geti(), cvstati(), strncmp(), clgfil()
begin
# Allocate space on stack for string buffers
call smark (sp)
call salloc (de, SZ_FNAME, TY_CHAR)
call salloc (cvf, SZ_FNAME, TY_CHAR)
call salloc (fun, SZ_FNAME, TY_CHAR)
# Get list of the database names. The curfit information is retrieved
# from the first file in the list, the list is then rewound.
fd = clpopni ("database")
junk = clgfil (fd, Memc[cvf], SZ_FNAME)
call clprew (fd)
# Get coefficients of common fit from cv_file
db = ddb_map (Memc[cvf], READ_ONLY)
rec = ddb_locate (db, "cv")
nsave = ddb_geti (db, rec, "save")
call salloc (save, nsave, TY_REAL)
call ddb_gar (db, rec, "save", Memr[save], nsave, nsave)
call ddb_gstr (db, rec, "function", Memc[fun], SZ_LINE)
call cvrestore (cv, Memr[save])
ncoeff = cvstati (cv, CVNCOEFF)
call salloc (ps_coeff, ncoeff, TY_REAL)
if (strncmp (Memc[fun], "power", 5) == 0)
call cvcoeff (cv, Memr[ps_coeff], ncoeff)
else
call cvpower (cv, Memr[ps_coeff], ncoeff)
do i = 1, ncoeff {
call eprintf ("%d %.7g\n")
call pargi (i)
call pargr (Memr[ps_coeff+i-1])
}
call ddb_unmap (db)
nfile = 0
while (clgfil (fd, Memc[de], SZ_FNAME) != EOF) {
db = ddb_map (Memc[de], READ_ONLY)
call hds_read (db, den, exp, nvalues)
call hds_calc (den, exp, nvalues, Memr[ps_coeff], ncoeff, a0)
nfile = nfile + 1
if (nfile == 1)
ref_a0 = a0
a0 = a0 - ref_a0
call printf ("file %s: subtracting zero point a0 = %.7g\n")
call pargstr (Memc[de])
call pargr (a0)
# Write new log exposure information to database
db = ddb_map (Memc[de], APPEND)
call hds_wdb (db, exp, nvalues, a0)
call mfree (den, TY_REAL)
call mfree (exp, TY_REAL)
call ddb_unmap (db)
}
call clpcls (fd)
call sfree (sp)
end
# HDS_READ -- Read the density and exposure values from the database file.
# The density above fog and log exposure values are returned, as well as
# the number of data pairs read.
procedure hds_read (db, den, exp, nvalues)
pointer db # Pointer to input database file
pointer den # Pointer to density array - returned
pointer exp # Pointer to exposure array - returned
int nvalues # Number of data pairs read - returned
real fog
int nden, nexp, rec
int ddb_locate(), ddb_geti()
real ddb_getr()
begin
# Get fog value to be subtracted from density
rec = ddb_locate (db, "fog")
fog = ddb_getr (db, rec, "density")
# Get density array
rec = ddb_locate (db, "density")
nden = ddb_geti (db, rec, "den_val")
call malloc (den, nden, TY_REAL)
call ddb_gar (db, rec, "den_val", Memr[den], nden, nden)
call asubkr (Memr[den], fog, Memr[den], nden)
# Get exposure array
rec = ddb_locate (db, "exposure")
nexp = ddb_geti (db, rec, "log_exp")
call malloc (exp, nexp, TY_REAL)
call ddb_gar (db, rec, "log_exp", Memr[exp], nexp, nexp)
nvalues = min (nden, nexp)
end
# HDS_CALC -- Calculate the individual shift, a0.
procedure hds_calc (den, exp, nvalues, ps_coeff, ncoeff, a0)
pointer den
pointer exp
int nvalues
real ps_coeff[ARB]
int ncoeff
real a0
int i
real yavg, ycalc, xavg
begin
# Calculate average density and log exposure values
xavg = 0.0
yavg = 0.0
do i = 1, nvalues {
xavg = xavg + Memr[den+i-1]
yavg = yavg + Memr[exp+i-1]
}
xavg = xavg / real (nvalues)
yavg = yavg / real (nvalues)
ycalc = 0.0
do i = 2, ncoeff
ycalc = ycalc + ps_coeff[i] * (xavg ** real (i-1))
# Subtraction yields the zero point shift in question
a0 = yavg - ycalc
end
# HDS_WDB -- Write shifted log exposure values to database.
procedure hds_wdb (db, exp, nvalues, a0)
pointer db # Pointer to database
pointer exp # Pointer to array of exposure values
int nvalues # Number of exposure values in sample
real a0 # Shift to be subtracted
pointer sp, expsub
begin
call smark (sp)
call salloc (expsub, nvalues, TY_REAL)
call ddb_ptime (db)
call ddb_prec (db, "exposure")
call eprintf ("a0 = %g\n")
call pargr (a0)
call asubkr (Memr[exp], a0, Memr[expsub], nvalues)
call ddb_par (db, "log_exp", Memr[expsub], nvalues)
call ddb_putr (db, "A0 shift", a0)
call sfree (sp)
end
|
