← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/mos/src/moscalib.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/src/moscalib.c
 
1
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
 
2
/* .COPYRIGHT   (C) 1993-2011 European Southern Observatory     */
 
3
/* .IDENT       moscalib.c                                  */
 
4
/* .AUTHORS     Pascal Ballester (ESO/Garching)            */
 
5
/*              Cristian Levin   (ESO/La Silla)            */
 
6
/*              Sabine Moehler   (LSW)                     */
 
7
/* .KEYWORDS    Spectroscopy, Long-Slit, MOS               */
 
8
/* .PURPOSE     Wavelength Calibration for MOS Spectra     */
 
9
/* .VERSION     1.0  Package Creation  17-MAR-1993         */
 
10
/*              1.1  modified for MOS  20-AUG-1993         */
 
11
/*                                                         */
 
12
/*
 
13
   .INPUT/OUTPUT:
 
14
       IN_A         :  line table   (def. linpos.tbl)
 
15
       IN_B         :  line catalog 
 
16
       MOS          :  mos table (slitlet positions) (def. mos)
 
17
       OUT_A        :  output coefficients table (def. coerbr.tbl)
 
18
       INPUTC       :  mode (Ident, Linear or Fors)(Constant/Variable fit)
 
19
       INPUTD/D/1/6 :  start(1),start(2),step(1),step(2),startw,pixel
 
20
       INPUTR/R/1/3 :  alpha,maxdev,tol
 
21
       INPUTI/I/1/5 :  degx,degl,miniter,maxiter,ystart
 
22
       INPUTI(10)   :  Debug level
 
23
 
 
24
.VERSION
 
25
 
 
26
 111025         last modif
 
27
 
 
28
   -------------------------------------------------------
 
29
*/
 
30
 
 
31
 
 
32
 
 
33
#include <tbldef.h>
 
34
#include <midas_def.h>
 
35
#include <ok.h>
 
36
#include <moscalib.h>
 
37
#include <proto_mos.h>
 
38
#include <proto_longmos.h>
 
39
#include <math.h>
 
40
 
 
41
#include <stdio.h>
 
42
#include <ctype.h>
 
43
 
 
44
 
 
45
int       inull;   /* Representation of NULL values */
 
46
float     rnull;
 
47
double    dnull;
 
48
char      mode[2];
 
49
int       recall;
 
50
#define DEGY 3
 
51
#define DEGXY 2
 
52
/*---------------------declaration of prototypes-----------------------------*/
 
53
#ifdef __STDC__
 
54
int fit_select(double [], double [], double [], int, double, int [],
 
55
               double [], double [], int, double [], int, int);
 
56
double comp_dif (double [], double [], double [], int);
 
57
double mode_init(char, double [], double [], double [], int, int);
 
58
void auto_id( int, int, int, double [], int[], float [], int [], double [], 
 
59
             int [], int [], double, int [], int [], int *, double [], int);
 
60
void write_icol(int, int, int [], int, int []);
 
61
int read_select( int, int, int []);
 
62
void read_ident(double *, double *, int, double *, double *, int *);
 
63
int fit_select_2D(double [], double [], double [], double [], int, double, 
 
64
             int [], double [], double [], double [], int, double [], int);
 
65
void auto_2d( int, int, int, double [], int[], float [], int [], double [], 
 
66
             int [], int [], double, int [], int [], int *, double [], int);
 
67
void read_ident_2D(double *, double *, double *, int, double *, double *, 
 
68
             double *, int *);
 
69
#else
 
70
int fit_select();
 
71
double comp_dif ();
 
72
double mode_init();
 
73
void auto_id();
 
74
void write_icol();
 
75
int read_select();
 
76
void read_ident();
 
77
int fit_select_2D();
 
78
void auto_2d();
 
79
void read_ident_2D();
 
80
#endif
 
81
 
 
82
extern void set_zero(), mos_eval_disp();
 
83
 
 
84
 
 
85
 
 
86
 
 
87
int main()
 
