← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/spec/src/spematch.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/spec/src/spematch.c
 
1
/*===========================================================================
 
2
  Copyright (C) 1991-2209 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       spematch.c                                       */
 
30
/* .AUTHOR      Pascal Ballester,  ESO - Garching                */
 
31
/* .KEYWORDS    Spectroscopy, Echelle,                           */
 
32
/* .PURPOSE     associates line catalog and computed wavelengths */
 
33
/* .VERSION     1.0    Creation    15-JAN-1991  PB               */
 
34
 
 
35
/* 090904               last modif                               */
 
36
/* ------------------------------------------------------------- */
 
37
 
 
38
#include <tbldef.h>
 
39
#include <midas_def.h>
 
40
#include <math.h>
 
41
#include <stdlib.h>
 
42
#include <stdio.h>
 
43
#define  ipos(col,row,nrow)   (col-1)*nrow + row - 1
 
44
 
 
45
void stdev(), open_acc(), close_acc(), accumulate(), findmax(); 
 
46
 
 
47
 
 
48
/* All Global Definitions for the accumulation array */
 
49
 
 
50
#define  ipos2D(col,row,npix) row*npix[0]+col
 
51
#define  MLEN 255
 
52
 
 
53
char     name[60];
 
54
char     Msg[MLEN];
 
55
float    *pntr, maxval;
 
56
int      imno, npix;
 
57
double   start, step, stend;
 
58
 
 
59
double          *pntr_buf;
 
60
int             nbrow;
 
61
 
 
62
int    tmp_pos[2000];
 
63
double tmp_iso[2000];
 
64
double tmp_res[2000];
 
65
 
 
66
/* Starting the main code */
 
67
 
 
68
int main()
 
