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
|
/*===========================================================================
Copyright (C) 2001 European Southern Observatory (ESO)
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of
the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public
License along with this program; if not, write to the Free
Software Foundation, Inc., 675 Massachusetss Ave, Cambridge,
MA 02139, USA.
Corresponding concerning ESO-MIDAS should be addressed as follows:
Internet e-mail: midas@eso.org
Postal address: European Southern Observatory
Data Management Division
Karl-Schwarzschild-Strasse 2
D 85748 Garching bei Muenchen
GERMANY
===========================================================================*/
/* Program : opterrors.c */
/* Author : G. Mulas - ITAL_FLAMES Consortium */
/* Date : */
/* */
/* Purpose : Missing */
/* */
/* */
/* Input: see interface */
/* */
/* Output: */
/* */
/* DRS Functions called: */
/* none */
/* */
/* Pseudocode: */
/* Missing */
/* */
/* Version : */
/* Last modification date: 2002/08/05 */
/* Who When Why Where */
/* AMo 02-08-05 Add header header */
/*-------------------------------------------------------------------------*/
/* C functions include files */
/* 020716 KB */
#ifdef HAVE_CONFIG_H
# include <config.h>
#endif
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
/* MIDAS include files */
#include <flames_midas_def.h>
/* FLAMES-UVES include files */
#include <flames_uves.h>
#include <flames_opterrors.h>
#include <flames_newmatrix.h>
flames_err opterrors(flames_frame *ScienceFrame, allflats *SingleFF,
orderpos *Order, int32_t ix, frame_mask **mask,
double **aa, double **xx, int32_t *fibrestosolve,
int32_t *orderstosolve, int32_t numslices,
int32_t arraysize)
{
/* char output[70]; */ /* only for testing purposes */
double **covariance;
/* some array buffers to speed up access to multidimensional arrays */
frame_data *vecsigmaflatj, *vecflatk, *vecflatp;
int32_t m, n, i, j, k, p, ilowj, ihighj, ilowk, ihighk, ilowp, ihighp;
int32_t ilow, ihigh, ilow2, ihigh2;
int32_t orderj, orderk, orderp;
int32_t framej, framek, framep;
int32_t fibrem, fibrej, fibrek, fibrep;
int32_t fibrelow, fibrehigh, fibrelow2, fibrehigh2, fibrelow3, fibrehigh3;
int32_t fibrelow4, fibrehigh4;
double overlap;
frame_data sigma=0, sigma2=0;
double dbuf1=0, dbuf2=0;
double *dvecbuf1=0;
double *dvecbuf2=0;
double *dvecbuf3=0;
int32_t *lvecbuf1=0;
int32_t *lvecbuf2=0;
frame_data *fdvecbuf1=0;
frame_data *fdvecbuf2=0;
frame_data *fdvecbuf3=0;
frame_data *fdvecbuf4=0;
frame_mask *fmvecbuf1=0;
int32_t moffset=0;
int32_t moffset1=0;
int32_t moffset2=0;
int32_t noffset=0;
int32_t noffset1=0;
int32_t noffset2=0;
int32_t mnindex1=0;
int32_t mnindex2=0;
int32_t joffset=0;
int32_t mjindex2=0;
int32_t njindex2=0;
int32_t koffset=0;
int32_t mkindex2=0;
int32_t nkindex2=0;
int32_t poffset=0;
int32_t mpindex2=0;
int32_t npindex2=0;
int32_t ordermfibremindex=0;
int32_t orderjfibrejixindex=0;
int32_t orderkfibrekixindex=0;
int32_t orderpfibrepixindex=0;
int32_t iixindex=0;
int32_t ksubcols;
/* Here the variance of the spectra is computed too. This calculation is
carried out at first order of approximation.
i.e. the following formulae hold
A*X=B
which gives the spectra and
(A+dA)*(X+dX)=B+dB
which with the approximation
dA*dX=0
and the first formula becomes
A*dX=dB+dA*X
and so we have the variances for the each spectrum.
Moreover, we assume the variance in the FF frames to be negligible with
respect to that in the Science Frame, which in turn means that dA is
negligible with respect to dB. If all of this holds, it can be shown
that, to a very decent first order approximation,
cov(Xi,Xj) = (Inv(A))ij
The above quantity is known already as a side product of the optimal
extraction previously done, and is taken as an input parameter here.
A more accurate estimate involves some heavier loops which add some
more complicated terms to this first estimate.
*/
covariance = dmatrix(1, numslices, 1, numslices);
dvecbuf1 = covariance[1];
dvecbuf2 = aa[1];
dvecbuf3 = xx[1]+1;
lvecbuf1 = SingleFF->lowfibrebounds[0][0]+ix;
lvecbuf2 = SingleFF->highfibrebounds[0][0]+ix;
fdvecbuf1 = ScienceFrame->frame_array[0]+ix;
fdvecbuf2 = ScienceFrame->frame_sigma[0]+ix;
fmvecbuf1 = mask[0]+ix;
fdvecbuf3 = ScienceFrame->specsigma[ix][0];
fdvecbuf4 = ScienceFrame->speccovar[ix][0];
ksubcols = ScienceFrame->subcols;
/* since the covariances matrix is symmetric, compute only the half
above (and including) the diagonal */
moffset1 = 0;
moffset2 = 0;
for (m=1; m<=numslices; m++)
{
moffset = m-1;
fibrelow = moffset>1 ? moffset : 1;
fibrehigh = (m+1)<numslices ? m+1 : numslices;
fibrem = fibrestosolve[m];
ordermfibremindex = (orderstosolve[m]*ScienceFrame->maxfibres)+fibrem;
/* this is overkill right now, since we only save the m,m+1 covariances,
but I prefer to write this in the Right Way now and comment out
superfluous stuff later */
/* after some profiling, it is clear that this part is a major bottleneck;
therefore, just do the excruciatingly detailed and expensive
calculations for the variances and the covariances of adjacent fibres.
Moreover, assume without checking that non-adjacent fibres do not
overlap */
noffset1 = (m-1)*numslices;
noffset2 = (m-1)*arraysize;
for (n=m; n<=numslices; n++)
{
noffset = n-1;
mnindex1 = moffset1+n;
mnindex2 = moffset2+n;
/* start with the first, presumably largest, component */
dvecbuf1[mnindex1] = dvecbuf2[mnindex2];
/* here check whether to spend additional time on this covariance */
if ((n>=fibrelow) && (n<=fibrehigh))
{
fibrelow2 = noffset>fibrelow ? noffset : fibrelow;
fibrehigh2 = (n+1)<fibrehigh ? n+1 : fibrehigh;
/* here goes the double internal loop for the second term */
for (j=fibrelow2; j<=fibrehigh2; j++)
{
joffset = j-1;
mjindex2 = moffset2+j;
njindex2 = noffset2+j;
fibrelow3 = joffset>fibrelow2 ? joffset : fibrelow2;
fibrehigh3 = (j+1)<fibrehigh2 ? j+1: fibrehigh2;
fibrej = fibrestosolve[j];
orderj = orderstosolve[j];
framej = SingleFF->fibre2frame[fibrej];
vecsigmaflatj = SingleFF->flatdata[framej].sigma[0]+ix;
orderjfibrejixindex =
(((orderj*SingleFF->maxfibres)+fibrej)*SingleFF->subcols);
ilowj = lvecbuf1[orderjfibrejixindex];
ihighj = lvecbuf2[orderjfibrejixindex];
for (k=fibrelow3; k<=fibrehigh3; k++)
{
koffset = k-1;
mkindex2 = moffset2+k;
nkindex2 = noffset2+k;
fibrelow4 = koffset>fibrelow3 ? koffset : fibrelow3;
fibrehigh4 = (k+1)<fibrehigh3 ? k+1 : fibrehigh3;
fibrek = fibrestosolve[k];
orderk = orderstosolve[k];
framek = SingleFF->fibre2frame[fibrek];
vecflatk = SingleFF->flatdata[framek].data[0]+ix;
/* this term can only be non zero if the jth and kth slices
have some overlap, i.e. if and only if the corresponding
fibres have some overlap */
orderkfibrekixindex =
(((orderk*SingleFF->maxfibres)+fibrek)*SingleFF->subcols);
ilowk = lvecbuf1[orderkfibrekixindex];
ihighk = lvecbuf2[orderkfibrekixindex];
/* set the y interval to the overlap between the jth and kth
slice */
ilow = (ilowj>ilowk) ? ilowj : ilowk;
ihigh = (ihighj<ihighk) ? ihighj : ihighk;
if (ihigh>=ilow)
{ /* there is some overlap */
overlap = 0;
iixindex = ilow*ksubcols;
for (i=ilow; i<=ihigh; i++)
{
if (fmvecbuf1[iixindex]==0)
{ /* this is a good pixel, include it */
sigma = fdvecbuf2[iixindex];
sigma2 = sigma*sigma;
overlap +=
(double) (fdvecbuf1[iixindex]*vecsigmaflatj[iixindex]*
vecflatk[iixindex]/sigma2);
}
iixindex += ksubcols;
} /* end of `i' loop */
dvecbuf1[mnindex1] += overlap*
(2*dvecbuf2[mjindex2]*dvecbuf2[njindex2]*dvecbuf3[koffset]+
(dvecbuf2[mjindex2]*dvecbuf2[nkindex2]+
dvecbuf2[mkindex2]*dvecbuf2[njindex2])*dvecbuf3[joffset]);
} /* end of if statement */
/* add here the third loop for last contribution to the
covariance */
for (p=fibrelow4; p<=fibrehigh4; p++)
{
poffset = p-1;
mpindex2 = moffset2+p;
npindex2 = noffset2+p;
fibrep = fibrestosolve[p];
orderp = orderstosolve[p];
framep = SingleFF->fibre2frame[fibrep];
vecflatp = SingleFF->flatdata[framep].data[0]+ix;
/* this term can only be non zero if the jth, kth and pth
slices have some overlap, i.e. if and only if the
corresponding fibres have some overlap */
orderpfibrepixindex =
(((orderp*SingleFF->maxfibres)+fibrep)*SingleFF->subcols);
ilowp = lvecbuf1[orderpfibrepixindex];
ihighp = lvecbuf2[orderpfibrepixindex];
/* set the y interval to the overlap between the jth, the
kth and the pth slice */
ilow2 = (ilowp>ilow) ? ilowp : ilow;
ihigh2 = (ihighp<ihigh) ? ihighp : ihigh;
if (ihigh2>=ilow2)
{ /* there is some overlap */
overlap = 0;
iixindex = ilow*ksubcols;
for (i=ilow; i<=ihigh; i++)
{
if (fmvecbuf1[iixindex]==0)
{ /* this is a good pixel, include it */
sigma = fdvecbuf2[iixindex];
sigma2 = sigma*sigma;
overlap +=
(double) (vecsigmaflatj[iixindex]*vecflatk[iixindex]*
vecflatp[iixindex]/sigma2);
}
iixindex += ksubcols;
} /* end of `i' loop */
dbuf1 = dvecbuf3[joffset];
dbuf2 = dbuf1*dbuf1;
dvecbuf1[mnindex1] += overlap*
(dvecbuf2[mjindex2]*dvecbuf2[njindex2]*
dvecbuf3[koffset]*dvecbuf3[poffset]+
(dvecbuf2[mjindex2]*dvecbuf2[npindex2]+
dvecbuf2[mpindex2]*dvecbuf2[njindex2])*
dvecbuf3[joffset]*dvecbuf3[koffset]+
dvecbuf2[mkindex2]*dvecbuf2[npindex2]*dbuf2);
}
} /* end `p' loop */
} /* end `k' loop */
} /* end `j' loop */
} /* end if statement */
noffset1 += numslices;
noffset2 += arraysize;
} /* end `n' loop */
/* put the variance where it belongs */
fdvecbuf3[ordermfibremindex] = (frame_data) dvecbuf1[moffset1+m];
/* if the following covariance is between this fibre and the following
one, put it in speccovar */
if (((m+1)<=numslices) && (fibrestosolve[m+1]==fibrem+1))
fdvecbuf4[ordermfibremindex] = (frame_data) dvecbuf1[moffset1+m+1];
moffset1 += numslices;
moffset2 += arraysize;
} /* end `m' loop */
/* now the covariances matrix is not used anymore, free it */
free_dmatrix(covariance, 1, numslices, 1, numslices);
return NOERR;
}
|