88
{
 
89
  char   line[80], lincat[60], outtab[60];       /*mode[2]; */
 
90
  char   mod_save[1], mos[60], text[120];        /*Table names*/
 
91
  int     tidlin, tidcat, tidmos;                /* Table id numbers */
 
92
  int    nbcol, nbrow, nsort, allcol, allrow;    /* Descriptors of line */
 
93
  int    nbcolcat, nbrowcat, nsortcat, allcolcat, allrowcat; 
 
94
                                                 /* Descriptors of lincat */
 
95
  int    nbcolmos, nbrowmos, nsortmos, allcolmos, allrowmos; 
 
96
  int    catwav, linwav, linwavc, linslit, mosoff, mosslit;
 
97
  int    linx,  liny, lindif, numlin[8], numcat[3];  /* Column numbers */
 
98
  int    linrej, *reject, *select, *yrow, nyrow, nbsel;
 
99
  int    ipar[10], degx,  minter, maxter, ystart; 
 
100
  int    maxslit, prevslit=0, *slitrow, numslit, islit, islit1;
 
101
  int    actvals, kunit, null;
 
102
  int    i, j, k, slit_id, ystart_row, rowstart, rowend, direction;
 
103
  int    slitsel, iref, itest, disp;
 
104
  int    slit, y_index, row, shift_tol;
 
105
  int    tmps, ix, rownum, rownum1, idrow, nid, loop=0;
 
106
 
 
107
  float  rpar[10], alpha, maxdev, tol;                 /* Real parameters */
 
108
 
 
109
  double *colx, *off;  
 
110
  double mos_save[100], *coly, *yval, *slitval, *colslit; 
 
111
  double delta, mindelta, spar[11], pixel;
 
112
  double xpos[50], ident[50], dpar[11], start[2], step[2];
 
113
  double *xid, *yid, *lid, xide[50], lide[50], Xid[50], Lid[50];        
 
114
  double tmpo, xtmp[100], diffx;
 
115
               
 
116
 
 
117
mindelta = 0.0;
 
118
ystart_row = slit_id = 0;
 
119
 
 
120
  SCSPRO("moscalib");
 
121
 
 
122
/* Read name of table line in keyword IN_A, lincat in IN_B 
 
123
   and mode in INPUTC(1:2). Mode can be Ident, Linear, or Fors 
 
124
   and Constant or Variable linelist*/
 
125
 
 
126
  SCKGETC("IN_A",  1, 60, &actvals, line);
 
127
  SCKGETC("IN_B",  1, 60, &actvals, lincat);
 
128
  SCKGETC("IN_C",   1, 60, &actvals, mos);
 
129
  SCKGETC("OUT_A", 1, 60, &actvals, outtab);
 
130
  SCKGETC("INPUTC",1, 2,  &actvals, mode);
 
131
  
 
132
  SCKRDR("INPUTR", 1,  3, &actvals, rpar, &kunit, &null);
 
133
  SCKRDI("INPUTI", 1,  5, &actvals, ipar, &kunit, &null);
 
134
  SCKRDI("INPUTI", 10, 1, &actvals, &disp,&kunit, &null);  
 
135
                                                            /* Display level */
 
136
  SCKRDD("INPUTD", 1,  7, &actvals, dpar, &kunit, &null);
 
137
  mod_save[0] = toupper(mode[0]);
 
138
  if (toupper(mode[0]) == 'I')  
 
139
    {
 
140
      SCKRDD("XPOS",  1,  50, &actvals, xpos, &kunit, &null);
 
141
      SCKRDD("LID", 1,  50, &actvals, ident, &kunit, &null);
 
142
    }
 
143
 
 
144
  degx   = ipar[0];  
 
145
  minter = ipar[2]; 
 
146
  maxter = ipar[3]; 
 
147
  ystart = ipar[4];
 
148
 
 
149
  alpha = rpar[0]; 
 
150
  maxdev = rpar[1]; 
 
151
  tol = rpar[2];    
 
152
 
 
153
  start[0] = dpar[0]; 
 
154
  start[1] = dpar[1]; 
 
155
  step[0]  = dpar[2]; 
 
156
  step[1]  = dpar[3];
 
157
  dpar[10] = 0.;  /* linear dispersion */
 
158
 
 
159
 
 
160
/* Open table     mos     and search for column :xoffset. All 
 
161
   table are mapped to double precision arrays */
 
162
 
 
163
  if (TCTOPN(mos, F_I_MODE, &tidmos))
 
164
    SCTPUT("**** Error while opening table mos");
 
165
  TCIGET(tidmos, &nbcolmos, &nbrowmos, &nsortmos, &allcolmos, &allrowmos);
 
166
  numslit = nbrowmos;
 
167
/* Finds column :SLIT in table mos */
 
168
  TCCSER(tidmos, ":SLIT", &mosslit);
 
169
  if (mosslit == (-1)) SCTPUT("**** Column :SLIT not found");
 
170
/* Finds column :XOFFSET in table mos */
 
171
  TCCSER(tidmos, ":XOFFSET", &mosoff);
 
172
  if (mosoff == (-1)) SCTPUT("**** Column :XOFFSET not found");
 
173
 
 
174
/* Open table lincat and search for column :WAVE. All 
 
175
   table are mapped to double precision arrays */
 
176
 
 
177
  if (TCTOPN(lincat, F_I_MODE, &tidcat))
 
178
    SCTPUT("**** Error while opening line catalog");
 
179
  TCIGET(tidcat, &nbcolcat, &nbrowcat, &nsortcat, &allcolcat, &allrowcat);
 
180
    /* Finds column :WAVE in line catalog */
 
181
  TCCSER(tidcat, ":WAVE", &catwav);
 
182
  if (catwav == (-1)) SCTPUT("**** Column :WAVE not found");
 
183
  /* Storing adresses in array numcat[3] */
 
184
  numcat[0] = tidcat;
 
185
  numcat[1] = catwav;
 
186
  numcat[2] = nbrowcat;
 
187
 
 
188
/* Open table line and search for columns :X, :Y, :WAVE, :WAVEC, :RESIDUAL, 
 
189
    :REJECT. All table columns are mapped to double precision arrays */
 
190
 
 
191
  if (TCTOPN(line, F_IO_MODE, &tidlin))
 
192
    SCTPUT("**** Error while opening table line.tbl");
 
193
  TCIGET(tidlin, &nbcol, &nbrow, &nsort, &allcol, &allrow);
 
194
  TCCSER(tidlin, ":X", &linx);
 
195
  if (linx == (-1)) SCTPUT("**** Column :X not found");
 
196
  TCCSER(tidlin, ":Y", &liny);
 
197
  if (liny == (-1)) SCTPUT("**** Column :Y not found");
 
198
  TCCSER(tidlin, ":SLIT", &linslit);
 
199
  if (linslit == (-1)) SCTPUT("**** Column :SLIT not found");
 
200
  TCCSER(tidlin, ":WAVE", &linwav);
 
201
  if (linwav == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom",
 
202
                             "WAVE", &linwav);
 
203
  TCCSER(tidlin, ":WAVEC", &linwavc);
 
204
  if (linwavc == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom",
 
205
                              "WAVEC", &linwavc);
 
206
  TCCSER(tidlin, ":RESIDUAL", &lindif);
 
207
  if (lindif == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom",
 
208
                             "RESIDUAL", &lindif);
 
209
  TCCSER(tidlin, ":REJECT", &linrej);
 
210
  if (linrej == (-1)) TCCINI(tidlin, D_I4_FORMAT, 1, "I6", "Rejection Code",
 
211
                             "REJECT", &linrej);
 
212
    /* Storing adresses in array numlin[8] */
 
213
  numlin[0] = tidlin;
 
214
  numlin[1] = linwav;
 
215
  numlin[2] = linwavc;
 
216
  numlin[3] = lindif;
 
217
  numlin[4] = linrej;
 
218
  numlin[5] = linx;
 
219
  numlin[6] = liny;
 
220
  numlin[7] = nbrow;
 
221
 
 
222
/* Read other parameters: 
 
223
    INPUTI(1)   Degree in X
 
224
    INPUTI(2)   Degree in Y
 
225
    INPUTI(3)   Minimal number of iterations
 
226
    INPUTI(4)   Maximal number of iterations
 
227
   
 
228
    INPUTR(1)   Parameter Alpha (tolerance) for matching 
 
229
                (allowed residual for line identification = 
 
230
                alpha*(minimum residual)(between the matched and the 
 
231
                next catalog line OR between the searched and the next
 
232
                computed line))
 
233
    INPUTR(2)   Maximal standard deviation of residuals (pixels)
 
234
    INPUTR(3)   Tolerance of individual lines for fitting
 
235
                (pixels if >0, wavel. units if <0
 
236
*/
 
237
 
 
238
  if (disp >= 50) 
 
239
    {
 
240
      SCTPUT("   ");
 
241
      sprintf(text,"Line table                 : %s ",line); 
 
242
      SCTPUT(text);
 
243
      sprintf(text,"Line catalog               : %s ",lincat); 
 
244
      SCTPUT(text);
 
245
      sprintf(text,"Mode                       : %s ",mode); 
 
246
      SCTPUT(text);
 
247
      sprintf(text,"Nb of iterations (min,max) : %d , %d",minter,maxter); 
 
248
      SCTPUT(text);
 
249
      sprintf(text,"Degree                     : %d ",degx); 
 
250
      SCTPUT(text);
 
251
      sprintf(text,"Tolerance for individual lines (pixels) : %f",tol); 
 
252
      SCTPUT(text);
 
253
      sprintf(text,"Rejection parameter Alpha               : %f",alpha); 
 
254
      SCTPUT(text);
 
255
      sprintf(text,"Maximum mean deviation (pixels)         : %f",maxdev); 
 
256
      SCTPUT(text);
 
257
    }
 
258
 
 
259
/* Initialisations */
 
260
 
 
261
  /* Get representation of NULL values and allocate memory for auxiliary 
 
262
     arrays */
 
263
  TCMNUL(&inull, &rnull, &dnull);
 
264
  /* Load line table and line catalog in memory. Table columns are stored in 
 
265
     double although they can be of any type (I, R*4 or R*8) */
 
266
  select = (int *) osmmget((nbrow+1)*sizeof(int));
 
267
  nbsel = read_select(tidlin, nbrow, select);
 
268
  sprintf(text,"   Number of lines (total, selected) : %d, %d", nbrow, nbsel);
 
269
  SCTPUT(text);
 
270
  colx = (double *) osmmget((nbrow+1)*sizeof(double)); 
 
271
  read_col(tidlin, nbrow, linx, colx, dnull);  
 
272
  coly = (double *) osmmget((nbrow+1)*sizeof(double));
 
273
  read_col(tidlin, nbrow, liny,   coly, dnull);
 
274
  colslit = (double *) osmmget((nbrow+1)*sizeof(double));
 
275
  read_col(tidlin, nbrow, linslit, colslit, dnull);
 
276
  /* find out how many slitlets are selected */
 
277
  slitsel = 0;
 
278
  iref = -1;
 
279
  for (i=1; i<=nbsel; i++)
 
280
    {
 
281
      itest = colslit[i];
 
282
      if (itest != iref) 
 
283
        {
 
284
          slitsel++;
 
285
          iref = colslit[i];
 
286
        }
 
287
    }
 
288
  maxslit = slitsel;
 
289
  sprintf(text,"   number of slitlets                : %i",maxslit);
 
290
  SCTPUT(text);
 
291
  off = (double *) osmmget((nbrowmos+1)*sizeof(double)); 
 
292
  /* For table MOS the selection flag is ignored */
 
293
  /* Get offsets as function of slit number */
 
294
  for (i = 1; i <= nbrowmos; i++) off[i] = 0;
 
295
  for (i = 1; i <= nbrowmos; i++) 
 
296
    {
 
297
      TCERDD(tidmos, i, mosoff, &tmpo, &null); 
 
298
      TCERDI(tidmos, i, mosslit, &tmps, &null);
 
299
      off[tmps] = tmpo;   
 
300
    }
 
301
 
 
302
  slitval = (double *) osmmget((nbrow+1)*sizeof(double));
 
303
  yval    = (double *) osmmget((nbrow+1)*sizeof(double));
 
304
  reject  = (int *) osmmget((nbrow+1)*sizeof(int));
 
305
  slitrow = (int *) osmmget((nbrowmos+2)*sizeof(int));
 
306
  yrow    = (int *) osmmget((nbrow+1)*sizeof(int)); 
 
307
  xid     = (double *)   osmmget((nbrow+1)*sizeof(double));
 
308
  yid     = (double *)   osmmget((nbrow+1)*sizeof(double));
 
309
  lid     = (double *)   osmmget((nbrow+1)*sizeof(double));
 
310
 
 
311
  /* Go through the line table to identify the starting row number 
 
312
     for the different y values */
 
313
  y_index = 1, row=0;
 
314
  while (y_index <= nbsel) 
 
315
    {
 
316
      /* Aendern: yval direkt aus tidlin einlesen und dann nur wenn yval =
 
317
         coly yrow beschreiben. so muesset Einhalten der Selektion moeglich
 
318
         sein */
 
319
      yval[++row] = coly[y_index];
 
320
      yrow[row]   = y_index;
 
321
      while (y_index <= nbsel && coly[y_index] == yval[row]) ++y_index;
 
322
    }
 
323
  nyrow = row;
 
324
  yrow[nyrow+1] = nbsel+1;
 
325
 
 
326
  /* Go through the line table to identify the starting row number 
 
327
     for the different slitlets */
 
328
  row = 1;
 
329
  slit = 0;
 
330
  while (row <= nyrow) 
 
331
    {
 
332
      y_index = yrow[row];
 
333
      slitval[++slit] = colslit[y_index];
 
334
      iref = colslit[y_index];
 
335
      slitrow[iref] = row;  
 
336
      while (y_index <= nbsel && colslit[y_index] == slitval[slit]) 
 
337
        {
 
338
          ++row;
 
339
          y_index = yrow[row];
 
340
        }
 
341
    }
 
342
  slitval[slitsel+1] = iref+1; /* y-position of the selected 
 
343
                                  slitlets -               - TS*/
 
344
  slitrow[iref+1] = nyrow + 1; /* coerbr-table: row number in the table
 
345
                                  to write the coefs!     - TS*/
 
346
  /* Note: if the coefs are calculated for slitlet 1 in six rows, the data
 
347
     of the second slitlet starts in row slitrow[2] = 7   
 
348
     array counts from 1 to (number of slitlets + 1)      - TS*/
 
349
 
 
350
 
 
351
 
 
352
 
 
353
  /* For mode IDENT: Find the position ystart_row corresponding to the 
 
354
     position of ystart 
 
355
     For all other modes: Set ystart = first row of first selected slitlet */
 
356
  if (toupper(mode[0]) == 'I') 
 
357
    {
 
358
      /* ystart is in world coordiantes */
 
359
      for (row=1; row<=nyrow; row++) 
 
360
        {
 
361
          delta = yval[row] - ystart;
 
362
          if (delta < 0)  delta = delta*(-1.);
 
363
          if (row == 1)  mindelta = delta;
 
364
          if (delta <= mindelta)
 
365
            {
 
366
              mindelta = delta; 
 
367
              ystart_row = row;
 
368
              y_index = yrow[row];
 
369
              slit_id = colslit[y_index];
 
370
              off[0] = off[slit_id];
 
371
            }
 
372
        }
 
373
    }
 
374
  else 
 
375
    {
 
376
      ystart = coly[1];
 
377
      for (row=1; row<=nyrow; row++) 
 
378
        {
 
379
          delta = yval[row] - ystart;
 
380
          if (delta < 0)
 
381
            {
 
382
            delta = delta*(-1.);
 
383
            }
 
384
          if (row == 1)  mindelta = delta;
 
385
          if (delta <= mindelta)
 
386
            {
 
387
              mindelta = delta; 
 
388
              ystart_row = row;
 
389
              y_index = yrow[row];
 
390
              slit_id = colslit[y_index];
 
391
              off[0] = off[slit_id];
 
392
              dpar[6] = off[slit_id];
 
393
            }
 
394
        }
 
395
    }   
 
396
  /* Estimate a first dispersion relation in mode ident or read it
 
397
     in modes linear and fors */
 
398
  recall = 0;
 
399
  if (toupper(mode[1]) == 'T')
 
400
    setrefdeg_2D(degx,DEGY,DEGXY); /*DEGY=4 DEGXY=1 from the header*/
 
401
  pixel = mode_init(mode[0], xpos, ident, dpar, degx, nbrow); 
 
402
  if (pixel == -1) goto ciao;
 
403
  if (toupper(mode[0]) == 'R')  recall--;
 
404
 
 
405
  /* save display value */
 
406
  ipar[9] = disp;
 
407
  
 
408
  /* identify the selected slitlet only */ 
 
409
  for (islit = 1; islit <= maxslit; islit++) /*loop all slitlets to select*/
 
410
    {
 
411
      ipar[5] = slitval[islit];
 
412
      if (slitval[islit] == slit_id)   /*select here*/
 
413
        {
 
414
          slitsel--;
 
415
          for (direction=1; direction>=-1; direction-=2) /*two directions*/ 
 
416
            {
 
417
              if (direction == 1) 
 
418
                {
 
419
                  setrefdeg(degx);            
 
420
                  mos_initdisp(outtab,"NEW",1);
 
421
                  /* islit1:   ID-number of next slitlet  -TS0896*/
 
422
                  /* rowstart, rowend: 1st and last row in a slitlet */
 
423
                  rowstart = ystart_row;
 
424
                  islit1   = slitval[islit+1]; 
 
425
                  rowend   = slitrow[islit1]-1;
 
426
                  if (toupper(mode[1]) == 'V')
 
427
                    {
 
428
                      auto_id(rowstart, rowend, direction, yval, yrow, rpar, 
 
429
                            ipar, dpar, numlin, numcat, pixel, select, reject,
 
430
                            &loop, mos_save, nyrow);
 
431
                      finishdisp();  
 
432
                    }
 
433
                  else
 
434
                    {
 
435
                      auto_2d(rowstart, rowend, direction, yval, yrow, rpar, 
 
436
                            ipar, dpar, numlin, numcat, pixel, select, reject,
 
437
                            &loop, mos_save, nyrow);
 
438
                      setrefdeg(degx);            
 
439
                      finishdisp();  
 
440
                    }
 
441
                }
 
442
              if (direction == -1 && ystart_row > 1)  
 
443
                {
 
444
                  setrefdeg(degx);            
 
445
                  mos_initdisp(outtab,"OLD",1);
 
446
                  rowstart = ystart_row;
 
447
                  islit1   = slitval[islit];
 
448
                  rowend   = slitrow[islit1];
 
449
                  if (toupper(mode[1]) == 'V')
 
450
                    {
 
451
                      auto_id(rowstart, rowend, direction, yval, yrow, rpar, 
 
452
                              ipar, dpar, numlin, numcat, pixel, select, 
 
453
                              reject, &loop, mos_save, nyrow);
 
454
                      finishdisp();  
 
455
                    }
 
456
                  else
 
457
                    { /*here I need the mode_init to start with 1D-Coefs*/
 
458
                      pixel = mode_init(mode[0],xpos,ident,dpar,degx,nbrow); 
 
459
                      auto_2d(rowstart, rowend, direction, yval, yrow, rpar, 
 
460
                              ipar,dpar, numlin, numcat, pixel, select, reject,
 
461
                              &loop, mos_save, nyrow);
 
462
                      setrefdeg(degx);            
 
463
                      finishdisp();
 
464
                    }
 
465
                }
 
466
            } /* End of loop on direction */
 
467
        } /*end of (slitval[islit] == slit_id)*/
 
468
    } /*end of for all slitlets*/
 
469
  if (mos_save[1] > 0 && toupper(mode[0]) == 'R') 
 
470
    {
 
471
      recall = 1; /* successful fit -- mode 'R' is active now */
 
472
      /*  mos_savedisp(mos_save);*//* save the first dispersion coeficients*/
 
473
    }
 
474
 
 
475
 
 
476
 
 
477
 
 
478
 
 
479
 
 
480
  /* initialize spar[10] */
 
481
  /* For mode LINEAR the initial parameters dpar[10] are kept */
 
482
  for (i = 1; i <= 10; i++) spar[i] = dpar[i];
 
483
 
 
484
 
 
485
 
 
486
 
 
487
 
 
488
 
 
489
 
 
490
  for (i = 1; i <= maxslit; i++)  /* Loop for all selected slitlets */
 
491
    { /* i <= maxslit */
 
492
      islit = (int) slitval[i];
 
493
      ipar[5] = islit;
 
494
      if ( (islit != slit_id && slitsel > 0) || toupper(mode[0]) == 'I' )
 
495
        { /* islit != slit_id && slitsel > 0 -- slit_id has been done before*/
 
496
          slitsel--;
 
497
          spar[6] = off[islit]; 
 
498
          /* Estimate a first dispersion relation */
 
499
          if (toupper(mode[0]) == 'L')
 
500
            { /*linear no recall - the coefs are reinitialized*/  
 
501
              setdisp(degx,mos_save);
 
502
              pixel = mode_init(mode[0], xpos, ident, spar, degx, nbrow); 
 
503
              if (pixel == -1) goto ciao;
 
504
            } 
 
505
          if (toupper(mode[0]) == 'I' || mod_save[0] == 'I') /*interactive*/
 
506
            { /*interactive mode*/
 
507
              nid = 0;
 
508
              /* read identifications and correct for offset */
 
509
              for (j=0; j<50; j++) 
 
510
                {
 
511
                  if (ident[j] != 0)
 
512
                    {
 
513
                      Xid[nid] = xpos[j]-off[slit_id]+off[islit];
 
514
                      Lid[nid++] = ident[j];
 
515
                    }
 
516
                }
 
517
              /* read x-positions of lines in first row of resp. slitlet */
 
518
              ix=0;
 
519
              idrow = slitrow[islit];
 
520
              rownum = yrow[idrow];
 
521
              rownum1 = yrow[idrow+1];
 
522
              while (rownum < rownum1)
 
523
                xtmp[ix++] = colx[rownum++];
 
524
              /* identify lines */
 
525
              SCKRDI ("SHIFTTOL", 1, 1, &actvals, &shift_tol, &kunit,
 
526
                      &null);
 
527
              for (j = 0; j < 50; j++) lide[j] = 0;
 
528
              for (j = 0; j < 50; j++) xide[j] = 0;
 
529
              j = 0;
 
530
              k = 0;
 
531
              while (j < nid)
 
532
                {
 
533
                  for (k = 0; k < ix; k++)
 
534
                    {
 
535
                      diffx = xtmp[k]-Xid[j];
 
536
                      if (diffx < 0) diffx = -diffx;
 
537
                      if (diffx <= shift_tol)
 
538
                        {
 
539
                          lide[k] = Lid[j];
 
540
                          xide[k] = xtmp[k];
 
541
                        }
 
542
                    }
 
543
                  j++;
 
544
                }
 
545
              pixel = mode_init(mode[0], xide, lide, spar, degx, nbrow); 
 
546
            } /*end of interactive mode*/
 
547
 
 
548
          direction = 1; 
 
549
          setrefdeg(degx);
 
550
          mos_initdisp(outtab,"OLD",1);
 
551
          rowstart = slitrow[islit];
 
552
          islit1   = slitval[i+1];
 
553
          rowend   = slitrow[islit1]-1;   
 
554
          ipar[0] = degx;
 
555
          if (toupper(mode[1]) != 'V')
 
556
            {
 
557
              auto_2d(rowstart, rowend, direction, yval, yrow, rpar, 
 
558
                      ipar, spar, numlin, numcat, pixel, select, reject,
 
559
                      &loop, mos_save, nyrow);
 
560
              setrefdeg(degx);            
 
561
              finishdisp();
 
562
            }
 
563
          else
 
564
            {
 
565
              auto_id(rowstart, rowend, direction, yval, yrow, rpar, ipar, 
 
566
                      spar, numlin, numcat, pixel, select, reject, &loop, 
 
567
                      mos_save, nyrow);
 
568
              finishdisp();
 
569
            }
 
570
          prevslit = islit;
 
571
        } /*end of islit != slitid*/
 
572
      else  prevslit = slitval[i-1]; 
 
573
    } /*end of i<=maxslit*/ 
 
574
ciao:
 
575
 
 
576
  TCSINI(tidcat);
 
577
  TCTCLO(tidcat);
 
578
  TCSINI(tidlin);
 
579
  TCTCLO(tidlin);
 
580
 
 
581
  osmmfree((char *) off);
 
582
  
 
583
  osmmfree((char *)slitval );
 
584
  osmmfree((char *)yval    );
 
585
  osmmfree((char *)reject  );
 
586
  osmmfree((char *)slitrow );
 
587
  osmmfree((char *)yrow    );
 
588
  osmmfree((char *)select  );
 
589
  osmmfree((char *)colx    );
 
590
  osmmfree((char *)coly    );
 
591
  osmmfree((char *)colslit );
 
592
 
 
593
  osmmfree((char *) xid);
 
594
  osmmfree((char *) yid);
 
595
  osmmfree((char *) lid);
 
596
  SCSEPI();
 
597
return 0;
 
598
}
 