69
{
 
70
 
 
71
   char            line[60], lincat[60];
 
72
 
 
73
   int             linwav, catwav, linid, lindif, liny, linerr, liniso;
 
74
   int             row, rowcat, null;
 
75
   int             order, ordref, bufsize, imno;
 
76
   int             status, kunit, actvals, tidlin, tidcat;
 
77
   int             nbcol, nsort, allcol, allrow;
 
78
   int             nbcolcat, nbrowcat, nsortcat, allcolcat, allrowcat;
 
79
   int             rowlow, rowupp, rowref, index, number;
 
80
 
 
81
   double          resref, reslow, resupp, residual, finres;
 
82
   double          resref_abs, reslow_abs, resupp_abs, res_abs;
 
83
   double          lambda_ref, lambda_low, lambda_upp;
 
84
   double          lambda_lin, lambda_cat, error, error_min;
 
85
 
 
86
   float           thres[3], fraction;
 
87
 
 
88
 
 
89
 
 
90
   SCSPRO("echmatch");
 
91
   rowlow = rowupp = 0;
 
92
   resupp = 0.0;
 
93
   
 
94
   status = SCKRDR("INPUTR", 1, 1, &actvals, thres, &kunit, &null);
 
95
   status = SCKGETC("IN_A", 1, 60, &actvals, line);
 
96
   status = SCKGETC("IN_B", 1, 60, &actvals, lincat);
 
97
 
 
98
   /* Read all keywords for the accumulator */
 
99
   status = SCKGETC("OUT_A",1, 60, &actvals, name);
 
100
   status = SCKRDI("INPUTI", 1, 1, &actvals, &npix, &kunit, &null);
 
101
   status = SCKRDD("INPUTD", 1, 1, &actvals, &start,&kunit, &null);
 
102
   status = SCKRDD("INPUTD", 2, 1, &actvals, &step, &kunit, &null);
 
103
   status = SCKRDR("INPUTR", 2, 1, &actvals, &fraction, &kunit, &null);
 
104
 
 
105
   if (TCTOPN(line, F_IO_MODE, &tidlin))
 
106
      SCTPUT("**** Error while opening table line.tbl");
 
107
 
 
108
   TCIGET(tidlin, &nbcol, &nbrow, &nsort, &allcol, &allrow);
 
109
 
 
110
   if (TCTOPN(lincat, F_I_MODE, &tidcat))
 
111
      SCTPUT("**** Error while opening line catalog");
 
112
 
 
113
   TCIGET(tidcat, &nbcolcat, &nbrowcat, &nsortcat, &allcolcat, &allrowcat);
 
114
 
 
115
 
 
116
   /* Finds column :WAVE in line catalog                    */
 
117
 
 
118
   TCCSER(tidcat, ":WAVE", &catwav);
 
119
      if (catwav == (-1)) SCTPUT("**** Column :WAVE not found");
 
120
 
 
121
 
 
122
   /* Verification loop */
 
123
   /* If the threshold parameter thres<0, the line catalog is explored to
 
124
      make sure that wavelengths are sorted by increasing value and no 
 
125
      wavelength is duplicated */
 
126
 
 
127
   if (thres[0] < 0) {
 
128
   thres[0] = 1.;
 
129
   TCERDD(tidcat, 1, catwav, &lambda_low, &null);
 
130
   for (row=2; row<=nbrowcat; row++) {
 
131
       TCERDD(tidcat, row, catwav, &lambda_upp, &null);
 
132
       if (lambda_low > lambda_upp) {
 
133
          printf("Warning: Column :WAVE of the table %s is not sorted by increasing wavelength\n",lincat);
 
134
          thres[0] = -1.;
 
135
          break;
 
136
        }
 
137
       if (lambda_low == lambda_upp) {
 
138
          printf("Error: Column :WAVE of the table %s contains duplicated wavelength : %f\n",
 
139
                  lincat, lambda_low);
 
140
          SCKWRD("OUTPUTD",&lambda_low,1,1,&kunit);
 
141
          thres[0] = -2.;
 
142
          break;
 
143
        }
 
144
       lambda_low = lambda_upp;
 
145
     }
 
146
   status = SCKWRR("INPUTR", &thres[0], 1, 1, &kunit);
 
147
   SCSEPI();
 
148
   exit(0);
 
149
 }
 
150
 
 
151
   /* Locates columns :IDENT, :WAVEC, :RESIDUAL in table LINE */
 
152
 
 
153
   TCCSER(tidlin, ":IDENT", &linid);
 
154
      if (linid == (-1)) SCTPUT("**** Column :IDENT not found");
 
155
 
 
156
   TCCSER(tidlin, ":WAVEC", &linwav);
 
157
      if (linwav == (-1)) SCTPUT("**** Column :WAVEC not found");
 
158
 
 
159
   TCCSER(tidlin, ":Y", &liny);
 
160
      if (liny == (-1)) SCTPUT("**** Column :Y not found");
 
161
 
 
162
   TCCSER(tidlin, ":ERROR", &linerr);
 
163
   TCCSER(tidlin, ":DISTANCE", &liniso);
 
164
 
 
165
   TCCSER(tidlin, ":RESIDUAL", &lindif);
 
166
      if (lindif == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom",
 
167
                               "RESIDUAL", &lindif);
 
168
 
 
169
 
 
170
   bufsize = 11*nbrow;
 
171
 
 
172
   SCFCRE("buffer.bdf",D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,bufsize,&imno);
 
173
   SCFMAP(imno,F_X_MODE,1, bufsize, &actvals, (char **)&pntr_buf);
 
174
 
 
175
   /* Matching loop */
 
176
 
 
177
   SCTPUT ("Identifying lines...\n");
 
178
 
 
179
   rowcat = 2L;
 
180
   for (row = 1; row <= nbrow; row++) {
 
181
 
 
182
      /* Reads current wavelength to be identified and relative order number */
 
183
      TCERDD(tidlin, row, linwav, &lambda_lin, &null);
 
184
      TCERDI(tidlin, row, liny,   &ordref, &null);
 
185
 
 
186
      /* Read the error on wavelength if the column exists */
 
187
 
 
188
      if (linerr > 0) 
 
189
        TCERDD(tidlin, row, linerr,   &error, &null);
 
190
      else
 
191
        error = -1.;
 
192
 
 
193
      /* printf ("Reading line.tbl row = %d line = %f\n",row,lambda_lin); */
 
194
 
 
195
      pntr_buf[ipos(6,row,nbrow)] = 0.;
 
196
 
 
197
      if (lambda_lin != 0.) {
 
198
 
 
199
          pntr_buf[ipos(6,row,nbrow)] = 1.;
 
200
          reslow_abs = 99999.9;
 
201
          resupp_abs = 99999.9;
 
202
 
 
203
          /* Finds the residual to the nearest lower line in the same order */
 
204
          if (row > 1) {
 
205
            rowlow = row - 1;
 
206
            TCERDI(tidlin, rowlow, liny,  &order, &null);
 
207
            if (order == ordref) {
 
208
              TCERDD(tidlin, rowlow, linwav, &lambda_low, &null);
 
209
              reslow = lambda_lin - lambda_low;
 
210
              reslow_abs = (reslow < 0) ? -reslow : reslow; /* abs. value */
 
211
            }
 
212
          }
 
213
 
 
214
          /* Finds the residual to the nearest upper line in the same order */
 
215
          if (row < nbrow) {
 
216
            rowupp = row + 1;
 
217
            TCERDI(tidlin, rowupp, liny,  &order, &null);
 
218
            if (order == ordref) {
 
219
              TCERDD(tidlin, rowupp, linwav, &lambda_upp, &null);
 
220
              resupp = lambda_lin - lambda_upp;
 
221
              resupp_abs = (resupp < 0) ? -resupp : resupp; /* abs. value */
 
222
            }
 
223
          }
 
224
 
 
225
          /* Compares both residuals and selects the minimum */
 
226
          if (row == 1) reslow_abs = resupp_abs;
 
227
          if (row == nbrow) resupp_abs = reslow_abs;
 
228
          /* printf("Table line: low,high %f %f\n",reslow_abs, resupp_abs); */
 
229
 
 
230
          if (reslow_abs < resupp_abs)
 
231
             pntr_buf[ipos(1,row,nbrow)] = reslow_abs;
 
232
          else
 
233
             pntr_buf[ipos(1,row,nbrow)] = resupp_abs;
 
234
 
 
235
          
 
236
          /* Reads the line catalog at the first entry point and initializes variables */
 
237
          /* printf ("Reading catalog thar.tbl rowcat = %d\n",rowcat); */
 
238
          if (rowcat > 1) {
 
239
          rowlow = rowcat - 1;
 
240
          TCERDD(tidcat, rowlow, catwav, &lambda_low, &null);
 
241
          rowupp = rowcat;
 
242
          TCERDD(tidcat, rowupp, catwav, &lambda_upp, &null);
 
243
          reslow = lambda_lin - lambda_low;
 
244
          reslow_abs = (reslow < 0) ? -reslow : reslow; /* abs. value */
 
245
          resupp = lambda_lin - lambda_upp;
 
246
          resupp_abs = (resupp < 0) ? -resupp : resupp; /* abs. value */
 
247
        }
 
248
 
 
249
          /* The line catalog is assumed to be ordered. Depending on which direction
 
250
          is found the minimum residual, the catalog will be searched further in this
 
251
          direction */
 
252
 
 
253
          if (reslow_abs < resupp_abs) {
 
254
              index = -1;
 
255
              resref = reslow;
 
256
              resref_abs = reslow_abs;
 
257
              lambda_ref = lambda_low;
 
258
              rowcat = rowlow;
 
259
            }
 
260
          else {
 
261
              index = 1;
 
262
              resref = resupp;
 
263
              resref_abs = resupp_abs;
 
264
              lambda_ref = lambda_upp;
 
265
              rowcat = rowupp;
 
266
            }
 
267
 
 
268
          res_abs = resref_abs;
 
269
          residual = resref;
 
270
          lambda_cat = lambda_ref;
 
271
 
 
272
          /* Explores the catalog in the direction indicated by index and finds the
 
273
          position of the minimum residual */
 
274
          while (res_abs <= resref_abs) {
 
275
             resref = residual;
 
276
             resref_abs = res_abs;
 
277
             lambda_ref = lambda_cat;
 
278
             rowcat += index;
 
279
             if (rowcat<=1 || rowcat >= (nbrowcat-1)) {
 
280
                 printf ("Error: reached limit of line catalog\n");
 
281
                 printf ("Wavelength: %f Position: %d\n",lambda_lin,rowcat);
 
282
                 break;
 
283
               }
 
284
             else {
 
285
             TCERDD(tidcat, rowcat, catwav, &lambda_cat, &null);
 
286
  /* printf ("Matching loop at rowcat = %d cat = %f\n",rowcat,lambda_cat); */
 
287
             residual = lambda_lin - lambda_cat;
 
288
             res_abs = (residual < 0) ?  -residual : residual; /* abs. value */
 
289
           }
 
290
           }
 
291
 
 
292
             rowref = rowcat - index;
 
293
             /* printf("Ref: Read wavelength %f at position %d\n",lambda_ref,rowref); */
 
294
             /* rowref is now the position of the minimum residual */
 
295
             pntr_buf[ipos(2,row,nbrow)] = resref_abs;
 
296
             pntr_buf[ipos(3,row,nbrow)] = lambda_ref;
 
297
             pntr_buf[ipos(4,row,nbrow)] = lambda_lin;
 
298
 
 
299
          if (rowref > 1) {
 
300
            rowlow = rowref - 1;
 
301
            TCERDD(tidcat, rowlow, catwav,  &lambda_low, &null);
 
302
             /* printf("Low: Read wavelength %f at position %d\n",lambda_low,rowlow); */
 
303
             reslow = lambda_ref - lambda_low;
 
304
             reslow_abs = (reslow < 0) ? -reslow : reslow; /* abs. value */
 
305
            }
 
306
 
 
307
          if (rowref < nbrowcat && rowref > 1) {
 
308
            rowupp = rowref + 1;
 
309
            TCERDD(tidcat, rowupp, catwav,  &lambda_upp, &null);
 
310
            /* printf("Upp: Read wavelength %f at position %d\n",lambda_upp,rowupp); */
 
311
            resupp = lambda_ref - lambda_upp;
 
312
            resupp_abs = (resupp < 0) ? -resupp : resupp; /* abs. value */
 
313
            }
 
314
 
 
315
          if (rowref == 1) reslow_abs = resupp_abs;
 
316
          if (rowref == nbrowcat) resupp_abs = reslow_abs;
 
317
          /* printf("Catalog: low,high %f %f\n",reslow_abs, resupp_abs); */
 
318
 
 
319
          if (reslow_abs < resupp_abs)
 
320
             pntr_buf[ipos(5,row,nbrow)] = reslow_abs;
 
321
          else
 
322
             pntr_buf[ipos(5,row,nbrow)] = resupp_abs;
 
323
 
 
324
 
 
325
          if (pntr_buf[ipos(1,row,nbrow)] < pntr_buf[ipos(5,row,nbrow)])
 
326
             reslow = pntr_buf[ipos(1,row,nbrow)];
 
327
          else
 
328
             reslow = pntr_buf[ipos(5,row,nbrow)];
 
329
          
 
330
          reslow = reslow * thres[0];
 
331
          if (error < 0.) error = reslow;
 
332
          finres = pntr_buf[ipos(2,row,nbrow)];
 
333
 
 
334
          pntr_buf[ipos(7,row,nbrow)] = lambda_lin;
 
335
          pntr_buf[ipos(8,row,nbrow)] = reslow;
 
336
          pntr_buf[ipos(9,row,nbrow)] = error;
 
337
 
 
338
          pntr_buf[ipos(10,row,nbrow)] = lambda_ref;
 
339
          pntr_buf[ipos(11,row,nbrow)] = resref;
 
340
 
 
341
 
 
342
      }                         /* End of if lambda_lin != 0.  */
 
343
 
 
344
/*    printf("Row %d Order %d - %f  %f  %f  %f  %f\n",row,ordref, */
 
345
/*    pntr_buf[ipos(1,row,nbrow)],pntr_buf[ipos(2,row,nbrow)],  */
 
346
/*    pntr_buf[ipos(3,row,nbrow)],pntr_buf[ipos(4,row,nbrow)],  */
 
347
/*    pntr_buf[ipos(5,row,nbrow)]);  */
 
348
 
 
349
   }                            /* End of for row = */
 
350
 
 
351
 
 
352
   /* Determines the standard error on a fraction of the most isolated lines           TCEWRD(tidlin, row, liniso, &reslow); */
 
353
 
 
354
   number = (int) ((float) nbrow / 10); 
 
355
   if (number < 10 ) number = 10;
 
356
   stdev(number,&error_min);
 
357
   sprintf(Msg,"Std Dev on %d lines: %g wav. units \n",number,error_min);
 
358
   SCTPUT (Msg);
 
359
 
 
360
   open_acc();
 
361
 
 
362
   for (row = 1; row <= nbrow; row++) {
 
363
     lambda_lin = pntr_buf[ipos(7,row,nbrow)];
 
364
 
 
365
     if (lambda_lin != 0.) {
 
366
       finres = pntr_buf[ipos(2,row,nbrow)];
 
367
       reslow = pntr_buf[ipos(8,row,nbrow)];
 
368
 
 
369
       error  = pntr_buf[ipos(9,row,nbrow)];
 
370
 
 
371
       accumulate(log(finres/error),log(reslow/error));
 
372
     }}
 
373
 
 
374
   findmax(fraction); 
 
375
   close_acc();
 
376
 
 
377
   for (row = 1; row <= nbrow; row++) {
 
378
 
 
379
   lambda_lin = pntr_buf[ipos(7,row,nbrow)];
 
380
 
 
381
   if (lambda_lin != 0.) {
 
382
 
 
383
    finres = pntr_buf[ipos(2,row,nbrow)];
 
384
    reslow = pntr_buf[ipos(8,row,nbrow)];
 
385
    error  = pntr_buf[ipos(9,row,nbrow)]*exp(maxval);
 
386
    /* error = 2.5*error_min; */
 
387
 
 
388
    lambda_ref = pntr_buf[ipos(10,row,nbrow)];
 
389
    resref     = pntr_buf[ipos(11,row,nbrow)];
 
390
 
 
391
   /* printf("Comparing line %f with catalog %f.\n ",lambda_ref,lambda_lin);
 
392
   printf("Residual is %f and maximum residual is %f\n",finres,reslow); */
 
393
 
 
394
   TCEWRD(tidlin, row, liniso, &reslow);
 
395
   TCEWRD(tidlin, row, lindif, &resref);
 
396
   if (finres < error && error <= reslow) { 
 
397
          TCEWRD(tidlin, row, linid,  &lambda_ref);
 
398
          /* printf("**** Identified line: %f  %f\n",lambda_ref,resref); */
 
399
        }
 
400
 
 
401
  }} /* End of if lambda_lin and for (row = */
 
402
 
 
403
   TCTCLO(tidlin);
 
404
   TCTCLO(tidcat);
 
405
   SCSEPI();
 
406
return 0;
 
407
}
 
