← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to stdred/mos/libsrc/mosdisp.c

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
stdred/mos/libsrc/mosdisp.c
 
1
/*===========================================================================
 
2
  Copyright (C) 1995-2009 European Southern Observatory (ESO)
 
3
 
 
4
  This program is free software; you can redistribute it and/or 
 
5
  modify it under the terms of the GNU General Public License as 
 
6
  published by the Free Software Foundation; either version 2 of 
 
7
  the License, or (at your option) any later version.
 
8
 
 
9
  This program is distributed in the hope that it will be useful,
 
10
  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
  GNU General Public License for more details.
 
13
 
 
14
  You should have received a copy of the GNU General Public 
 
15
  License along with this program; if not, write to the Free 
 
16
  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
17
  MA 02139, USA.
 
18
 
 
19
  Correspondence concerning ESO-MIDAS should be addressed as follows:
 
20
        Internet e-mail: midas@eso.org
 
21
        Postal address: European Southern Observatory
 
22
                        Data Management Division 
 
23
                        Karl-Schwarzschild-Strasse 2
 
24
                        D 85748 Garching bei Muenchen 
 
25
                        GERMANY
 
26
===========================================================================*/
 
27
 
 
28
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
 
29
/* .IDENT       mosdisp.c                                  */
 
30
/* .AUTHORS     Pascal Ballester (ESO/Garching)            */
 
31
/*              Cristian Levin   (ESO/La Silla)            */
 
32
/* .KEYWORDS    Spectroscopy, Long-Slit, MOS               */
 
33
/* .PURPOSE                                                */
 
34
/* .VERSION     1.0  Package Modification 21-SEP-1993      */
 
35
 
 
36
/*  090507              last modif                         */
 
37
/* ------------------------------------------------------- */
 
38
 
 
39
/* Purpose:
 
40
   This module is related to the dispersion relation and its representation.
 
41
   The upper modules (lncalib.c, lnrebin.c) do not know which representation is
 
42
   used (Chebyshev polynomials, or any other), nor how it is stored (which
 
43
   descriptor, meaning of coefficients). This design enables to replace 
 
44
   this module lndisp.c by any other based on a different representation.
 
45
 */
 
46
 
 
47
/* Note:
 
48
 
 
49
   Because the module is based on Numerical Recipes routines, which arrays
 
50
   start at index 1, the general rule is that arrays start at index 1. However
 
51
   the module can accomodate any starting index, which is specified at 
 
52
   initialisation time (function initdisp).
 
53
 */
 
54
 
 
55
/* Algorithms:
 
56
 
 
57
   The module provides ten top-level routines related to the evaluation of the
 
58
   dispersion relation. Only those six routines can be used in applications 
 
59
   related to the dispersion relation (wavelength calibration, rebinning, ...).
 
60
   The lower level routines must be considered as private.
 
61
 
 
62
   --- void mos_savedisp (save) ---
 
63
 
 
64
   Stores dispersion coefficients in array <save>.
 
65
 
 
66
 
 
67
   --- double mos_fit_disp (ndata, degree, x, l, chi) ---
 
68
 
 
69
   This routine evaluates the dispersion relation on the basis of ndata 
 
70
   couples of d values (x,l). Array indexes start at 0. The degree is 
 
71
   indicative. If ndata is too small, a smaller degree will be fitted 
 
72
   to the couples. Minimum value for ndata is 2. The procedure returns 
 
73
   a rough estimate of the average pixel size and the error of the 
 
74
   dispersion relation.
 
75
 
 
76
   --- void mos_eval_disp (x, l, ndata) ---
 
77
 
 
78
   Applies the last previously estimated (or loaded) dispersion relation 
 
79
   to ndata values of x and estimates the corresponding wavelength l. 
 
80
   Arrays start at index 0.
 
81
 
 
82
   --- void mos_initdisp (name, mode, start) ---
 
83
 
 
84
   Open (mode = OLD) or Create (mode = NEW) a storage table (named <name>).
 
85
   In addition to initdisp also the column :SLIT is created.
 
86
 
 
87
 
 
88
   --- void mos_writedisp(line, slit, pix_y, world_y, nyrow) ---
 
89
 
 
90
   Store in the storage table the dispersion relation corresponding to the 
 
91
   slitlet <slit>, row number <line>, and world coordinate <y>.
 
92
 
 
93
   --- double mos_readdisp(line) ---
 
94
 
 
95
   Retrieves from the storage table the dispersion relation corresponding 
 
96
   to the row number <line>.  The procedure returns a rough estimate of 
 
97
   the average pixel size.
 
98
 
 
99
 
 
100
 */
 
