aboutsummaryrefslogtreecommitdiff
path: root/noao/onedspec/sensfunc/sfstds.x
blob: 072197295ad8642f0e425a8ce9e10ea90c78b85f (plain) (blame)
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
include "sensfunc.h"


# SF_STDS -- Get the standard observations for the specified apertures.
# If ignoring aperture set all apertures to 1.
# This routine knows the output of the task STANDARD.

procedure sf_stds (standards, aps, ignoreaps, stds, nstds)

char	standards	# Standard star data file
pointer	aps		# Aperture list
bool	ignoreaps	# Ignore apertures?
pointer	stds		# Pointer to standard observations (returned)
int	nstds		# Number of standard observations (returned)

int	i, j, fd, beam, npts, nwaves, nalloc
real	exptime, airmass, wstart, wend
real	wavelength, flux, dwave, count
pointer	sp, image, title, std
pointer	waves, fluxes, dwaves, counts, sens, fit, wts, iwts, x, y

bool	rng_elementi()
int	open(), fscan(), nscan(), stridxs()
errchk	open, malloc, realloc

begin
	call smark (sp)
	call salloc (image, SZ_STDIMAGE, TY_CHAR)
	call salloc (title, SZ_STDTITLE, TY_CHAR)

	# Open the standard observation data file.
	fd = open (standards, READ_ONLY, TEXT_FILE)

	# Read the standard observations and create a structure for each one.
	# The beginning of a new star is found by a line whose first word
	# begins with the character '['.  Otherwise the line is interpreted
	# as a data line.  All unrecognized formats are skipped.

	nwaves = 0
	nstds = 0
	while (fscan (fd) != EOF) {
	    call gargwrd (Memc[image], SZ_STDIMAGE)
	    if (Memc[image] == '[') {
	        call gargi (beam)
	        call gargi (npts)
		call gargr (exptime)
	        call gargr (airmass)
	        call gargr (wstart)
	        call gargr (wend)
	        call gargstr (Memc[title], SZ_STDTITLE)
	        if (nscan() < 7)
		    next
		if (!rng_elementi (aps, beam))
		    next
		if (IS_INDEF (exptime) || exptime <= 0.) {
		    call eprintf (
		"%s: Warning - exposure time missing or zero, using 1 second\n")
			call pargstr (Memc[image])
		    exptime = 1.
		}

		# For the first one create the pointer to the array of
		# structures.  For the following stars increase the size
		# of the pointer array and finish up the previous standard
		# star.

	        if (nstds == 0) {
		    nstds = nstds + 1
		    call calloc (stds, nstds, TY_INT)
		    call calloc (std, LEN_STD, TY_STRUCT)
		    Memi[stds+nstds-1] = std
	        } else {
	    	    if (nwaves > 0) {
			call realloc (waves, nwaves, TY_REAL)
			call realloc (fluxes, nwaves, TY_REAL)
			call realloc (dwaves, nwaves, TY_REAL)
			call realloc (counts, nwaves, TY_REAL)
			call realloc (wts, nwaves, TY_REAL)
			call malloc (sens, nwaves, TY_REAL)
			call malloc (fit, nwaves, TY_REAL)
			call malloc (iwts, nwaves, TY_REAL)
			call malloc (x, nwaves, TY_REAL)
			call malloc (y, nwaves, TY_REAL)
			call amovr (Memr[wts], Memr[iwts], nwaves)
			STD_NWAVES(std) = nwaves
			STD_WAVES(std) = waves
			STD_FLUXES(std) = fluxes
			STD_DWAVES(std) = dwaves
			STD_COUNTS(std) = counts
			STD_SENS(std) = sens
			STD_FIT(std) = fit
			STD_WTS(std) = wts
			STD_IWTS(std) = iwts
			STD_X(std) = x
			STD_Y(std) = y

			nstds = nstds + 1
			call realloc (stds, nstds, TY_INT)
			call calloc (std, LEN_STD, TY_STRUCT)
			Memi[stds+nstds-1] = std
		    }
	        }

		# Start a new standard star.
	        std = Memi[stds+nstds-1]
		if (ignoreaps)
	            STD_BEAM(std) = 1
		else
	            STD_BEAM(std) = beam
		STD_NPTS(std) = npts
	        STD_EXPTIME(std) = exptime
	        STD_AIRMASS(std) = airmass
		STD_WSTART(std) = wstart
		STD_WEND(std) = wend
		STD_SHIFT(std) = 0.
		STD_NWAVES(std) = 0

		# Decode the image and sky strings.
	        call strcpy (Memc[title], STD_TITLE(std), SZ_STDTITLE)
		i = stridxs ("]", Memc[image])
		if (Memc[image+i] == ']')
		    i = i + 1
		Memc[image+i-1] = EOS
	        call strcpy (Memc[image+1], STD_IMAGE(std), SZ_STDIMAGE)
		if (Memc[image+i] == '-') {
		    i = i + 2
		    j = stridxs ("]", Memc[image+i]) + i
		    Memc[image+j-1] = EOS
	            call strcpy (Memc[image+i], STD_SKY(std), SZ_STDIMAGE)
		} else
		    STD_SKY(std) = EOS
	        nwaves = 0

	    # Interprete the line as standard star wavelength point.
	    } else if (nstds > 0) {
		call reset_scan()
		call gargr (wavelength)
		call gargr (flux)
		call gargr (dwave)
		call gargr (count)
		if (nscan() < 3)
		    next
		if (wavelength < min (wstart, wend) ||
		    wavelength > max (wstart, wend) ||
		    flux<=0. || dwave<=0. || count<=0.)
		    next
		if (!rng_elementi (aps, beam))
		    next
		nwaves = nwaves + 1

		# Allocate in blocks to minimize the number of reallocations.
		if (nwaves == 1) {
		    nalloc = 100
		    call malloc (waves, nalloc, TY_REAL)
		    call malloc (fluxes, nalloc, TY_REAL)
		    call malloc (dwaves, nalloc, TY_REAL)
		    call malloc (counts, nalloc, TY_REAL)
		    call malloc (wts, nalloc, TY_REAL)
		} else if (nwaves > nalloc) {
		    nalloc = nalloc + 100
		    call realloc (waves, nalloc, TY_REAL)
		    call realloc (fluxes, nalloc, TY_REAL)
		    call realloc (dwaves, nalloc, TY_REAL)
		    call realloc (counts, nalloc, TY_REAL)
		    call realloc (wts, nalloc, TY_REAL)
		}

		# Record the data and compute the sensitivity.
		Memr[waves+nwaves-1] = wavelength
		Memr[fluxes+nwaves-1] = flux
		Memr[dwaves+nwaves-1] = dwave
		Memr[counts+nwaves-1] = count
		Memr[wts+nwaves-1] = 1.
	    }
	}

	# Finish up the last standard star and close the file.
	if (nstds > 0) {
	    STD_NWAVES(std) = nwaves
    	    if (nwaves > 0) {
		call realloc (waves, nwaves, TY_REAL)
		call realloc (fluxes, nwaves, TY_REAL)
		call realloc (dwaves, nwaves, TY_REAL)
		call realloc (counts, nwaves, TY_REAL)
		call realloc (wts, nwaves, TY_REAL)
		call malloc (sens, nwaves, TY_REAL)
		call malloc (fit, nwaves, TY_REAL)
		call malloc (iwts, nwaves, TY_REAL)
		call malloc (x, nwaves, TY_REAL)
		call malloc (y, nwaves, TY_REAL)
		call amovr (Memr[wts], Memr[iwts], nwaves)
		STD_WAVES(std) = waves
		STD_FLUXES(std) = fluxes
		STD_DWAVES(std) = dwaves
		STD_COUNTS(std) = counts
		STD_SENS(std) = sens
		STD_FIT(std) = fit
		STD_WTS(std) = wts
		STD_IWTS(std) = iwts
		STD_X(std) = x
		STD_Y(std) = y
	    }
	}
	call close (fd)
	call sfree (sp)

	# Add standard stars for any added points and composite points.
	nstds = nstds + 2
	call realloc (stds, nstds, TY_INT)
	call calloc (std, LEN_STD, TY_STRUCT)
	Memi[stds+nstds-2] = std
	call strcpy ("Added", STD_IMAGE(std), SZ_STDIMAGE)
	STD_BEAM(std) = STD_BEAM(Memi[stds])
	STD_NPTS(std) = STD_NPTS(Memi[stds])
	STD_EXPTIME(std) = 1.
	STD_AIRMASS(std) = 1.
	STD_WSTART(std) = STD_WSTART(Memi[stds])
	STD_WEND(std) = STD_WEND(Memi[stds])
	STD_SHIFT(std) = 0.
	STD_NWAVES(std) = 0
	call calloc (std, LEN_STD, TY_STRUCT)
	Memi[stds+nstds-1] = std
	call strcpy ("Composite", STD_IMAGE(std), SZ_STDIMAGE)
	STD_BEAM(std) = STD_BEAM(Memi[stds])
	STD_NPTS(std) = STD_NPTS(Memi[stds])
	STD_EXPTIME(std) = 1.
	STD_AIRMASS(std) = 1.
	STD_WSTART(std) = STD_WSTART(Memi[stds])
	STD_WEND(std) = STD_WEND(Memi[stds])
	STD_SHIFT(std) = 0.
	STD_NWAVES(std) = 0
end


# SF_FREE -- Free the standard observations and aperture array.

procedure sf_free (stds, nstds, apertures, napertures)

pointer	stds		# Pointer to standard observations
int	nstds		# Number of standard observations
pointer	apertures	# Pointer to apertures array
int	napertures	# Number of apertures

int	i
pointer	std

begin
	do i = 1, nstds {
	    std = Memi[stds+i-1]
	    if (STD_NWAVES(std) > 0) {
		call mfree (STD_WAVES(std), TY_REAL)
		call mfree (STD_FLUXES(std), TY_REAL)
		call mfree (STD_DWAVES(std), TY_REAL)
		call mfree (STD_COUNTS(std), TY_REAL)
		call mfree (STD_SENS(std), TY_REAL)
		call mfree (STD_FIT(std), TY_REAL)
		call mfree (STD_WTS(std), TY_REAL)
		call mfree (STD_IWTS(std), TY_REAL)
		call mfree (STD_X(std), TY_REAL)
		call mfree (STD_Y(std), TY_REAL)
	    }
	    call mfree (std, TY_STRUCT)
	}
	call mfree (stds, TY_INT)
	call mfree (apertures, TY_INT)
end