599
 
 
600
 
 
601
/*--------------------------------------------------------------------------*/
 
602
#ifdef __STDC__
 
603
double mode_init(char mod1, double x[], double id[], 
 
604
                 double linpar[], int degx, int nbrow)
 
605
#else
 
606
double mode_init( mod1, x, id, linpar, degx, nbrow )
 
607
  char   mod1;
 
608
  double x[],id[],linpar[];
 
609
  int    degx, nbrow;
 
610
#endif
 
611
 
 
612
/* Estimate a first dispersion relation in mode ident or read it
 
613
   in modes linear and fors */
 
614
 
 
615
{
 
616
   char   text[120];
 
617
   int    i, nid;
 
618
   double pixel, *xid, *lid, chi;
 
619
 
 
620
   xid     = (double *)   osmmget((nbrow+1)*sizeof(double));
 
621
   lid     = (double *)   osmmget((nbrow+1)*sizeof(double));
 
622
 
 
623
   switch (toupper(mod1)) 
 
624
     {
 
625
      case ('I') : 
 
626
        {
 
627
         nid = 0;
 
628
         for (i=0; i<50; i++) 
 
629
           {
 
630
           if (id[i] != 0)
 
631
              {
 
632
               xid[++nid] = x[i];
 
633
               lid[nid] = id[i];
 
634
              } 
 
635
           }
 
636
         if (nid < 2) 
 
637
           {
 
638
            sprintf (text,"Not enough identifications... Exiting.\n"); 
 
639
            SCTPUT(text); 
 
640
            goto ciao_m;    
 
641
           } 
 
642
         set_zero(degx);
 
643
         pixel = mos_fit_disp (&nid, &degx, xid, lid, &chi); 
 
644
         osmmfree((char *) xid);
 
645
         osmmfree((char *) lid);
 
646
         return (pixel);   
 
647
 
 
648
         break;
 
649
        }
 
650
 
 
651
      case ('L') : 
 
652
        {
 
653
          double coefs[2];  
 
654
          coefs[0] = linpar[4]-linpar[6]*linpar[5];/* x = lam0-offs*disp */ 
 
655
          coefs[1] = linpar[5];/*disp*/
 
656
          /*for (i=2; i<=6; i++) coefs[i] = 0;*//* these are degrees in Y*/
 
657
          setdisp(1,coefs);/* degree=1, coef[i] = coefs[i-1]*/
 
658
          pixel = linpar[5];   
 
659
          osmmfree((char *) xid);
 
660
          osmmfree((char *) lid);
 
661
          return (pixel);
 
662
          break;
 
663
        }
 
664
 
 
665
      case ('R') : 
 
666
        {
 
667
          double coefs[2];  
 
668
          coefs[0] = linpar[4]-linpar[6]*linpar[5];/* x = lam0-offs*disp */ 
 
669
          coefs[1] = linpar[5];/*disp*/
 
670
          /*for (i=2; i<=6; i++) coefs[i] = 0;*/
 
671
          setdisp(1,coefs); /* degree=1, coef[i] = coefs[i-1]*/
 
672
          pixel = linpar[5];   
 
673
          osmmfree((char *) xid);
 
674
          osmmfree((char *) lid);
 
675
          return (pixel);
 
676
          break;
 
677
        }
 
678
 
 
679
      default: 
 
680
        {
 
681
        osmmfree((char *) xid);
 
682
        osmmfree((char *) lid);
 
683
        (void) sprintf(text,
 
684
               "Error in moscalib.c: Unknown calibration method %c\n",mod1); 
 
685
         SCETER(9,text);                /* already leaves this function... */
 
686
         break;
 
687
        }
 
688
     }
 
689
ciao_m:
 
690
           pixel = -1.;
 
691
 
 
692
           osmmfree((char *) xid);
 
693
           osmmfree((char *) lid);
 
694
           return (pixel);
 
695
}
 