101
 
 
102
 
 
103
#include <stdio.h>
 
104
#include <ctype.h>
 
105
 
 
106
#include <dispersion.h>
 
107
 
 
108
 
 
109
 
 
110
/* -------------------------------------------------------------------------- */
 
111
 
 
112
#ifdef __STDC__
 
113
void 
 
114
mos_savedisp (double save[])
 
115
#else
 
116
void 
 
117
mos_savedisp (save)
 
118
     double save[];
 
119
#endif
 
120
 
 
121
{
 
122
 
 
123
  int i;
 
124
 
 
125
  /* printf("Dispersion Relation. Degree: %d. Refdeg: %d. MaxCoef:%d\n",
 
126
     fdeg,refdeg,maxcoef);
 
127
     printf("Coefficients: ");  */
 
128
 
 
129
  for (i = 0; i <= ncoef-1; i++)
 
130
    {
 
131
      save[i] = coef[i+1];
 
132
      /*     printf("save[%d] =  %f ",i, save[i]);  */
 
133
    }
 
134
 
 
135
}
 
136
 
 
137
/* -------------------------------------------------------------------------- */
 
138
 
 
139
#ifdef __STDC__
 
140
double 
 
141
mos_fit_disp (int *ndata, int *deg, double x[],
 
142
              double l[], double *chi)
 
143
#else
 
144
double 
 
145
mos_fit_disp (ndata, deg, x, l, chi)
 
146
     int *ndata, *deg;
 
147
     double x[], l[], *chi;
 
148
#endif
 
149
 
 
150
{
 
151
 
 
152
  double *sig, *chisq, **covar;
 
153
  int    *lista;
 
154
  int    actvals, maxdeg, i; 
 
155
  char   poltyp[12];
 
156
 
 
157
/* printf("NData : %d   ---   Deg value : %d\n",*ndata,*deg); */
 
158
 
 
159
  refdeg = *deg; 
 
160
  maxdeg = *ndata - 1;
 
161
  if (refdeg > maxdeg)
 
162
    fdeg = maxdeg;
 
163
  else
 
164
    fdeg = refdeg;
 
165
  ncoef = fdeg + 1;
 
166
  maxcoef = refdeg + 1;
 
167
 
 
168
  if (*ndata < 2)
 
169
    {
 
170
      printf ("Not enough lines (minimum is 2). \nNo dispersion relation computed\n");
 
171
      return (-2.);
 
172
    }
 
173
 
 
174
  else if (fdeg < 1)
 
175
    {
 
176
      printf ("Degree : %d. No dispersion relation fitted\n", *deg);
 
177
      return (-1.);
 
178
    }
 
179
 
 
180
  else
 
181
    {
 
182
      covar = (double **) dmatrix (1, *ndata, 1, *ndata);
 
183
      chisq = (double *) dvector (0, *ndata);
 
184
      sig = (double *) dvector (1, *ndata);
 
185
 
 
186
      lista = (int *) ivector (1, ncoef);
 
187
      for (i = 1; i <= ncoef; i++)
 
188
        lista[i] = i;
 
189
      for (i = 1; i <= *ndata; i++)
 
190
        sig[i] = 1.;
 
191
 
 
192
      SCKGETC("POLTYP", 1, 8, &actvals, poltyp);
 
193
      if (toupper(poltyp[0]) == 'L')
 
194
        lfit2(x,l,sig,*ndata,coef,ncoef,lista,ncoef,covar,chisq,fleg);
 
195
      else
 
196
        lfit2(x,l,sig,*ndata,coef,ncoef,lista,ncoef,covar,chisq,fcheb);
 
197
 
 
198
      *chi = *chisq;
 
199
 
 
200
      free_dmatrix (covar, 1, *ndata, 1, *ndata);
 
201
      free_dvector (chisq, 0, *ndata);
 
202
      free_dvector (sig, 1, *ndata);
 
203
      free_ivector (lista, 1, ncoef);
 
204
 
 
205
      pixbin = coef[2];
 
206
      return (pixbin);
 
207
    }
 
208
}
 
