← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/feros/libsrc/lnmatch.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/feros/libsrc/lnmatch.c
 
1
/* 
 
2
 
 
3
   Pascal Ballester (ESO/Garching)
 
4
 
 
5
   lnmatch.c
 
6
 
 
7
   line identification
 
8
 
 
9
.VERSION
 
10
 090728         last modif
 
11
*/
 
12
 
 
13
/* system includes */
 
14
 
 
15
#include <stdio.h>
 
16
 
 
17
/* general Midas includes */
 
18
 
 
19
#include <midas_def.h>
 
20
 
 
21
/* FEROS specific includes */
 
22
 
 
23
#include <lnmatch.h>
 
24
 
 
25
int match
 
26
#ifdef __STDC__
 
27
(
 
28
 int verif, double linid[], double linwav[], double liny[], double lindif[],
 
29
 int nbrow, double lcat[], int nbrowcat, double alpha, double *stdres,
 
30
 double dnull, int reject[]
 
31
 )
 
32
#else
 
33
     (
 
34
      verif, linid, linwav, liny, lindif, nbrow, 
 
35
      lcat, nbrowcat, alpha, stdres, dnull, reject
 
36
      )
 
37
     
 
38
     double    lcat[], linid[], linwav[], liny[], lindif[], alpha;
 
39
     double    dnull, *stdres;
 
40
     int       reject[], verif, nbrow, nbrowcat;
 
41
#endif
 