696
/*---------------------------------------------------------------------------*/
 
697
#ifdef __STDC__
 
698
void auto_id( int rstart, int rend, int direc, double val_y[], int row_y[], 
 
699
         float rlist[], int ilist[], double dlist[], int nlin[], 
 
700
         int ncat[], double pixel, int sel[], int rej[], 
 
701
         int *loop, double mos_save[], int nyrow)
 
702
#else
 
703
void auto_id( rstart, rend,  direc, val_y,  row_y, rlist, ilist, dlist, nlin, 
 
704
         ncat, pixel, sel, rej, loop, mos_save, nyrow)
 
705
  int rstart,rend,direc,row_y[],ilist[],nlin[],ncat[],sel[],rej[],*loop,nyrow;
 
706
  float rlist[];
 
707
  double val_y[],pixel,mos_save[],dlist[];  
 
708
#endif
 
709
 
 
710
{
 
711
  char   text[120];
 
712
  int    tidlin, linwavc, lindif, linrej, linwav, linx, liny, nbrow;
 
713
  int    *select, *sline, itmax, itmin, tidcat, nbrowcat, catwav, dgx;
 
714
  int    nblines, nbsel, y_index, iter, prevmatch, end_of_iter, slit;
 
715
  int    ungraceful_exit, linb=0, nmatch, verif=0, cal=0;
 
716
  int    lindeg = 1, nid, row, kunit, disp;
 
717
  double *colcat, *coldif, *colx, *coly, *colwav, *colwavc;
 
718
  double *xline, *yline, *xid, *lid, *xtmp, start[2], step[2];
 
719
  double slit_save[100], chi, current_y, stdres, stdpix, tolwav;
 
720
  float  alpha, maxdev, tol;
 
721
 
 
722
#ifdef __STDC__
 
723
  int match(int, double [], double [], double [], double [], int,
 
724
            double [], int, double, double *, double, int []);
 
725
#else
 
726
  int match();
 
727
#endif
 
728
 
 
729
/* Read parameters */
 
730
  dgx   = ilist[0];
 
731
  itmin = ilist[1];
 
732
  itmax = ilist[2];
 
733
  disp  = ilist[9];
 
734
  start[1] = dlist[1];
 
735
  step[1]  = dlist[3];
 
736
  alpha    = rlist[0];
 
737
  maxdev   = rlist[1];
 
738
  tol      = rlist[2];
 
739
 
 
740
 
 
741
  /* Get representation of NULL values and allocate memory for auxiliary 
 
742
     arrays */
 
743
  TCMNUL(&inull, &rnull, &dnull);
 
744
  
 
745
 
 
746
  /* read stored table adresses:  nlin = numlin  --  ncat = numcat*/
 
747
  tidlin  = nlin[0];
 
748
  linwav  = nlin[1]; 
 
749
  linwavc = nlin[2];
 
750
  lindif  = nlin[3];
 
751
  linrej  = nlin[4];
 
752
  linx    = nlin[5];
 
753
  liny    = nlin[6];
 
754
  nbrow   = nlin[7];
 
755
  tidcat  = ncat[0];
 
756
  catwav  = ncat[1];
 
757
  nbrowcat = ncat[2];
 
758
 
 
759
 
 
760
  /* Load line table and line catalog in memory. Table columns are stored in 
 
761
     double although they can be of any type (I, R*4 or R*8) */
 
762
  select = (int *) osmmget((nbrow+1)*sizeof(int));
 
763
  nbsel = read_select(tidlin, nbrow, select);
 
764
  colcat = (double *) osmmget((nbrowcat+1)*sizeof(double));
 
765
  read_col(tidcat, nbrowcat, catwav, colcat, dnull); 
 
766
  colx = (double *) osmmget((nbsel+1)*sizeof(double));
 
767
  read_col(tidlin, nbrow, linx,   colx, dnull);  
 
768
  coly = (double *) osmmget((nbsel+1)*sizeof(double));
 
769
  read_col(tidlin, nbrow, liny,   coly, dnull);
 
770
  coldif = (double *) osmmget((nbsel+1)*sizeof(double));
 
771
  read_col(tidlin, nbrow, lindif, coldif, dnull); 
 
772
  xline   = (double *)   osmmget((nbsel+1)*sizeof(double));
 
773
  yline   = (double *)   osmmget((nbsel+1)*sizeof(double));  
 
774
  sline   = (int *) osmmget((nbsel+1)*sizeof(int)); 
 
775
  colwav  = (double *)   osmmget((nbsel+1)*sizeof(double)); 
 
776
  colwavc = (double *)   osmmget((nbsel+1)*sizeof(double)); 
 
777
  xid     = (double *)   osmmget((nbsel+1)*sizeof(double));
 
778
  lid     = (double *)   osmmget((nbsel+1)*sizeof(double));
 
779
  rej  = (int *) osmmget((nbsel+1)*sizeof(int));
 
780
  xtmp     = (double *)   osmmget((nbsel+1)*sizeof(double)); /*xtmp=xline-off*/ 
 
781
 
 
782
 
 
783
  /* Write null values in all output columns */
 
784
  if (loop == 0) 
 
785
    {
 
786
      for (row=0; row<=nbrow; row++) 
 
787
        {
 
788
          sline[row] = row;
 
789
          colwav[row]  = dnull; 
 
790
          colwavc[row] = dnull;
 
791
          coldif[row]  = dnull;
 
792
          rej[row]  = inull;
 
793
        }
 
794
      write_dcol(tidlin, nbrow, sline, linwav, colwav);
 
795
      write_dcol(tidlin, nbrow, sline, linwavc,colwavc);
 
796
      write_dcol(tidlin, nbrow, sline, lindif, coldif);
 
797
      write_icol(tidlin, nbrow, sline, linrej, rej);
 
798
      *loop=1;
 
799
    }
 
800
 
 
801
 
 
802
  /* Read a set of values at first row */
 
803
  current_y = val_y[rstart];
 
804
  /* Read slitlet number */
 
805
  slit = ilist[5];
 
806
  
 
807
  nblines = 0;
 
808
  for (y_index=row_y[rstart]; y_index<row_y[rstart+1]; y_index++) 
 
809
    {
 
810
      if (colx[y_index] != dnull) /* check that line position is valid */
 
811
        {      
 
812
          xline[++nblines]  =  colx[y_index];
 
813
          yline[nblines]    =  coly[y_index];
 
814
          sline[nblines]   =  sel[y_index];
 
815
          xtmp[nblines]  =  colx[y_index]-dlist[6];
 
816
        }
 
817
    }
 
818
 
 
819
 
 
820
  /* matching lines -- Iteration Loop -- until the sample is unchaged*/
 
821
  iter=0, prevmatch=0;
 
822
  end_of_iter     = FALSE;
 
823
  ungraceful_exit = FALSE;
 
824
  if (disp >= 50)  
 
825
    {
 
826
      sprintf(text,"Variable dispersion for slit nr. %d ystart = %7.1f", 
 
827
              ilist[5], current_y); 
 
828
      SCTPUT(text);
 
829
    }
 
830
  while (!end_of_iter && !ungraceful_exit) 
 
831
    { /* iterative match loop */
 
832
      /* Compute wavelength for all lines: */
 
833
      if (recall != 0 && iter == 0)    /*recall disp... if (mode[0] = 'R')*/
 
834
        {
 
835
          if (recall == -1)          /* first linear fit for first slitlet*/
 
836
            {
 
837
              mos_eval_disp(xline, colwavc, nblines);
 
838
            }
 
839
          else                        /* recall disp... for other slitlets*/
 
840
            {
 
841
              setdisp(dgx,mos_save);
 
842
              mos_eval_disp(xtmp, colwavc, nblines);
 
843
            }
 
844
        }
 
845
      else 
 
846
        mos_eval_disp(xline, colwavc, nblines);
 
847
      /* match (nblines) lines against the catalog: the matched lines of the
 
848
       first row are recalled in the loop*/
 
849
      nmatch =  match(verif, colwav, colwavc, yline, coldif, nblines, 
 
850
                      colcat, nbrowcat, alpha, &stdres, dnull, rej);
 
851
      stdpix = stdres/pixel;
 
852
      if (disp >= 100)
 
853
        {
 
854
          sprintf(text,"   row Y = %4d: matching %2d lines out of %2d", (int)current_y, nmatch, nblines); 
 
855
          SCTPUT(text); 
 
856
        }
 
857
      iter++;
 
858
             /* Test end of iteration: */
 
859
      end_of_iter = ((iter >= itmax) || (iter > itmin &&  nmatch == prevmatch));
 
860
      prevmatch = nmatch;
 
861
      if (!(ungraceful_exit = (stdpix > maxdev)))       /*not crash condition*/
 
862
        {                 /* Read the matched lines  --  if !ungraceful_exit */
 
863
          read_ident(xline, colwav, nblines, xid, lid, &nid);
 
864
          pixel =  mos_fit_disp (&nid, &dgx, xid, lid, &chi);
 
865
          if (pixel < 0.) ungraceful_exit = TRUE;
 
866
        }
 
867
    } /* end_of_iter */
 
868
 
 
869
 
 
870
 
 
871
 
 
872
 
 
873
  /* Fit the single rows...      of LINPOS-table*/
 
874
  if (ungraceful_exit) 
 
875
    { /*ungraceful_exit*/
 
876
      sprintf(text,"Sorry, wrong identifications...\n"); 
 
877
      SCTPUT(text);
 
878
      set_zero(dgx);
 
879
      stdres = -1.;
 
880
      for (row=rstart; row <= rend+direc ; row+=direc) 
 
881
        {
 
882
          linb = (int) ((current_y - start[1])/step[1] + 1.5); 
 
883
          mos_writedisp (row, -1, linb, current_y, nyrow, stdres);
 
884
          cal = -1;
 
885
          SCKWRI("CAL", &cal, slit, 1, &kunit); 
 
886
        }
 
887
      if (recall != 0)
 
888
        { /*reset coefs to the start values*/ 
 
889
          setdisp(1,mos_save);
 
890
          setrefdeg(dgx);
 
891
        }
 
892
    }
 
893
  else 
 
894
    { /*not ungraceful_exit*/      
 
895
      mos_savedisp(slit_save); /* ... slitsave[i] = coef[i+1]*/
 
896
      cal = 1;
 
897
      SCKWRI("CAL", &cal, slit, 1, &kunit); 
 
898
      /* Fit dispersion coefficients for all rows: */
 
899
      /* changed from  !=  to  <= for the case rstart = rend  -- TS*/ 
 
900
      for (row=rstart; row != rend+direc ; row+=direc) 
 
901
        { /*loop on the rows*/  /* Read a set of values at constant y */
 
902
          current_y = val_y[row];
 
903
          nblines = 0;
 
904
          for (y_index=row_y[row]; y_index<row_y[row+1]; y_index++) 
 
905
            {
 
906
              /* do not check here for (colx[y_index] != dnull)if it's valid 
 
907
                 or your line match is mis-matched - fit_select does the job*/ 
 
908
              xline[++nblines]  =  colx[y_index];
 
909
              yline[nblines]    =  coly[y_index];
 
910
              sline[nblines]   =  sel[y_index];
 
911
            }
 
912
          setdisp(dgx, slit_save);        /*coef[i] = slt_save[i-1]*/
 
913
          mos_eval_disp(xline, colwavc, nblines);
 
914
          stdres = comp_dif(colwav, colwavc, coldif, nblines); 
 
915
          tolwav = (tol>0.) ? tol*pixel : -tol;
 
916
          nid = fit_select (xline, colwav, coldif, nblines, tolwav, 
 
917
                            rej, xid, lid, nid, colwavc, dgx, (int) current_y);
 
918
          if (nid >= 2)
 
919
            pixel = mos_fit_disp   (&nid, &dgx, xid, lid, &chi);
 
920
          if (pixel > 0 && nid >= 2)
 
921
            {
 
922
              mos_eval_disp(xline, colwavc, nblines);
 
923
              stdres = comp_dif(colwav, colwavc, coldif, nblines);
 
924
              if (disp >= 100) 
 
925
                {
 
926
                  sprintf(text,"      Final selection for Y = %4d: %2d lines out of %2d", (int)current_y, nid, nblines); 
 
927
                  SCTPUT(text); 
 
928
                }
 
929
              if (disp >= 50) 
 
930
                {
 
931
                  sprintf(text,"      RMS = %6.2f - Tolerance = %6.2f (wav. units)", stdres, tolwav); 
 
932
                  SCTPUT(text);
 
933
                  if (disp < 100) disp = 40; 
 
934
                }
 
935
              linb = (int) ((current_y - start[1])/step[1] + 1.5); 
 
936
              mos_writedisp (row, slit, linb, current_y, nyrow, stdres); 
 
937
            }
 
938
          if (pixel <= 0 || nid <= 1) 
 
939
            {
 
940
              stdres = -1.;
 
941
              set_zero(dgx);
 
942
              mos_writedisp (row, -1, linb, current_y, nyrow, stdres);
 
943
            }
 
944
          write_dcol(tidlin, nblines, sline, linwav, colwav);
 
945
          write_dcol(tidlin, nblines, sline, linwavc,colwavc);
 
946
          write_dcol(tidlin, nblines, sline, lindif, coldif);
 
947
          write_icol(tidlin, nblines, sline, linrej, rej);
 
948
          if (pixel > 0 && nid > dgx && recall == -1)
 
949
            {
 
950
              mos_savedisp(slit_save);
 
951
              mos_savedisp(mos_save);
 
952
              recall = recall + 2;
 
953
              if (disp >= 50) 
 
954
                {
 
955
                  printf("   save mos_disp : ");
 
956
                  printdisp();
 
957
                }
 
958
            }
 
959
        }  /* End of loop on row */
 
960
 
 
961
 
 
962
          /* Perform linear fit to get a mean linear dispersion for the 
 
963
             transformation*/
 
964
      if (pixel > 0 && nid >= dgx && recall == 0) 
 
965
        {
 
966
          mos_fit_disp   (&nid, &lindeg, xid, lid, &chi);
 
967
          mos_savedisp(mos_save); /* ...mos_save[i] = coef[i+1]*/ 
 
968
        }
 
969
      setrefdeg(dgx);            
 
970
    } /* end of (un)graceful_exit*/ 
 
971
 
 
972
  /*  mos_savedisp(mos_save); *//* ...mos_save[i] = coef[i+1]*/ 
 
973
 
 
974
  osmmfree((char *) select);
 
975
  osmmfree((char *) colcat);
 
976
  osmmfree((char *) colx);
 
977
  osmmfree((char *) coly);
 
978
  osmmfree((char *) coldif);
 
979
  osmmfree((char *) xline);
 
980
  osmmfree((char *) yline   );
 
981
  osmmfree((char *) sline   );
 
982
  osmmfree((char *) colwav  );
 
983
  osmmfree((char *) colwavc );
 
984
  osmmfree((char *) xid     );
 
985
  osmmfree((char *) lid     );
 
986
  osmmfree((char *) rej);
 
987
  osmmfree((char *) xtmp);  
 
988
}/* end of AUTO_ID */
 
