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
|
include <gset.h>
include <error.h>
include <math/curfit.h>
include "sensfunc.h"
define NCURVE 50 # Number of vectors in curve
# SF_GRAPH -- Graph the desired data on the output graphics device.
# This entry procedure determines the graph types, sets the device viewports
# for each graph, and calls a procedure to make the graph.
procedure sf_graph (gp, stds, nstds, cv, wextn, extn, nextn, ecv)
pointer gp # Graphics structure
pointer stds[nstds] # Standard star data
int nstds # Number of standard stars
pointer cv # Sensitivity function curve
real wextn[nextn] # Extinction table wavelengths
real extn[nextn] # Extinction table values
int nextn # Number of values in extinction table
pointer ecv # Residual extinction curve.
int i, image, ngraphs, strlen()
real fa[8]
pointer sp, id, gio
data fa/0.,1.,1.,0.,0.,0.,1.,1./
begin
# Clear the graphs, write the title, and set the viewports based on
# the number of graphs.
call smark (sp)
call salloc (id, SZ_LINE, TY_CHAR)
call sysid (Memc[id], SZ_LINE)
gio = GP_GIO(gp)
call gclear (gio)
call gseti (gio, G_FACOLOR, 0)
call gseti (gio, G_NMINOR, 0)
ngraphs = strlen (GP_GRAPHS(gp,1))
switch (ngraphs) {
case 1:
call gseti (gio, G_WCS, 1)
call gsview (gio, .10, .97, .1, .9)
GP_SZMARK(gp) = 2.
GP_SZMDEL(gp) = 2.
case 2:
call gseti (gio, G_WCS, 1)
call gsview (gio, .10, .97, .55, .9)
call gseti (gio, G_WCS, 2)
call gsview (gio, .10, .97, .10, .45)
GP_SZMARK(gp) = 2.
GP_SZMDEL(gp) = 2.
case 3:
call gseti (gio, G_WCS, 1)
call gsview (gio, .10, .97, .55, .9)
call gseti (gio, G_WCS, 2)
call gsview (gio, .10, .50, .1, .45)
call gseti (gio, G_WCS, 3)
call gsview (gio, .57, .97, .1, .45)
GP_SZMARK(gp) = 2.
GP_SZMDEL(gp) = 2.
default:
call gseti (gio, G_WCS, 1)
call gsview (gio, .10, .50, .55, .9)
call gseti (gio, G_WCS, 2)
call gsview (gio, .57, .97, .55, .9)
call gseti (gio, G_WCS, 3)
call gsview (gio, .10, .50, .1, .45)
call gseti (gio, G_WCS, 4)
call gsview (gio, .57, .97, .1, .45)
GP_SZMARK(gp) = .01
GP_SZMDEL(gp) = .01
}
# For each graph select the viewport and call a procedure to make
# the graph.
image = 0
for (i = 1; GP_GRAPHS(gp,i) != EOS; i = i + 1) {
call gseti (gio, G_WCS, i)
if (i > 1)
call gfill (gio, fa, fa[5], 4, GF_SOLID)
switch (GP_GRAPHS(gp,i)) {
case 'a', 's', 'r':
call sf_data (stds, nstds, GP_GRAPHS(gp,i))
call sf_graph1 (gp, stds, nstds, GP_GRAPHS(gp,i))
case 'e':
call sf_egraph (gp, wextn, extn, nextn, ecv)
case 'c':
call sf_cgraph (gp, stds, nstds, cv)
case 'i', 'l':
if (GP_GRAPHS(gp,i) == 'i')
GP_LOG(gp) = NO
else
GP_LOG(gp) = YES
image = image + 1
iferr (call sf_image (gp, image, stds, nstds, cv, wextn, extn,
nextn, ecv))
call erract (EA_WARN)
}
}
call gseti (gio, G_WCS, 0)
call gtext (gio, 0.5, 1., Memc[id], "h=c,v=t,f=b")
call gtext (gio, 0.5, 0.97, GP_TITLE(gp), "h=c,v=t,f=b")
call sfree (sp)
end
# SF_GRAPH1 -- Make graph of sensitivity or residual sensitivity vs wavelength.
procedure sf_graph1 (gp, stds, nstds, graph)
pointer gp # Graphics structure
pointer stds[nstds] # Standard star data
int nstds # Number of standard stars
char graph # Graph type
int i, j, n, mark, mdel, cdel, color
real szmark, szmdel, ymin, ymax, y1, y2
pointer x, y, w, gio
begin
gio = GP_GIO(gp)
# Autoscale the included data in Y and set wavelength range.
j = 0
do i = 1, nstds {
if (STD_FLAG(stds[i]) != SF_INCLUDE)
next
j = j + 1
n = STD_NWAVES(stds[i])
x = STD_X(stds[i])
y = STD_Y(stds[i])
if (j == 1)
call alimr (Memr[y], n, ymin, ymax)
else {
call alimr (Memr[y], n, y1, y2)
ymin = min (ymin, y1)
ymax = max (ymax, y2)
}
}
y1 = 0.05 * (ymax - ymin)
ymin = ymin - y1
ymax = ymax + y1
# Draw axes and title based on type of graph.
switch (graph) {
case 'a':
call gswind (gio, GP_AIRMASS(gp,1), GP_AIRMASS(gp,2), ymin, ymax)
call glabax (gio, "Sensitivity Residuals vs Airmass", "", "")
case 's':
call gswind (gio, GP_WSTART(gp), GP_WEND(gp), ymin, ymax)
call glabax (gio, "Sensitivity vs Wavelength", "", "")
case 'r':
call gswind (gio, GP_WSTART(gp), GP_WEND(gp), ymin, ymax)
call glabax (gio, "Sensitivity Residuals vs Wavelength", "", "")
}
# Plot the data with appropriate mark types and sizes.
mdel = GP_MDEL(gp)
cdel = GP_CDEL(gp)
szmdel = GP_SZMDEL(gp)
szmark = GP_SZMARK(gp)
do i = 1, nstds {
if (STD_FLAG(stds[i]) != SF_INCLUDE)
next
if (i != nstds-1) {
mark = GP_MARK(gp)
color = GP_CMARK(gp)
} else {
mark = GP_MADD(gp)
color = GP_CADD(gp)
}
n = STD_NWAVES(stds[i])
x = STD_X(stds[i]) - 1
y = STD_Y(stds[i]) - 1
w = STD_WTS(stds[i]) - 1
do j = 1, n {
if (Memr[w+j] == 0.) {
call gseti (gio, G_PLCOLOR, cdel)
call gmark (gio, Memr[x+j], Memr[y+j], mdel, szmdel, szmdel)
} else {
call gseti (gio, G_PLCOLOR, color)
call gmark (gio, Memr[x+j], Memr[y+j], mark, szmark, szmark)
}
}
}
end
# SF_EGRAPH -- Graph extinction curves with and without residual extinction
# correction.
procedure sf_egraph (gp, wextn, extn, nextn, ecv)
pointer gp # Graphics structure
real wextn[nextn] # Extinction table wavelengths
real extn[nextn] # Extinction table values
int nextn # Number of extinction table values
pointer ecv # Residual extinction curve
int i, j
real xmin, xmax, dx, x, cveval()
pointer sp, ext, extnew, gio
begin
call smark (sp)
call salloc (ext, NCURVE, TY_REAL)
# Interpolate extinction table to a grid of wavelengths within
# the range of the sensitivity data.
gio = GP_GIO(gp)
xmin = GP_WSTART(gp)
xmax = GP_WEND(gp)
dx = (xmax - xmin) / (NCURVE - 1)
x = xmin
do i = 0, NCURVE-1 {
call intrp (1, wextn, extn, nextn, x, Memr[ext+i], j)
x = x + dx
}
call gswind (gio, xmin, xmax, INDEF, INDEF)
call gascale (gio, Memr[ext], NCURVE, 2)
# If there is a residual extinction curve determine a new extinction
# curve.
if (ecv != NULL) {
call salloc (extnew, NCURVE, TY_REAL)
call amovr (Memr[ext], Memr[extnew], NCURVE)
x = xmin
do i = 0, NCURVE-1 {
Memr[extnew+i] = Memr[extnew+i] + cveval (ecv, x)
x = x + dx
}
call grscale (gio, Memr[extnew], NCURVE, 2)
}
# Draw the axes and title and extinction curves.
call glabax (gio, "Extinction vs Wavelength", "", "")
call gseti (gio, G_PLCOLOR, GP_PLCOLOR(gp))
call gvline (gio, Memr[ext], NCURVE, xmin, xmax)
if (ecv != NULL) {
call gseti (gio, G_PLTYPE, 2)
call gseti (gio, G_PLCOLOR, GP_PLCOLOR(gp)+1)
call gvline (gio, Memr[extnew], NCURVE, xmin, xmax)
call gseti (gio, G_PLTYPE, 1)
}
call sfree (sp)
end
# SF_FITGRAPH -- Overplot the fitted sensitivity curve.
procedure sf_fitgraph (gp, cv)
pointer gp # Graphics structure
pointer cv # Sensitivity function curve
int i, j
real x1, x2, y1, y2, cveval()
pointer gio
begin
gio = GP_GIO(gp)
call gseti (gio, G_PLCOLOR, GP_PLCOLOR(gp))
# Only plot on sensitivity curve graph types.
for (j = 1; GP_GRAPHS(gp,j) != EOS; j = j + 1) {
if (GP_GRAPHS(gp,j) != 's')
next
call gseti (gio, G_WCS, j)
call ggwind (gio, x1, x2, y1, y2)
x2 = (x2 - x1) / NCURVE
y1 = cveval (cv, x1)
call gamove (gio, x1, y1)
do i = 1, NCURVE {
x1 = x1 + x2
y1 = cveval (cv, x1)
call gadraw (gio, x1, y1)
}
}
end
|