209
/* -------------------------------------------------------------------------- */
 
210
#ifdef __STDC__
 
211
void mos_eval_disp (double x[], double l[], int n)
 
212
#else
 
213
void mos_eval_disp (x, l, n)
 
214
     double x[], l[];
 
215
     int n;
 
216
#endif
 
217
 
 
218
{
 
219
 
 
220
  int actvals, i, icoef;
 
221
  double xp[MAXNCOE];
 
222
  char   poltyp[12];
 
223
 
 
224
  SCKGETC("POLTYP", 1, 8, &actvals, poltyp);
 
225
  for (i = start_index; i < (n + start_index); i++)
 
226
    {
 
227
      l[i] = 0.;
 
228
      if (toupper(poltyp[0]) == 'L')
 
229
        fleg (x[i], xp, ncoef);
 
230
      else 
 
231
        fcheb (x[i], xp, ncoef);
 
232
      for (icoef = 1; icoef <= ncoef; icoef++)
 
233
        l[i] += coef[icoef] * xp[icoef];
 
234
    }
 
235
}
 
236
 
 
237
/* -------------------------------------------------------------------------- */
 
238
 
 
239
#ifdef __STDC__
 
240
void mos_initdisp (char name[], char mode[], int start)
 
241
#else
 
242
void mos_initdisp (name, mode, start)
 
243
     char name[], mode[];
 
244
     int start;
 
245
#endif
 
246
 
 
247
{
 
248
 
 
249
  int icoef;
 
250
  char colnam[30], numb[10];
 
251
  int actvals, null;
 
252
  int kunit;
 
253
  int ncol, nrow, nsort, allcol, allrow;
 
254
 
 
255
 
 
256
  start_index = start;
 
257
 
 
258
  if (toupper (mode[0]) == 'N')
 
259
    {                           /*  NEW table */
 
260
      if (TCTINI (name, F_TRANS, F_IO_MODE, 5, 10, &tide))
 
261
        SCTPUT ("**** Error while creating output table");
 
262
      nbline = 0;
 
263
    }
 
264
  else
 
265
    {
 
266
      if (TCTOPN (name, F_IO_MODE, &tide))
 
267
        SCTPUT ("**** Error while opening output table");
 
268
      SCDRDD (tide, "LNPIX", 1, 1, &actvals, &pixbin, &kunit, &null);
 
269
      SCDRDI (tide, "LNDEG", 1, 1, &actvals, &refdeg, &kunit, &null);
 
270
      SCDRDI (tide, "LNCOE", 1, 1, &actvals, &maxcoef, &kunit, &null);
 
271
      /* fdeg=refdeg, ncoef=maxcoef; */
 
272
      TCIGET (tide, &ncol, &nrow, &nsort, &allcol, &allrow);
 
273
      nbline = nrow;
 
274
    }
 
275
 
 
276
  TCCSER (tide, ":SLIT", &colslit);
 
277
  if (colslit == (-1))
 
278
    TCCINI (tide, D_I4_FORMAT, 1, "I6", "Slit Number",
 
279
            "SLIT", &colslit);
 
280
 
 
281
  TCCSER (tide, ":ROW", &colline);
 
282
  if (colline == (-1))
 
283
    TCCINI (tide, D_I4_FORMAT, 1, "I6", "Row Number",
 
284
            "ROW", &colline);
 
285
 
 
286
  TCCSER (tide, ":Y", &coly);
 
287
  if (coly == (-1))
 
288
    TCCINI (tide, D_R8_FORMAT, 1, "F8.2", "Y Value",
 
289
            "Y", &coly);
 
290
 
 
291
  TCCSER (tide, ":RMS", &colrms);
 
292
  if (colrms == (-1))
 
293
    TCCINI (tide, D_R8_FORMAT, 1, "F8.4", "Angstrom",
 
294
            "RMS", &colrms);
 
295
 
 
296
  for (icoef = 1; icoef <= maxcoef; icoef++)
 
297
    {
 
298
      strcpy (colnam, ":COEF_");
 
299
      sprintf (numb, "%d", icoef);
 
300
      strcat (colnam, numb);
 
301
      TCCSER (tide, colnam, &colcoef[icoef]);
 
302
      if (colcoef[icoef] == (-1))
 
303
        TCCINI (tide, D_R8_FORMAT, 1, "F16.10", "Coefficients",
 
304
                colnam, &colcoef[icoef]);
 
305
    }
 
306
}
 