989
 
 
990
/*---------------------------------------------------------------------------*/
 
991
#ifdef __STDC__
 
992
void auto_2d( int rstart, int rend, int direc, double val_y[], int row_y[], 
 
993
         float rlist[], int ilist[], double dlist[], int nlin[], 
 
994
         int ncat[], double pixel, int sel[], int rej[], 
 
995
         int *loop, double mos_save[], int nyrow)
 
996
#else
 
997
void auto_2d( rstart, rend,  direc, val_y,  row_y, rlist, ilist, dlist, nlin, 
 
998
         ncat, pixel, sel, rej, loop, mos_save, nyrow)
 
999
  int rstart,rend,direc,row_y[],ilist[],nlin[],ncat[],sel[],rej[],*loop,nyrow;
 
1000
  float rlist[];
 
1001
  double val_y[],pixel,mos_save[],dlist[];  
 
1002
#endif
 
1003
 
 
1004
{
 
1005
  char   text[120];
 
1006
  int    tidlin, linwavc, lindif, linrej, linwav, linx, liny, nbrow;
 
1007
  int    *select, *sline, itmax, itmin, tidcat, nbrowcat, catwav, dgx;
 
1008
  int    nblines, nbsel, y_index, iter, prevmatch, end_of_iter, slit;
 
1009
  int    ungraceful_exit, linb, nmatch, verif=0, cal=0;
 
1010
  int    lindeg = 1, i, nid, row=0, kunit, disp,ypixel[9999],numypix;
 
1011
  double *colcat, *coldif, *colx, *coly, *colwav, *colwavc;
 
1012
  double *xline, *yline, *xid, *yid, *lid, *xtmp, *ytmp, start[2], step[2];
 
1013
  double slit_save[100], chi, current_y, stdres, stdpix, tolwav;
 
1014
  float  alpha, maxdev, tol, offset;
 
1015
 
 
1016
#ifdef __STDC__
 
1017
  int match(int, double [], double [], double [], double [], int,
 
1018
            double [], int, double, double *, double, int []);
 
1019
#else
 
1020
  int match();
 
1021
#endif
 
1022
 
 
1023
/* Read parameters */
 
1024
  dgx   = ilist[0];
 
1025
  itmin = ilist[1];
 
1026
  itmax = ilist[2];
 
1027
  disp  = ilist[9];
 
1028
  start[1] = dlist[1];
 
1029
  step[1]  = dlist[3];
 
1030
  offset   = dlist[6];
 
1031
  alpha    = rlist[0];
 
1032
  maxdev   = rlist[1];
 
1033
  tol      = rlist[2];
 
1034
  
 
1035
 
 
1036
  /* Get representation of NULL values and allocate memory for auxiliary 
 
1037
     arrays */
 
1038
  TCMNUL(&inull, &rnull, &dnull);
 
1039
  /* read stored table adresses :  nlin=numlin ncat=numcat -- from main{}*/
 
1040
  tidlin  = nlin[0];
 
1041
  linwav  = nlin[1]; 
 
1042
  linwavc = nlin[2];
 
1043
  lindif  = nlin[3];
 
1044
  linrej  = nlin[4];
 
1045
  linx    = nlin[5];
 
1046
  liny    = nlin[6];
 
1047
  nbrow   = nlin[7];
 
1048
  tidcat  = ncat[0];
 
1049
  catwav  = ncat[1];
 
1050
  nbrowcat = ncat[2];
 
1051
 
 
1052
 
 
1053
  /* Load line table and line catalog in memory. Table columns are stored in 
 
1054
     double although they can be of any type (I, R*4 or R*8) */
 
1055
  select = (int *) osmmget((nbrow+1)*sizeof(int));
 
1056
  nbsel = read_select(tidlin, nbrow, select);
 
1057
  colcat = (double *) osmmget((nbrowcat+1)*sizeof(double));
 
1058
  read_col(tidcat, nbrowcat, catwav, colcat, dnull); 
 
1059
  colx = (double *) osmmget((nbsel+1)*sizeof(double));
 
1060
  read_col(tidlin, nbrow, linx,   colx, dnull);  
 
1061
  coly = (double *) osmmget((nbsel+1)*sizeof(double));
 
1062
  read_col(tidlin, nbrow, liny,   coly, dnull);
 
1063
  coldif = (double *) osmmget((nbsel+1)*sizeof(double));
 
1064
  read_col(tidlin, nbrow, lindif, coldif, dnull); 
 
1065
  xline   = (double *)   osmmget((nbsel+1)*sizeof(double));
 
1066
  yline   = (double *)   osmmget((nbsel+1)*sizeof(double));  
 
1067
  sline   = (int *) osmmget((nbsel+1)*sizeof(int)); 
 
1068
  colwav  = (double *)   osmmget((nbsel+1)*sizeof(double)); 
 
1069
  colwavc = (double *)   osmmget((nbsel+1)*sizeof(double)); 
 
1070
  xid     = (double *)   osmmget((nbsel+1)*sizeof(double));
 
1071
  yid     = (double *)   osmmget((nbsel+1)*sizeof(double));
 
1072
  lid     = (double *)   osmmget((nbsel+1)*sizeof(double));
 
1073
  rej  = (int *) osmmget((nbsel+1)*sizeof(int));
 
1074
  xtmp   = (double *)   osmmget((nbsel+1)*sizeof(double));
 
1075
  ytmp   = (double *)   osmmget((nbsel+1)*sizeof(double));
 
1076
 
 
1077
 
 
1078
  /* Write null values in all output columns */
 
1079
  if (loop == 0) 
 
1080
    {
 
1081
      for (row=0; row<=nbrow; row++) 
 
1082
        {
 
1083
          sline[row] = row;
 
1084
          colwav[row]  = dnull; 
 
1085
          colwavc[row] = dnull;
 
1086
          coldif[row]  = dnull;
 
1087
          rej[row]  = inull;
 
1088
        }
 
1089
      write_dcol(tidlin, nbrow, sline, linwav, colwav);
 
1090
      write_dcol(tidlin, nbrow, sline, linwavc,colwavc);
 
1091
      write_dcol(tidlin, nbrow, sline, lindif, coldif);
 
1092
      write_icol(tidlin, nbrow, sline, linrej, rej);
 
1093
      *loop=1;
 
1094
    }
 
1095
  current_y = val_y[rstart];   /* Read a set of values at first row */
 
1096
  slit = ilist[5];             /* Read slitlet number */
 
1097
  nblines = 0;
 
1098
  for (y_index=row_y[rstart]; y_index<row_y[rend+1]; y_index++) 
 
1099
    {
 
1100
      if (colx[y_index] != dnull)      /* check that line position is valid */
 
1101
        {       
 
1102
          xline[++nblines]  =  colx[y_index];
 
1103
          yline[nblines]    =  coly[y_index];
 
1104
          sline[nblines]    =  sel[y_index];
 
1105
          xtmp[nblines]  =  colx[y_index] - offset;
 
1106
          ytmp[nblines] = coly[y_index] - coly[row_y[rstart]] + coly[row_y[1]];
 
1107
          /*  ytmp[nblines] = coly[y_index] - coly[row_y[rstart]];*/
 
1108
        }
 
1109
    }
 
1110
  numypix = rend - rstart + 1;
 
1111
 
 
1112
  /* matching lines -- Iteration Loop -- until the sample is unchaged*/
 
1113
  iter=0, prevmatch=0;
 
1114
  end_of_iter     = FALSE;
 
1115
  ungraceful_exit = FALSE;
 
1116
  while (!end_of_iter && !ungraceful_exit) 
 
1117
    { /* iterative match loop */
 
1118
      if (toupper(mode[1]) == 'T' && numypix > 5) 
 
1119
        { /* the real 2D-fit: */
 
1120
          if (disp >= 50 && iter == 0)  
 
1121
            {
 
1122
              sprintf(text,"2D-Dispersion for slit nr. %d ystart = %7.1f numypix = %6d", ilist[5], current_y, numypix); 
 
1123
              SCTPUT(text);
 
1124
            }
 
1125
          if  (toupper(mode[0]) == 'I' && iter == 0)
 
1126
            {
 
1127
              mos_savedisp(mos_save);                       /*read 1-D coefs*/
 
1128
              for (i=dgx+1;i<=dgx+6;i++) mos_save[i] = 0.;
 
1129
              setdisp_2D(dgx,mos_save);                    /*write 2-D coefs*/
 
1130
            }
 
1131
          /* Compute wavelength for all lines: */
 
1132
          if (recall != 0 && iter == 0)  /*recall disp... if (mode[0] = 'F')*/
 
1133
            {
 
1134
              if (recall == -1)        /* first linear fit for first slitlet*/
 
1135
                {
 
1136
                  mos_eval_disp_2D(xtmp, yline, colwavc, nblines);
 
1137
                }
 
1138
              else                      /* recall disp... for other slitlets*/
 
1139
                {
 
1140
                  setdisp_2D(dgx,mos_save);
 
1141
                  mos_eval_disp_2D(xtmp, ytmp, colwavc, nblines);
 
1142
                }
 
1143
            }
 
1144
          else 
 
1145
            mos_eval_disp_2D(xline, yline, colwavc, nblines);
 
1146
 
 
1147
          /* match all lines against the catalog */
 
1148
          nmatch =  match(verif, colwav, colwavc, yline, coldif, nblines, 
 
1149
                          colcat, nbrowcat, alpha, &stdres, dnull, rej);
 
1150
          /*      for (i=1;i<=nblines;i++) printf("xline %6.1f yline[i] %7.1f colwavc %8.3f colwav %8.3f \n",xline[i],yline[i],colwavc[i],colwav[i]);*/
 
1151
          stdpix = stdres/pixel;
 
1152
          if (disp >= 100)
 
1153
            {
 
1154
              sprintf(text,"   row Y = %4d: matching %4d lines out of %4d", (int)current_y, nmatch, nblines); 
 
1155
              SCTPUT(text); 
 
1156
            }
 
1157
          /* Test end of iteration: */
 
1158
          iter++;
 
1159
          end_of_iter = ((iter >= itmax) || (iter > itmin &&  nmatch == prevmatch));
 
1160
          prevmatch = nmatch;
 
1161
          if (!(ungraceful_exit = (stdpix > maxdev))) 
 
1162
            {   /* Estimate the new dispersion relation on the matched lines */ 
 
1163
              read_ident_2D(xline, yline, colwav, nblines, xid, yid, lid,&nid);
 
1164
              pixel = mos_fit_disp_2D (&nid, &dgx, xid, yid, lid, &chi);
 
1165
              if (pixel < 0.) ungraceful_exit = TRUE;
 
1166
            }
 
1167
        }
 
1168
      else  
 
1169
        { /* constant dispersion over the slitlet */
 
1170
          if (disp >= 50 && iter == 0)  
 
1171
            {
 
1172
              sprintf(text,"constant-dispersion for slit nr. %d ystart = %7.1f numypix = %6d", ilist[5], current_y, numypix); 
 
1173
              SCTPUT(text);
 
1174
            }
 
1175
          /* Compute wavelength for all lines: */
 
1176
          if (recall != 0 && iter == 0) /*recall disp... if (mode[0] = 'F')*/
 
1177
            {
 
1178
              if (recall == -1) /* first linear fit for first slitlet*/
 
1179
                {
 
1180
                  mos_eval_disp(xtmp, colwavc, nblines);
 
1181
                }
 
1182
              else              /* recall disp... for other slitlets*/
 
1183
                {
 
1184
                  if (toupper(mode[1]) == 'T')
 
1185
                    { /* the coefs are still in 2D - but not if iter>0 */
 
1186
                      setdisp_2D(dgx, mos_save);
 
1187
                      mos_eval_disp_2D(xtmp, ytmp, colwavc, nblines);
 
1188
                    }
 
1189
                  else
 
1190
                    {
 
1191
                      setdisp(dgx,mos_save);
 
1192
                      mos_eval_disp(xtmp, colwavc, nblines);
 
1193
                    }
 
1194
                }
 
1195
            }
 
1196
          else 
 
1197
            mos_eval_disp(xline, colwavc, nblines);
 
1198
          /* match all lines against the catalog: */
 
1199
          nmatch =  match(verif, colwav, colwavc, yline, coldif, nblines, 
 
1200
                          colcat, nbrowcat, alpha, &stdres, dnull, rej);
 
1201
          stdpix = stdres/pixel;
 
1202
          if (disp >= 100)
 
1203
            {
 
1204
              sprintf(text,"   row Y = %4d: matching %4d lines out of %4d", (int)current_y, nmatch, nblines); 
 
1205
              SCTPUT(text); 
 
1206
            }
 
1207
          /* Test end of iteration: */
 
1208
          iter++;
 
1209
          end_of_iter = ((iter>=itmax) || (iter>itmin &&  nmatch==prevmatch));
 
1210
          prevmatch = nmatch;
 
1211
          if (!(ungraceful_exit = (stdpix > maxdev))) 
 
1212
            {    /* Estimate the new dispersion relation on the matched lines */
 
1213
              read_ident(xline, colwav, nblines, xid, lid,&nid);
 
1214
              pixel = mos_fit_disp(&nid, &dgx, xid, lid, &chi);
 
1215
              if (pixel < 0.) ungraceful_exit = TRUE;
 
1216
            }
 
1217
        }
 
1218
    } /* end_of_iter */
 
1219
 
 
1220
 
 
1221
 
 
1222
 
 
1223
 
 
1224
  /* Fit the data in two dimensions...     data from LINPOS-table*/
 
1225
  /* reject the lines with residuals greater then the limit*/
 
1226
  if (ungraceful_exit) 
 
1227
    { /*ungraceful_exit*/
 
1228
      sprintf(text,"\nSlitlet %3d: Sorry, wrong identifications...\n",ilist[5]); 
 
1229
      SCTPUT(text);
 
1230
      set_zero(dgx);
 
1231
      stdres = -1.;
 
1232
      for (row=rstart; row != rend+direc ; row+=direc) 
 
1233
        {
 
1234
          linb = (int) ((current_y - start[1])/step[1] + 1.5); 
 
1235
          mos_writedisp (row, -1, linb, val_y[row], nyrow, stdres);
 
1236
          cal = -1;
 
1237
          SCKWRI("CAL", &cal, slit, 1, &kunit); 
 
1238
        }
 
1239
      if (recall == -1)
 
1240
        { /*reset coefs to the start values*/ 
 
1241
          setdisp(1,mos_save);
 
1242
          setrefdeg(dgx);
 
1243
        }
 
1244
    }
 
1245
  else 
 
1246
    { /*(not un)graceful_exit*/      
 
1247
      /* Fit dispersion coefficients for all rows: */ 
 
1248
      mos_savedisp(slit_save);                /* ...slit_save[i] = coef[i+1]*/
 
1249
      cal = 1;
 
1250
      SCKWRI("CAL", &cal, slit, 1, &kunit); 
 
1251
      current_y = val_y[row];
 
1252
      nblines = 0; 
 
1253
      /* loop on the rows  -- Read values at all y for 2D-fit and constant: */
 
1254
      for (y_index=row_y[rstart]; y_index<row_y[rend+1]; y_index++) 
 
1255
        {                              /* check that line position is valid */
 
1256
          if (colx[y_index] != dnull)
 
1257
            {
 
1258
              xline[++nblines]  =  colx[y_index];
 
1259
              yline[nblines]    =  coly[y_index];
 
1260
              sline[nblines]   =  sel[y_index];
 
1261
            }
 
1262
        }
 
1263
      if (toupper(mode[1]) == 'T' && numypix > 5) 
 
1264
        { /* the real 2D-fit */
 
1265
          setdisp_2D(dgx,slit_save);              /* coef[i] = slit_save[i-1]*/
 
1266
          mos_eval_disp_2D(xline, yline, colwavc, nblines);
 
1267
          comp_dif(colwav, colwavc, coldif, nblines); 
 
1268
          tolwav = (tol>0.) ? tol*pixel : -tol;
 
1269
          nid = fit_select_2D (xline, yline, colwav, coldif, nblines, tolwav, 
 
1270
                               rej, xid, yid, lid, nid, colwavc, dgx);    
 
1271
          if (nid >= 2*numypix)
 
1272
            pixel = mos_fit_disp_2D(&nid, &dgx, xid, yid, lid, &chi);
 
1273
          if (pixel > 0 && nid >= 2*numypix)
 
1274
            {                           /* here we have achieved convergence*/
 
1275
              mos_eval_disp_2D(xline, yline, colwavc, nblines);
 
1276
              stdres = comp_dif(colwav, colwavc, coldif, nblines);
 
1277
              if (disp >= 100) 
 
1278
                {
 
1279
                  sprintf(text,"   Final selection for Y = %4d: %2d lines out of %2d", (int)current_y, nid, nblines); 
 
1280
                  SCTPUT(text); 
 
1281
                }
 
1282
              if (disp >= 50) 
 
1283
                {
 
1284
                  sprintf(text,"   Slit = %3d RMS = %6.2f - Tolerance = %6.2f (wav. units)", ilist[5], stdres, tolwav); 
 
1285
                  SCTPUT(text);
 
1286
                  if (disp < 100) disp = 40; 
 
1287
                }
 
1288
              for (row = rstart; row != rend+direc ; row += direc) 
 
1289
                {                                    /* save the coeficients*/
 
1290
                  ypixel[row] = (int) ((val_y[row]-start[1])/step[1] + 1.5); 
 
1291
                  mos_writedisp_2D(row,slit, ypixel[row],val_y[row],nyrow,stdres); 
 
1292
                }
 
1293
            }
 
1294
          if (pixel <= 0 || nid <= numypix) 
 
1295
            {                                     /* here the real desaster */
 
1296
              stdres = -1.;
 
1297
              set_zero(dgx);
 
1298
              for (row = rstart; row != rend+direc ; row += direc)
 
1299
                {                                  /* save crash coeficients*/
 
1300
                  ypixel[row] = (int) ((val_y[row]-start[1])/step[1] + 1.5); 
 
1301
                  mos_writedisp_2D (row, -1, ypixel[row], val_y[row], nyrow, stdres);
 
1302
                }
 
1303
              setdisp(1,mos_save);
 
1304
              setrefdeg(dgx);
 
1305
            }
 
1306
          write_dcol(tidlin, nblines, sline, linwav, colwav);
 
1307
          write_dcol(tidlin, nblines, sline, linwavc,colwavc);
 
1308
          write_dcol(tidlin, nblines, sline, lindif, coldif);
 
1309
          write_icol(tidlin, nblines, sline, linrej, rej);
 
1310
 
 
1311
          /* Perform linear fit to get a mean linear dispersion for the 
 
1312
             transformation or set the initial parameters of the recall: */
 
1313
          if (pixel > 0 && nid > dgx*numypix && toupper(mode[0]) == 'L')    
 
1314
            {                                    /* only for linear */
 
1315
              mos_fit_disp(&nid, &lindeg, xid, lid, &chi);
 
1316
              mos_savedisp(mos_save);
 
1317
            }                             /* method 2D-fit - set initial disp*/
 
1318
          if (pixel > 0 && nid >= dgx*numypix && recall == -1) 
 
1319
            {
 
1320
              mos_savedisp(mos_save);
 
1321
              recall = recall + 2;
 
1322
            }
 
1323
          setrefdeg_2D(dgx,DEGY,DEGXY);
 
1324
        } /* end of the real 2D-fit*/
 
1325
      else
 
1326
        { /* constant dispersion over the slitlet: */
 
1327
          setdisp(dgx,slit_save);/* coef[i] = slit_save[i-1]*/   
 
1328
          mos_eval_disp(xline, colwavc, nblines);
 
1329
          comp_dif(colwav, colwavc, coldif, nblines); 
 
1330
          tolwav = (tol>0.) ? tol*pixel : -tol;
 
1331
          nid = fit_select(xline, colwav, coldif, nblines, tolwav, 
 
1332
                           rej, xid, lid, nid, colwavc, dgx, (int) current_y);
 
1333
          if (nid >= 2*numypix)
 
1334
            pixel = mos_fit_disp(&nid, &dgx, xid, lid, &chi);
 
1335
          if (pixel > 0 && nid >= 2*numypix)
 
1336
            {
 
1337
              mos_eval_disp(xline, colwavc, nblines);
 
1338
              stdres = comp_dif(colwav, colwavc, coldif, nblines);
 
1339
              if (disp >= 100 && iter == 0) 
 
1340
                {
 
1341
                  sprintf(text,"   Final selection for Y = %4d: %2d lines out of %2d", (int)current_y, nid, nblines); 
 
1342
                  SCTPUT(text); 
 
1343
                }
 
1344
              if (disp >= 50) 
 
1345
                {
 
1346
                  sprintf(text,"   Slit = %3d  RMS = %6.2f - Tolerance = %6.2f (wav. units)", ilist[5], stdres, tolwav); 
 
1347
                  SCTPUT(text);
 
1348
                  if (disp < 100) disp = 40; 
 
1349
                }
 
1350
              for (row = rstart; row != rend+direc ; row += direc) 
 
1351
                {
 
1352
                  ypixel[row] = (int) ((val_y[row]-start[1])/step[1] + 1.5); 
 
1353
                  mos_writedisp(row,slit,ypixel[row],val_y[row],nyrow,stdres); 
 
1354
                }
 
1355
            }
 
1356
          if (pixel <= 0 || nid <= numypix) 
 
1357
            {
 
1358
              stdres = -1.;
 
1359
              set_zero(dgx);
 
1360
              for (row = rstart; row != rend+direc ; row += direc)
 
1361
                {
 
1362
                  ypixel[row] = (int) ((val_y[row]-start[1])/step[1] + 1.5);
 
1363
                  mos_writedisp(row, -1, ypixel[row],val_y[row],nyrow,stdres);
 
1364
                }
 
1365
            }
 
1366
          write_dcol(tidlin, nblines, sline, linwav, colwav);
 
1367
          write_dcol(tidlin, nblines, sline, linwavc,colwavc);
 
1368
          write_dcol(tidlin, nblines, sline, lindif, coldif);
 
1369
          write_icol(tidlin, nblines, sline, linrej, rej);
 
1370
          if (pixel > 0 && nid >= dgx*numypix && toupper(mode[0]) == 'L') 
 
1371
            {                                 /* method linear */ 
 
1372
              mos_fit_disp (&nid, &lindeg, xid, lid, &chi);
 
1373
              mos_savedisp(mos_save);
 
1374
            }
 
1375
          if (pixel > 0 && nid >= dgx*numypix && recall == -1)
 
1376
            {                         /* method constant - set initial coefs*/
 
1377
              mos_savedisp(mos_save);
 
1378
              recall = recall + 2;
 
1379
            }
 
1380
          setrefdeg(dgx);
 
1381
        }/* end of method constant*/
 
1382
    } /* end of (un)graceful_exit*/ 
 
1383
 
 
1384
  osmmfree((char *) select);
 
1385
  osmmfree((char *) colcat);
 
1386
  osmmfree((char *) colx);
 
1387
  osmmfree((char *) coly);
 
1388
  osmmfree((char *) coldif);
 
1389
  osmmfree((char *) xline);
 
1390
  osmmfree((char *) xtmp);
 
1391
  osmmfree((char *) ytmp);
 
1392
  osmmfree((char *) yline   );
 
1393
  osmmfree((char *) sline   );
 
1394
  osmmfree((char *) colwav  );
 
1395
  osmmfree((char *) colwavc );
 
1396
  osmmfree((char *) xid     );
 
1397
  osmmfree((char *) yid     );
 
1398
  osmmfree((char *) lid     );
 
1399
  osmmfree((char *) rej);
 
1400
  
 
1401
}/* end of AUTO_2D */
 