408
 
 
409
 
 
410
 
 
411
void open_acc()
 
412
 
 
413
{
 
414
 
 
415
char  ident[81], cunit[33];
 
416
int   i;
 
417
 
 
418
strcpy(ident,"Error proportionality coefficient accumulator");
 
419
strcpy(cunit,"Frequency                       ");
 
420
 
 
421
SCIPUT(name,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,1,&npix,
 
422
       &start,&step,ident,cunit, (char **)&pntr, &imno);
 
423
 
 
424
stend = start + (npix-1)*step;
 
425
 
 
426
for (i=0; i<npix; i++) pntr[i] = 0.;
 
427
 
 
428
}
 
429
 
 
430
 
 
431
void accumulate(min,max)
 
432
 
 
433
double min,max;
 
434
 
 
435
{
 
436
 
 
437
int    pixmin, pixmax, i;
 
438
double vmin, vmax;
 
439
 
 
440
if (max < min) return;
 
441
 
 
442
if (min < start) vmin = start;
 
443
else             vmin = min;
 
444
 
 
445
if (max > stend)   vmax = stend;
 
446
else             vmax = max;
 
447
 
 
448
pixmin = (vmin-start)/step;
 
449
pixmax = (vmax-start)/step;
 
450
 
 
451
for (i=pixmin; i<=pixmax; i++) pntr[i] += 1;
 
452
 
 
453
}
 