42
     
 
43
{
 
44
  int             row, rowcat;
 
45
  int             ordref;
 
46
  int             rowlow, rowupp, rowref, index;
 
47
  int             nmatch=0, ldir, ldir_ref;
 
48
  int             not_bottom, not_top, debug=0, nbump=0;
 
49
  double          resref, reslow, resupp, residual;
 
50
  double          resref_abs, reslow_abs, resupp_abs, res_abs;
 
51
  double          lambda_ref, lambda_low, lambda_upp, lambda_next;
 
52
  double          minlin, mincat;
 
53
  double          lambda_lin, lambda_cat;
 
54
  char            text[100];
 
55
 
 
56
  /* Verification loop */
 
57
  /* If the verification parameter verif<0, the line catalog is explored to
 
58
     make sure that wavelengths are sorted by increasing value and no 
 
59
     wavelength is duplicated. The value returned is >0 if no violation 
 
60
     was found */
 
61
  
 
62
mincat = minlin = lambda_lin = lambda_low = 0.0;
 
63
resupp_abs = reslow_abs = reslow = 0.0;
 
64
rowlow = 0;
 
65
 
 
66
  if (verif < 0) {
 
67
    verif = 1;
 
68
    lambda_low = lcat[1];
 
69
    lambda_next = lcat[2];
 
70
    
 
71
    ldir_ref = (lambda_next > lambda_low) ? 1 : -1;
 
72
    
 
73
    for (row=2; row<=nbrowcat; row++) 
 
74
      {
 
75
        lambda_upp = lcat[row];
 
76
        ldir = (lambda_upp > lambda_low) ? 1 : -1;
 
77
        if (ldir != ldir_ref) {
 
78
          sprintf(text,"Warning: Column :WAVE of the line catalog is not sorted by increasing wavelength");
 
79
          SCTPUT(text);
 
80
          verif = -1;
 
81
        }
 
82
        if (lambda_low == lambda_upp) {
 
83
          sprintf(text,"Error: Column :WAVE of the line catalog contains duplicate wavelength : %f",lambda_low);
 
84
          SCTPUT(text);
 
85
          verif = -2;
 
86
        }
 
87
        lambda_low = lambda_upp;
 
88
      }
 
89
    return(verif);
 
90
  }
 
91
  
 
92
  /* Matching loop */
 
93
  
 
94
  rowcat  = 1;
 
95
  *stdres = 0.;
 
96
  
 
97
  /* Some assumptions are made concerning the input data :
 
98
     - Index of arrays are in the range 1 to N
 
99
     - All N values in the input arrays are valid
 
100
     - Line catalog is sorted by increasing wavelength and contain no duplicated lines 
 
101
  */
 
102
  
 
103
  for (row = 1; row <= nbrow; row++) 
 
104
    {
 
105
    
 
106
      reject[row] = NOT_PROCESSED; 
 
107
      linid[row]  = dnull;
 
108
      lindif[row] = dnull;
 
109
    
 
110
      /* Reads current wavelength to be identified and relative order number */
 
111
      lambda_lin = linwav[row];
 
112
      ordref = liny[row];
 
113
      if (debug == 1) {
 
114
        sprintf(text,"Reading line.tbl row = %d line = %f\n",row,lambda_lin);
 
115
        SCTPUT(text); 
 
116
      }
 
117
    
 
118
      if ((not_bottom = (row > 1 && liny[row-1] == ordref)))  { 
 
119
        rowlow = row-1;
 
120
        lambda_low = linwav[rowlow];
 
121
        reslow = lambda_lin - lambda_low;
 
122
        reslow_abs = (reslow < 0) ? -reslow : reslow; /* abs. value */
 
123
      }
 
124
      
 
125
      if ((not_top = (row < nbrow && liny[row+1] == ordref))) {
 
126
        rowupp = row+1;
 
127
        lambda_upp = linwav[rowupp];
 
128
        resupp = lambda_lin - lambda_upp;
 
129
        resupp_abs = (resupp < 0) ? -resupp : resupp; /* abs. value */
 
130
      }
 
131
      
 
132
      if (not_bottom && not_top)  
 
133
        minlin = (reslow_abs<resupp_abs) ? reslow_abs : resupp_abs;
 
134
      if (! not_bottom) minlin = resupp_abs;
 
135
      if (! not_top)    minlin = reslow_abs;
 
136
    
 
137
      /* Reads the line catalog at the first entry point and initializes 
 
138
         variables */
 
139
      if (debug == 1) {
 
140
        sprintf(text,"Reading catalog thar.tbl rowcat = %d\n",rowcat);
 
141
        SCTPUT(text); 
 
142
      }
 
143
    
 
144
      rowupp = rowcat;
 
145
      lambda_upp = lcat[rowupp];
 
146
      resupp = lambda_lin - lambda_upp;
 
147
      resupp_abs = (resupp < 0) ? -resupp : resupp; /* abs. value */
 
148
    
 
149
      if (rowcat > 1) {
 
150
        rowlow = rowcat - 1;
 
151
        lambda_low = lcat[rowlow];
 
152
        reslow = lambda_lin - lambda_low;
 
153
        reslow_abs = (reslow < 0) ? -reslow : reslow; /* abs. value */
 
154
      }
 
155
      else reslow_abs = resupp_abs + 1.;
 
156
    
 
157
 
 
158
      /* The line catalog is assumed to be sorted by increasing wavelength. 
 
159
         Depending on which direction is found the minimum residual, the 
 
160
         catalog will be searched further in this direction */
 
161
    
 
162
      if (reslow_abs < resupp_abs) {
 
163
        index = -1;
 
164
        resref = reslow;
 
165
        resref_abs = reslow_abs;
 
166
        lambda_ref = lambda_low;
 
167
        rowcat = rowlow;
 
168
      }
 
169
      else {
 
170
        index = 1;
 
171
        resref = resupp;
 
172
        resref_abs = resupp_abs;
 
173
        lambda_ref = lambda_upp;
 
174
        rowcat = rowupp;
 
175
      }
 
176
    
 
177
      res_abs = resref_abs;
 
178
      residual = resref;
 
179
      lambda_cat = lambda_ref;
 
180
    
 
181
      /* Explores the catalog in the direction indicated by index and finds the
 
182
         position of the minimum residual */
 
183
      while (res_abs <= resref_abs) {
 
184
        resref = residual;
 
185
        resref_abs = res_abs;
 
186
        lambda_ref = lambda_cat;
 
187
        rowcat += index;
 
188
        if (rowcat < 1 || rowcat > nbrowcat) {
 
189
          
 
190
          ++nbump;
 
191
          /* reject[row] = REACHED_LIMIT; */
 
192
          break;
 
193
        }
 
194
        else {
 
195
          lambda_cat = lcat[rowcat];
 
196
          if (debug == 1) {
 
197
            sprintf(text,"Matching loop at rowcat = %d lambda_cat = %f, lambda_lin = %f\n", rowcat,lambda_cat,lambda_lin);
 
198
            SCTPUT(text); 
 
199
          }
 
200
          residual = lambda_lin - lambda_cat;
 
201
          res_abs = (residual < 0) ?  -residual : residual; /* abs. value */
 
202
        }
 
203
      }
 
204
    
 
205
      rowref = rowcat - index;
 
206
      if (debug == 1) {
 
207
        sprintf(text,"Minimum found at %d, lambda_cat = %7.1f, lambda_lin = %7.1f, Residual = %4.1f\n",rowref, lambda_ref, lambda_lin, resref_abs);
 
208
        SCTPUT(text); 
 
209
      }
 
210
    
 
211
      if ((not_bottom = (rowref > 1))) { 
 
212
        rowlow = rowref-1;
 
213
        lambda_low = lcat[rowlow];
 
214
        reslow = lambda_ref - lambda_low;
 
215
        reslow_abs = (reslow < 0) ? -reslow : reslow; /* abs. value */
 
216
      }
 
217
    
 
218
      if ((not_top = (rowref < nbrowcat))) {
 
219
        rowupp = rowref+1;
 
220
        lambda_upp = lcat[rowupp];
 
221
        resupp = lambda_ref - lambda_upp;
 
222
        resupp_abs = (resupp < 0) ? -resupp : resupp; /* abs. value */
 
223
      }
 
224
    
 
225
      if (not_bottom && not_top)  
 
226
        mincat = (reslow_abs<resupp_abs) ? reslow_abs : resupp_abs;
 
227
      if (! not_bottom) mincat = resupp_abs;
 
228
      if (! not_top)    mincat = reslow_abs;
 
229
    
 
230
      /* Find minimum of all minima */
 
231
    
 
232
      reslow = (minlin < mincat) ? minlin : mincat; 
 
233
      reslow *= alpha;
 
234
      
 
235
      if (debug == 1) {
 
236
        sprintf(text,"Comparing line %f with catalog %f.\nResidual is %f and maximum residual is %f\n",lambda_ref,lambda_lin,resref_abs,reslow);
 
237
        SCTPUT(text); 
 
238
      }
 
239
    
 
240
      if (reject[row] == NOT_PROCESSED) {
 
241
        
 
242
        if (resref_abs < reslow ) {
 
243
          linid[row] = lambda_ref, lindif[row] = resref;
 
244
          reject[row] = SUCCESSFULL_MATCH, *stdres += resref_abs, nmatch++;
 
245
          if (debug == 1) {
 
246
            sprintf(text,"**** Identified line: cat = %f  line = %f  Residual = %f\n",lambda_ref,lambda_lin,resref);
 
247
            SCTPUT(text); 
 
248
          }
 
249
        }
 
250
        else   reject[row] = MATCH_FAILED;
 
251
      
 
252
      } /* End of if (reject[row] == NOT_PROCESSED) */
 
253
      
 
254
      if (debug == 1) {
 
255
        sprintf(text,"Row %d Order %d - %f  %f  %f  %f  %f\n",row,ordref,minlin,resref_abs, lambda_ref, lambda_lin,mincat);
 
256
        SCTPUT(text); 
 
257
      }
 
258
      
 
259
    
 
260
    }                           /* End of for */
 
261
  
 
262
  if (nmatch > 0)  *stdres /= (double) nmatch;
 
263
  
 
264
  return(nmatch);
 
265
 
 
266
}