1402
/*---------------------------------------------------------------------------*/
 
1403
#ifdef __STDC__
 
1404
void write_icol(int tid, int nb, int select[], int colnb, int col[]) 
 
1405
#else
 
1406
void write_icol( tid, nb, select, colnb, col ) 
 
1407
  int tid,nb,select[],colnb,col[];
 
1408
#endif
 
1409
{
 
1410
  int i;
 
1411
 
 
1412
  for (i=1; i<=nb; i++)
 
1413
    TCEWRI(tid, select[i], colnb, &col[i]);
 
1414
}
 
1415
/*---------------------------------------------------------------------------*/
 
1416
#ifdef __STDC__
 
1417
double comp_dif ( double colwav[], double colwavc[], double coldif[], int nb)
 
1418
#else
 
1419
double comp_dif ( colwav, colwavc, coldif, nb)
 
1420
  double colwav[],colwavc[],coldif[];
 
1421
  int nb;
 
1422
#endif
 
1423
{
 
1424
  int  i, number;
 
1425
  double stdres;
 
1426
 
 
1427
  stdres = 0.;
 
1428
  number = 0;
 
1429
  for (i=1; i<=nb; i++) 
 
1430
    {
 
1431
     if (colwav[i] != dnull) 
 
1432
       {
 
1433
        coldif[i] = colwavc[i]-colwav[i];
 
1434
        number++;
 
1435
        stdres += coldif[i]*coldif[i];
 
1436
       }
 
1437
    }
 
1438
  stdres = sqrt(stdres/(double) number);
 
1439
  return stdres;
 
1440
}
 