454
 
 
455
 
 
456
void close_acc()
 
457
{SCFCLO(imno);}
 
458
 
 
459
 
 
460
void findmax(fraction)
 
461
 
 
462
float fraction;
 
463
 
 
464
{
 
465
 
 
466
int    pixmax, i, pixopt;
 
467
double vmax, vopt;
 
468
 
 
469
vmax   = pntr[0];
 
470
pixmax = 0;
 
471
 
 
472
for (i=1; i<npix; i++) {
 
473
    if (pntr[i] > vmax) {
 
474
        vmax = pntr[i];
 
475
        pixmax = i;
 
476
      }
 
477
  }
 
478
 
 
479
maxval = start + (pixmax-1)*step;
 
480
/* printf("Found optimum at %f . Value: %f\n",maxval,pntr[pixmax]); */
 
481
 
 
482
vopt = fraction*pntr[pixmax];
 
483
pixopt = pixmax;
 
484
 
 
485
/* printf("Fraction: %f , Ref level: %f\n",fraction,vopt); */
 
486
 
 
487
while (vopt < pntr[pixopt] && pixopt < npix)  {
 
488
      pixopt++;
 
489
      /* printf("Current pos. %d Level: %f\n",pixopt,pntr[pixopt]); */
 
490
    }
 
491
 
 
492
maxval = start + (pixopt-1)*step;
 
493
/* printf("Selection level at %f . Value: %f\n",maxval,pntr[pixopt]); */
 
494
 
 
495
}
 