307
 
 
308
/*----------------------------------------------------------------------------*/
 
309
#ifdef __STDC__
 
310
int mos_readdisp (int y, int slit)
 
311
#else
 
312
int mos_readdisp (y, slit)
 
313
     int y, slit;
 
314
#endif
 
315
 
 
316
/* long int  y;     Pixel number of the row */
 
317
{
 
318
 
 
319
  int line, yval, linext, ydif, null, icoef;
 
320
  int mindif, actslit;
 
321
 
 
322
  mindif = -1;
 
323
 
 
324
  for (line = 1; line <= nbline; line++)
 
325
    {
 
326
      TCERDI (tide, line, colline, &yval, &null);
 
327
      TCERDI (tide, line, colslit, &actslit, &null);
 
328
      if (!null && actslit == slit)
 
329
        {
 
330
          ydif = y - yval;
 
331
          if (ydif < 0)
 
332
            ydif = (-1) * ydif;
 
333
          if (mindif < 0)
 
334
            mindif = ydif;
 
335
          if (ydif <= mindif)
 
336
            mindif = ydif, linext = line;
 
337
        }
 
338
    }
 
339
 
 
340
  if (mindif == -1)
 
341
    return (-1);
 
342
 
 
343
  fdeg = refdeg, ncoef = maxcoef;
 
344
 
 
345
  for (icoef = 1; icoef <= ncoef; icoef++)
 
346
    TCERDD (tide, linext, colcoef[icoef], &coef[icoef], &null);
 
347
 
 
348
  return (0);
 
349
 
 
350
}
 
351
 
 
352
/* -------------------------------------------------------------------------- */
 
353
 
 
354
#ifdef __STDC__
 
355
void mos_writedisp (int line, int slit, int ypix,
 
356
               double y, int numrow, double rms)
 
357
#else
 
358
void mos_writedisp (line, slit, ypix, y, numrow, rms)
 
359
     int line, slit, ypix, numrow;
 
360
     double y, rms;
 
361
#endif
 
362
 
 
363
{
 
364
 
 
365
  int icoef;
 
366
 
 
367
  TCEWRI (tide, line, colslit, &slit);
 
368
  TCEWRI (tide, line, colline, &ypix);
 
369
  TCEWRD (tide, line, coly, &y);
 
370
  TCEWRD (tide, line, colrms, &rms);
 
371
 
 
372
/*      if (nbline<line) nbline=line;   */
 
373
  if (nbline < line)
 
374
    nbline = numrow;
 
375
 
 
376
  for (icoef = 1; icoef <= maxcoef; icoef++)
 
377
    TCEWRD (tide, line, colcoef[icoef], &coef[icoef]);
 
378
}
 
379