1441
/*---------------------------------------------------------------------------*/
 
1442
#ifdef __STDC__
 
1443
int fit_select (double xline[], double colwav[], double coldif[], int nb, 
 
1444
             double tolwav, int reject[], double xid[], double lid[], 
 
1445
             int nid, double colwavc[], int dgx, int rowpos)
 
1446
#else
 
1447
int fit_select ( xline, colwav, coldif, nb, tolwav, reject, xid, lid, 
 
1448
                 nid, colwavc, dgx, rowpos)
 
1449
  double xline[],colwav[],coldif[],tolwav,xid[],lid[],colwavc[];
 
1450
  int    nb,reject[],nid,dgx,rowpos;
 
1451
#endif
 
1452
{
 
1453
  char      text[120];
 
1454
  int       maxpos=0, i=0;
 
1455
  double    px, maxres, resabs, chi, *tmpcolwav;
 
1456
  /* double    stdres; */
 
1457
 
 
1458
  /* tmpcolwav construction to reject lines temporary row by row only 
 
1459
     note that colwav is used for any row of a slitlet and the match would
 
1460
     be lost with any line rejected by fit_select              -- TS0896*/
 
1461
  tmpcolwav = (double *) osmmget ((nb+1)*sizeof(double));
 
1462
  for (i=1; i<= nb; i++) tmpcolwav[i] = colwav[i];
 
1463
  maxres = tolwav;
 
1464
  while (maxres >= tolwav) 
 
1465
    {
 
1466
      maxres = 0.;
 
1467
      nid = 0;
 
1468
      for (i=1; i<=nb; i++) 
 
1469
        {
 
1470
          if (reject[i] !=  RESIDUAL_GT_TOL && tmpcolwav[i] > 0)
 
1471
            {
 
1472
              resabs = (coldif[i] < 0) ? -coldif[i] : coldif[i];
 
1473
              if (resabs > maxres) maxpos = i, maxres = resabs;
 
1474
            }
 
1475
        }
 
1476
      if (maxres > tolwav)
 
1477
        {
 
1478
          if ( tmpcolwav[maxpos] > 0.)
 
1479
            {
 
1480
              sprintf(text,"   bad line at %10.3f - residual: %8.3f (wav. units)",  tmpcolwav[maxpos], maxres);
 
1481
              SCTPUT(text);
 
1482
            }
 
1483
          tmpcolwav[maxpos] = dnull;
 
1484
          reject[maxpos] = RESIDUAL_GT_TOL;
 
1485
          read_ident(xline, tmpcolwav, nb, xid, lid, &nid);
 
1486
          px = mos_fit_disp(&nid, &dgx, xid, lid, &chi); 
 
1487
          if (px > 0)
 
1488
            {
 
1489
              mos_eval_disp(xline, colwavc, nb);
 
1490
              comp_dif(tmpcolwav, colwavc, coldif, nb);
 
1491
              maxres=tolwav;
 
1492
            }     
 
1493
        }
 
1494
      else
 
1495
        {
 
1496
          for (i=1; i<=nb; i++)
 
1497
            {
 
1498
              if (reject[i] !=  RESIDUAL_GT_TOL && tmpcolwav[i] != dnull && xline[i])
 
1499
                {
 
1500
                  xid[++(nid)] = xline[i];
 
1501
                  lid[nid] = tmpcolwav[i];
 
1502
                }
 
1503
            } 
 
1504
        }
 
1505
    } /* end of while */
 
1506
  osmmfree((char *) tmpcolwav);
 
1507
  return nid;
 
1508
}
 