496
 
 
497
void stdev(number, error)
 
498
 
 
499
int    number;
 
500
double *error;
 
501
 
 
502
{
 
503
 
 
504
int row, posval, nb, flag, i;
 
505
double maxval, val, res, lambda;
 
506
 
 
507
for (nb=1; nb<=number; nb++) {
 
508
 
 
509
maxval = -1.; /* Sort positive values */
 
510
posval = 0;
 
511
 
 
512
for (row=1;row<=nbrow;row++) {
 
513
 
 
514
    flag = 0;
 
515
 
 
516
    lambda = pntr_buf[ipos(7,row,nbrow)];
 
517
    if (lambda < 1.) flag = 1;
 
518
 
 
519
    for(i=1; i<nb; i++) {
 
520
        if (row == tmp_pos[i])  {
 
521
           flag = 1;
 
522
         }}
 
523
 
 
524
    if (flag == 0) {
 
525
      val = pntr_buf[ipos(8,row,nbrow)];
 
526
      res = pntr_buf[ipos(2,row,nbrow)];
 
527
      if (val > res && val > maxval){
 
528
        maxval = val; posval = row;
 
529
      }}
 
530
  }
 
531
 
 
532
tmp_pos[nb] = posval;
 
533
tmp_iso[nb] = maxval;
 
534
tmp_res[nb] = pntr_buf[ipos(2,posval,nbrow)];
 
535
}
 
536
 
 
537
/* Now sorting tmp array by decreasing residual value */ 
 
538
 
 
539
{
 
540
int flag, tmpp;
 
541
double tmpr, tmpi;
 
542
 
 
543
flag = 1;
 
544
 
 
545
while (flag != 0) {
 
546
  flag =0;
 
547
  for(row=2; row<=number; row++) {
 
548
      if (tmp_res[row-1] < tmp_res[row]) {
 
549
        flag += 1;
 
550
        tmpr = tmp_res[row-1], tmpi = tmp_iso[row-1], tmpp = tmp_pos[row-1];
 
551
        tmp_res[row-1] = tmp_res[row];
 
552
        tmp_iso[row-1] = tmp_iso[row];
 
553
        tmp_pos[row-1] = tmp_pos[row];
 
554
        tmp_res[row] = tmpr, tmp_iso[row] = tmpi, tmp_pos[row] = tmpp;
 
555
      }
 
556
    }
 
557
}}
 
558
 
 
559
/*
 
560
std = 0.;
 
561
for (row=1; row<=number; row++) {
 
562
  res = tmp_res[row]; 
 
563
  std += res*res;  
 
564
}
 
565
*error = sqrt(std/(nb-1));
 
566
*/
 
567
 
 
568
*error = tmp_res[number/2];
 
569
 
 
570
/* Writing the residuals in a small auxiliary ascii file */
 
571
 
 
572
{
 
573
FILE   *fp;
 
574
fp = fopen("dat.dat","w");
 
575
for (row=1; row<=number; row++) 
 
576
    fprintf(fp,"%f %f\n",tmp_res[row],tmp_iso[row]);
 
577
fclose(fp);
 
578
}
 
579
 
 
580
}