1509
/*---------------------------------------------------------------------------*/
 
1510
#ifdef __STDC__
 
1511
int fit_select_2D (double xline[],double yline[], double colwav[], 
 
1512
             double coldif[], int nb, double tolwav, int reject[],double xid[],
 
1513
             double yid[], double lid[], int nid, double colwavc[], int dgx)
 
1514
#else
 
1515
int fit_select_2D ( xline, yline, colwav, coldif, nb, tolwav, reject, xid, yid,
 
1516
                    lid, nid, colwavc, dgx) 
 
1517
  double xline[],yline[],colwav[],coldif[],tolwav,xid[],yid[],lid[],colwavc[];
 
1518
  int    nb,reject[],nid,dgx;
 
1519
#endif
 
1520
{
 
1521
  char      text[120];
 
1522
  int       maxpos=0, i=0;
 
1523
  double    px, maxres, resabs, chi, lidrej;
 
1524
  /* double    stdres; */
 
1525
 
 
1526
  maxres = tolwav;
 
1527
  while (maxres >= tolwav) 
 
1528
    {
 
1529
      maxres = 0.;
 
1530
      nid = 0;
 
1531
      lidrej = 0.;
 
1532
      for (i=1; i<=nb; i++) 
 
1533
        {
 
1534
          if (reject[i] !=  RESIDUAL_GT_TOL && colwav[i] != dnull)
 
1535
            {
 
1536
              resabs = (coldif[i] < 0) ? -coldif[i] : coldif[i];
 
1537
              if (resabs > maxres) maxpos = i, maxres = resabs;
 
1538
            }
 
1539
        }
 
1540
      if (maxres > tolwav)
 
1541
        {
 
1542
          /*warning for identified lines only: */                   
 
1543
          sprintf(text,"   bad line at %10.3f - residual: %8.3f (wav. units)",  colwav[maxpos], maxres);
 
1544
          SCTPUT(text);
 
1545
          colwav[maxpos] = dnull;
 
1546
          reject[maxpos] = RESIDUAL_GT_TOL;
 
1547
          read_ident_2D(xline, yline, colwav, nb, xid, yid, lid, 
 
1548
                        &nid);
 
1549
          px = mos_fit_disp_2D(&nid, &dgx, xid ,yid, lid, &chi); 
 
1550
          if (px > 0)
 
1551
            {
 
1552
              mos_eval_disp_2D(xline, yline, colwavc, nb);
 
1553
              comp_dif(colwav, colwavc, coldif, nb);
 
1554
              maxres = tolwav;
 
1555
            }
 
1556
        }
 
1557
      else
 
1558
        {      
 
1559
          for (i=1; i<=nb; i++) 
 
1560
            {
 
1561
              if (reject[i] !=  RESIDUAL_GT_TOL && colwav[i] != dnull)
 
1562
                {
 
1563
                  xid[++(nid)] = xline[i];
 
1564
                  yid[nid] = yline[i];
 
1565
                  lid[nid] = colwav[i];
 
1566
                }
 
1567
            }
 
1568
        } /* end of if */
 
1569
    } /* end of while */
 
1570
  return nid;
 
1571
} /* end of fit_select_2D */ 
 
1572
/*---------------------------------------------------------------------------*/
 
1573
#ifdef __STDC__
 
1574
int read_select(  int tid, int nb, int col[])
 
1575
#else
 
1576
int read_select( tid, nb, col )
 
1577
  int tid,nb,col[];
 
1578
#endif
 
1579
{
 
1580
  int i, k=0, sel;
 
1581
 
 
1582
  for (i=1; i<=nb; i++) 
 
1583
    {
 
1584
     TCSGET(tid, i, &sel);
 
1585
     if (sel)  col[++k] = i;
 
1586
    }
 
1587
  return(k);
 
1588
}
 
1589
/*---------------------------------------------------------------------------*/
 
1590
#ifdef __STDC__
 
1591
void read_ident(double *colx, double *colid, int nbsel, double *xid, double *lid, 
 
1592
           int *nid)
 
1593
#else
 
1594
void read_ident( colx, colid, nbsel, xid, lid, nid)
 
1595
  double *colx,*colid,*xid,*lid;
 
1596
  int    nbsel,*nid;
 
1597
#endif
 
1598
{
 
1599
 
 
1600
  int i;
 
1601
 
 
1602
  *nid = 0;
 
1603
 
 
1604
  for (i=1; i<=nbsel; i++) 
 
1605
    {
 
1606
     if (colid[i] != dnull && colx[i] !=dnull) 
 
1607
       {  
 
1608
        xid[++(*nid)]     = colx[i];
 
1609
        lid[*nid] = colid[i];
 
1610
       }
 
1611
    }
 
1612
 
 
1613
}/*---------------------------------------------------------------------------*/
 
1614
#ifdef __STDC__
 
1615
void read_ident_2D(double *colx, double *coly, double *colid, int nbsel, double *xid, 
 
1616
           double *yid, double *lid, int *nid)
 
1617
           
 
1618
#else
 
1619
void read_ident_2D( colx, coly, colid, nbsel, xid, yid, lid, nid)
 
1620
  double *colx,*coly,*colid,*xid,*yid,*lid;
 
1621
  int    nbsel,*nid;
 
1622
#endif
 
1623
{
 
1624
 
 
1625
  int i;
 
1626
 
 
1627
  *nid = 0;
 
1628
 
 
1629
  for (i=1; i<=nbsel; i++) 
 
1630
    {
 
1631
     if (colid[i] != dnull && colx[i] !=dnull) 
 
1632
       {  
 
1633
        xid[++(*nid)]     = colx[i];
 
1634
        yid[*nid]     = coly[i];
 
1635
        lid[*nid] = colid[i];
 
1636
       }
 
1637
    }
 
1638
}
 
1639
 
 
1640
 
 
1641
 
 
1642
 
 
1643
 
 
1644