← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/long/src/lnplot.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/long/src/lnplot.c
 
1
/*===========================================================================
 
2
  Copyright (C) 1993-2009 European Southern Observatory (ESO)
 
3
 
 
4
  This program is free software; you can redistribute it and/or 
 
5
  modify it under the terms of the GNU General Public License as 
 
6
  published by the Free Software Foundation; either version 2 of 
 
7
  the License, or (at your option) any later version.
 
8
 
 
9
  This program is distributed in the hope that it will be useful,
 
10
  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
  GNU General Public License for more details.
 
13
 
 
14
  You should have received a copy of the GNU General Public 
 
15
  License along with this program; if not, write to the Free 
 
16
  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
17
  MA 02139, USA.
 
18
 
 
19
  Correspondence concerning ESO-MIDAS should be addressed as follows:
 
20
        Internet e-mail: midas@eso.org
 
21
        Postal address: European Southern Observatory
 
22
                        Data Management Division 
 
23
                        Karl-Schwarzschild-Strasse 2
 
24
                        D 85748 Garching bei Muenchen 
 
25
                        GERMANY
 
26
===========================================================================*/
 
27
 
 
28
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
 
29
/* .IDENT       lnplot.c                                   */
 
30
/* .AUTHORS     Pascal Ballester (ESO/Garching)            */
 
31
/*              Cristian Levin   (ESO/La Silla)            */
 
32
/* .KEYWORDS    Spectroscopy, Long-Slit                    */
 
33
/* .PURPOSE                                                */
 
34
/* .VERSION     1.0  Package Creation  17-MAR-1993         */
 
35
 
 
36
/*  090728      last modif                                 */
 
37
/* ------------------------------------------------------- */
 
38
 
 
39
#include <stdio.h>
 
40
#include <math.h>
 
41
 
 
42
#include <midas_def.h>
 
43
 
 
44
#include <gl_defs.h>
 
45
#include <main_defs.h>
 
46
#include <spec_defs.h>
 
47
#include <agl_defs.h>
 
48
 
 
49
#define MINWAVE         -99999.0
 
50
#define MAXWAVE         99999.0
 
51
#define MINY            -99999.0
 
52
#define MAXY            99999.0
 
53
#define DBORDER         10.0
 
54
#define NPOINTS         500
 
55
#define EPS             0.001
 
56
 
 
57
/* Actions */
 
58
#define PLOT            0
 
59
#define DELETE          1
 
60
#define GETCUR          2
 
61
#define PRINT           3
 
62
 
 
63
#define MAXDEL          80              /* max. temp. deleted points */
 
64
 
 
65
/* Keywords */
 
66
char    Lincat[MAXLINE];
 
67
char    Wlc[21];
 
68
int     Wrang[2];
 
69
float   Imin;
 
70
int     Ystart;
 
71
int     Fitd;
 
72
int     PlotType;       /* type of current plot */
 
73
int     PlotAction;     /* Action to do */
 
74
 
 
75
/* NULL values */
 
76
int     Inull;
 
77
float   Rnull;
 
78
double  Dnull;
 
79
 
 
80
/* structures for allocating column values of 'LINTAB' */
 
81
float *X, *Peak, *Ident, *Wave, *Wavec, *Delta, *Deltac;
 
82
int *Row;       /* row number associated to each value */
 
83
 
 
84
int     NumLine;                /* Number of rows read from 'LINTAB' */
 
85
float   Wlimits[4];             /* window limits */
 
86
 
 
87
int     DelList[MAXDEL];        /* List of temporarily deleted points */
 
88
int     NumDel = 0;             /* DelList[] cardinality */
 
89
 
 
90
char    CurrFile[81];           /* Global variable used for plot spectra */
 
91
char    DevErase[80], DevNoErase[80];
 
92
 
 
93
char    PlotFile[MAXLINE], Lintab[MAXLINE], Coerbr[MAXLINE];
 
94
int     Ycoerbr;        /* row nearest to :ystart in 'COERBR' */
 
95
double  Rms, Dispersion;
 
96
double  Coef[20];
 
97
 
 
98
LCTAB * Lc = NULL;
 
99
 
 
100
int     InitGraphic = FALSE;
 
101
 
 
102
struct column_ids {
 
103
        int x, y, peak, ident;
 
104
        int wave, wavec, delta, deltac;
 
105
        int erased;
 
106
} Col;
 
107
 
 
108
extern void free_fvector(), free_dvector(), free_ivector();
 
109
extern void lfit(), free_catalog_table();
 
110
 
 
111
extern int file_exists(), graphwin_exists(), read_catalog_table();
 
112
 
 
113
void compute_calib_values();
 
114
void read_coefs_ystart();
 
115
void read_rebin_parameters();
 
116
void read_column_ids();
 
117
void free_data();
 
118
void init_arrays_data();
 
119
void read_line_table();
 
120
void read_parameters();
 
121
void init_midas();
 
122
void end_midas();
 
123
void get_agldev();
 
124
void end_graphic();
 
125
void init_graphic();
 
126
void read_image();
 
127
void init_viewport();
 
128
void read_limits();
 
129
void save_limits();
 
130
 
 
131
void getcur_wave();
 
132
void plot_wave();
 
133
void plot_ident();
 
134
void plot_splabel();
 
135
void init_ident();
 
136
void init_splabel();
 
137
void plot_curve();
 
138
void plot_delete();
 
139
void save_limits();
 
140
void read_limits();
 
141
 
 
142
void del_ident();
 
143
void del_point();
 
144
void undel_point();
 
145
void dpoly();
 
146
 
 
147
float *fvector();
 
148
 
 
149
int read_lincat_table(), point_deleted();
 
150
 
 
151
double *dvector();
 
152
double eval_dpoly();
 
153
 
 
154
int main()
 
155
{
 
156
    init_midas();
 
157
 
 
158
    read_parameters();
 
159
    read_lincat_table();
 
160
    read_rebin_parameters();
 
161
    read_coefs_ystart();
 
162
    read_line_table();
 
163
 
 
164
    if ( PlotAction == PLOT ) {
 
165
        init_graphic(DEV_ERASE);
 
166
        AG_MOPN( PlotFile );
 
167
        AG_SSET("FONT=1");     /* nice font */
 
168
        switch( PlotType ) {
 
169
            case PLOT_IDENT:    plot_ident(); break;
 
170
            case PLOT_WAVE:     plot_wave(); break;
 
171
            case PLOT_SPLABEL:  plot_splabel(); break;
 
172
        }
 
173
        AG_MCLS();
 
174
    }
 
175
    else if ( PlotAction == DELETE ) {
 
176
        init_graphic(DEV_NO_ERASE);
 
177
        AG_SSET("FONT=1");
 
178
        switch( PlotType ) {
 
179
            case PLOT_IDENT:    
 
180
            case PLOT_WAVE:     
 
181
            case PLOT_SPLABEL:  init_viewport(); break;
 
182
        }
 
183
    }
 
184
    else if ( PlotAction == GETCUR ) {
 
185
        init_graphic(DEV_NO_ERASE);
 
186
        AG_SSET("FONT=1");
 
187
        init_viewport();
 
188
        getcur_wave();
 
189
    }
 
190
 
 
191
    if ( PlotAction == DELETE && PlotType != PLOT_NONE )
 
192
        plot_delete();
 
193
 
 
194
    free_data();
 
195
    end_graphic();
 
196
    end_midas();
 
197
return 0;
 
198
}
 
199
 
 
200
void plot_delete()
 
201
{
 
202
    float cpx = 0.0, cpy = 0.0;
 
203
    int key, valpix;
 
204
    int i = 0, min = 0;
 
205
    float min_deltax, min_deltay;
 
206
    float x[2], y[2], xi, yi;
 
207
    char str[80];
 
208
 
 
209
    while (1) {
 
210
        min_deltax =  min_deltay = MAXWAVE;
 
211
 
 
212
        AG_VLOC( &cpx, &cpy, &key, &valpix );
 
213
        if ( key == MID_BUT )
 
214
            break;
 
215
 
 
216
        AG_SSET( RED_COLOR );
 
217
        switch( PlotType ) {
 
218
           case PLOT_IDENT:
 
219
                for ( i = 0; i < NumLine; i++ ) {
 
220
                    xi = Wavec[i]; 
 
221
                    yi = Deltac[i] + Wavec[i] - Ident[i];
 
222
                    if ( Wave[i] != Rnull &&
 
223
                         fabs((double) cpx - xi) < min_deltax &&
 
224
                         fabs((double) cpy - yi) < min_deltay ) {
 
225
                        min_deltax = fabs((double) cpx - xi);
 
226
                        min_deltay = fabs((double) cpx - yi);
 
227
                        min = i;
 
228
                    }
 
229
                }
 
230
                x[0] = Wavec[min]; 
 
231
                y[0] = Deltac[min] + Wavec[min] - Ident[min];
 
232
                del_ident(min);
 
233
                sprintf(str, "point deleted : %10.3f %10.3f", x[0], y[0]);
 
234
                SCTPUT(str);
 
235
                AG_GPLM( x, y, 1, BOX_MARKER );
 
236
                break;
 
237
           case PLOT_WAVE:
 
238
                for ( i = 0; i < NumLine; i++ ) {
 
239
                    xi = Wave[i]; 
 
240
                    yi = Delta[i];
 
241
                    if ( Wave[i] != Rnull &&
 
242
                         fabs((double) cpx - xi) < min_deltax &&
 
243
                         fabs((double) cpy - yi) < min_deltay ) {
 
244
                        min_deltax = fabs((double) cpx - xi);
 
245
                        min_deltay = fabs((double) cpx - yi);
 
246
                        min = i;
 
247
                    }
 
248
                }
 
249
                x[0] = Wave[min]; 
 
250
                y[0] = Delta[min];
 
251
                if ( point_deleted(min) ) {
 
252
                    AG_SSET( BLUE_COLOR );
 
253
                    undel_point(min);
 
254
                    sprintf(str, "point added   : %10.3f %10.3f", x[0], y[0]);
 
255
                }
 
256
                else {
 
257
                    del_point(min);
 
258
                    sprintf(str, "point deleted : %10.3f %10.3f", x[0], y[0]);
 
259
                }
 
260
                SCTPUT(str);
 
261
                AG_GPLM( x, y, 1, X_MARKER );
 
262
                break;
 
263
           case PLOT_SPLABEL:
 
264
                for ( i = 0; i < NumLine; i++ ) {
 
265
                    if ( Wave[i] != Rnull &&
 
266
                         fabs((double) cpx - X[i]) < min_deltax ) {
 
267
                        min_deltax = fabs((double) cpx - X[i]);
 
268
                        min = i;
 
269
                    }
 
270
                }
 
271
                x[0] = x[1] = X[min];
 
272
                AG_RGET( "WNDL", Wlimits );
 
273
                y[0] = Wlimits[YMIN];
 
274
                y[1] = Peak[min];
 
275
                if ( point_deleted(min) ) {
 
276
                    AG_SSET( ((Ident[min] != Rnull) ? GREEN_COLOR:BLUE_COLOR) );
 
277
                    undel_point(min);
 
278
                    sprintf(str, "point added   : %10.3f", x[0]);
 
279
                }
 
280
                else {
 
281
                    del_point(min);
 
282
                    sprintf(str, "point deleted : %10.3f", x[0]);
 
283
                }
 
284
                SCTPUT(str);
 
285
                AG_GPLL( x, y, 2 );
 
286
                break;
 
287
        }
 
288
        AG_VUPD();
 
289
        AG_SSET( DEF_COLOR );
 
290
    }
 
291
    NumDel = 0; /* no temporarily deleted points */
 
292
}
 
293
 
 
294
void del_ident(n)
 
295
int n;
 
296
{
 
297
    int id;
 
298
 
 
299
    TCTOPN( Lintab, F_IO_MODE, &id );
 
300
    TCEWRR( id, Row[n], Col.ident, &Rnull );
 
301
    TCTCLO(id);
 
302
}
 
303
 
 
304
void del_point(n)
 
305
int n;
 
306
{
 
307
    int id;
 
308
    char del_str[2];
 
309
 
 
310
    sprintf(del_str, "%c", VAL_ERASED);
 
311
    DelList[NumDel++] = n;
 
312
 
 
313
    TCTOPN( Lintab, F_IO_MODE, &id );
 
314
    TCEWRC( id, Row[n], Col.erased, del_str );
 
315
    TCTCLO(id);
 
316
}
 
317
 
 
318
void undel_point(n)
 
319
int n;
 
320
{
 
321
    int  select = TRUE;
 
322
    int i, id, min = -1;
 
323
    char undel_str[2];
 
324
    float min_delta = MAXWAVE;
 
325
 
 
326
    sprintf(undel_str, "%c", VAL_OK);
 
327
 
 
328
    TCTOPN( Lintab, F_IO_MODE, &id );
 
329
    for ( i = 0; i < NumDel; i++ )
 
330
        if ( n == DelList[i] ) {
 
331
            DelList[i] = DelList[NumDel-1];
 
332
            TCEWRC( id, Row[n], Col.erased, undel_str );
 
333
            NumDel--;
 
334
            break;
 
335
        }
 
336
    TCTCLO(id);
 
337
 
 
338
    TCTOPN( Lincat, F_IO_MODE, &id );
 
339
    for ( i = 0; i < Lc->nrows; i++ )
 
340
        if ( fabs(Lc->wave[i] - Wave[n]) < min_delta && ! Lc->sel[i] ) {
 
341
            min = i;
 
342
            min_delta = fabs(Lc->wave[i] - Wave[n]);
 
343
        }
 
344
    if ( min != -1 )
 
345
        TCSPUT(id, Lc->row[min], &select);
 
346
    TCTCLO(id);
 
347
}
 
348
 
 
349
int point_deleted(n)
 
350
int n;
 
351
{
 
352
    int i;
 
353
 
 
354
    for ( i = 0; i < NumDel; i++ )
 
355
        if ( n == DelList[i] )
 
356
            return( TRUE );
 
357
    return( FALSE );
 
358
}
 
359
 
 
360
void plot_wave()
 
361
{
 
362
    int i, n1 = 0, n2 = 0;
 
363
    int ncoef;
 
364
    int ndel = 0;
 
365
    float *xdel, *ydel;
 
366
    float xmin = MAXWAVE, xmax = MINWAVE;
 
367
    float ymin = MAXY, ymax = MINY;
 
368
    float *x1, *x2, *y1, *y2;
 
369
    char options[MAXOPTS];
 
370
    char text[80];
 
371
 
 
372
    /*** get the values to plot ***/
 
373
 
 
374
    x1 = fvector( 0, NumLine-1 );
 
375
    y1 = fvector( 0, NumLine-1 );
 
376
    x2 = fvector( 0, NumLine-1 );
 
377
    y2 = fvector( 0, NumLine-1 );
 
378
    xdel = fvector( 0, NumLine-1 );
 
379
    ydel = fvector( 0, NumLine-1 );
 
380
 
 
381
    for ( i = 0; i < NumLine; i++ ) {
 
382
        if ( Ident[i] != Rnull ) {      /* get :IDENT points */
 
383
            x1[n1] = Wavec[i];
 
384
            y1[n1] = Deltac[i] + Wavec[i] - Ident[i];
 
385
            if ( x1[n1] < xmin ) xmin = x1[n1];
 
386
            if ( x1[n1] > xmax ) xmax = x1[n1];
 
387
            if ( y1[n1] < ymin ) ymin = y1[n1];
 
388
            if ( y1[n1] > ymax ) ymax = y1[n1];
 
389
            n1++;
 
390
        }
 
391
        if ( point_deleted(i) ) {       /* get deleted points & continue */
 
392
            xdel[ndel] = Wave[i];
 
393
            ydel[ndel] = Delta[i];
 
394
            if ( xdel[ndel] < xmin ) xmin = xdel[ndel];
 
395
            if ( xdel[ndel] > xmax ) xmax = xdel[ndel];
 
396
            if ( ydel[ndel] < ymin ) ymin = ydel[ndel];
 
397
            if ( ydel[ndel] > ymax ) ymax = ydel[ndel];
 
398
            ndel++;
 
399
            continue;
 
400
        }
 
401
        if ( Wave[i] != Rnull ) {       /* get :WAVE points not deleted */
 
402
            x2[n2] = Wave[i];
 
403
            y2[n2] = Delta[i];
 
404
            if ( x2[n2] < xmin ) xmin = x2[n2];
 
405
            if ( x2[n2] > xmax ) xmax = x2[n2];
 
406
            if ( y2[n2] < ymin ) ymin = y2[n2];
 
407
            if ( y2[n2] > ymax ) ymax = y2[n2];
 
408
            n2++;
 
409
        }
 
410
    }
 
411
 
 
412
    /*** plot... ***/
 
413
 
 
414
    /* let's give a border for a nicer plot */
 
415
    xmin -= fabs( (double) (xmax - xmin) / DBORDER );
 
416
    xmax += fabs( (double) (xmax - xmin) / DBORDER );
 
417
    ymin -= fabs( (double) (ymax - ymin) / DBORDER );
 
418
    ymax += fabs( (double) (ymax - ymin) / DBORDER );
 
419
 
 
420
    strcpy( options, "LABY=Delta(Wave);LABX=Wavelength" );
 
421
    AG_AXES( xmin, xmax, ymin, ymax, options );
 
422
    
 
423
    /* Plot the :IDENT related points */
 
424
    if (n1 > 0) {
 
425
        AG_SSET( GREEN_COLOR );
 
426
        AG_GPLM( x1, y1, n1, BOX_MARKER );
 
427
        AG_VUPD();
 
428
        AG_SSET( DEF_COLOR );
 
429
      }
 
430
 
 
431
    /* Plot the :WAVEC related points */
 
432
    AG_SSET( BLUE_COLOR );
 
433
    AG_GPLM( x2, y2, n2, X_MARKER );
 
434
    AG_VUPD();
 
435
    AG_SSET( DEF_COLOR );
 
436
 
 
437
    /* Plot the deleted points */
 
438
    if ( ndel > 0 ) {
 
439
        AG_SSET( RED_COLOR );
 
440
        AG_GPLM( xdel, ydel, ndel, X_MARKER );
 
441
        AG_VUPD();
 
442
        AG_SSET( DEF_COLOR );
 
443
    }
 
444
 
 
445
    /* Plot horizontal line at y = 0 */
 
446
    AG_SSET( DASH_LINE );
 
447
    x1[0] = xmin; x1[1] = xmax;
 
448
    y1[0] = y1[1] = 0.0;
 
449
    AG_GPLL( x1, y1, 2 );
 
450
    AG_VUPD();
 
451
    AG_SSET( SOLID_LINE );
 
452
 
 
453
    ncoef = (n2 > Fitd) ? Fitd+1 : n2; /* coeffs. to obtain */
 
454
    plot_curve( x2, y2, n2, xmin, xmax, ncoef );
 
455
 
 
456
    AG_RGET( "WNDL", Wlimits );
 
457
    AG_SSET( LR_TEXT );
 
458
    sprintf( text, "RMS=%.3f   DISP=%.3f", Rms, Dispersion );
 
459
    AG_GTXT( Wlimits[XMAX], Wlimits[YMAX], text, 18 );
 
460
    AG_VUPD();
 
461
 
 
462
    free_fvector( x1, 0, NumLine-1 );
 
463
    free_fvector( y1, 0, NumLine-1 );
 
464
    free_fvector( x2, 0, NumLine-1 );
 
465
    free_fvector( y2, 0, NumLine-1 );
 
466
    free_fvector( xdel, 0, NumLine-1 );
 
467
    free_fvector( ydel, 0, NumLine-1 );
 
468
 
 
469
    save_limits( xmin, xmax, ymin, ymax ); /* for init_viewport() */
 
470
}
 
471
 
 
472
void plot_curve( x, y, n, xmin, xmax, ncoef )
 
473
float *x, *y;
 
474
float xmin, xmax;
 
475
int n, ncoef;
 
476
{
 
477
    float *px, *py, dx, step;
 
478
    double *a; /* coefficients for fitting */
 
479
    double *p, *xfit, *yfit;
 
480
    int np, intervals;
 
481
    int i;
 
482
 
 
483
    a = dvector( 1, ncoef );
 
484
    p = dvector( 1, ncoef );
 
485
    xfit = dvector( 1, n );
 
486
    yfit = dvector( 1, n );
 
487
 
 
488
    intervals = NPOINTS + 2;
 
489
    step = (xmax - xmin) / NPOINTS;
 
490
    px = fvector( 0, intervals-1 );
 
491
    py = fvector( 0, intervals-1 );
 
492
 
 
493
    for ( i = 1; i <= n; i++ ) {
 
494
        xfit[i] = x[i-1];
 
495
        yfit[i] = y[i-1];
 
496
    }
 
497
    lfit( xfit, yfit, n, a, ncoef, dpoly ); /* fitting */
 
498
   
 
499
    np = 0;
 
500
    for ( dx = xmin; dx <= xmax; dx += step ) {
 
501
        px[np] = dx;
 
502
        py[np] = eval_dpoly(dx, a, ncoef);
 
503
        np++;
 
504
    }
 
505
    AG_GPLL( px, py, np );
 
506
    AG_VUPD();
 
507
 
 
508
    free_fvector( px, 0, intervals-1 );
 
509
    free_fvector( py, 0, intervals-1 );
 
510
    free_dvector( a, 1, Fitd+1 );
 
511
    free_dvector( p, 1, Fitd+1 );
 
512
    free_dvector( xfit, 1, n );
 
513
    free_dvector( yfit, 1, n );
 
514
}
 
515
 
 
516
void plot_splabel()
 
517
{
 
518
    char text[80];
 
519
    int i;
 
520
    float xdel[2], ydel[2];
 
521
 
 
522
    if ( !file_exists( Wlc, ".bdf" ) ) {
 
523
        SCTPUT( "*** Calibration image doesn't exist ***" );
 
524
        end_midas();
 
525
    }
 
526
    read_image( Ystart, Wlc );
 
527
 
 
528
    AG_RGET( "WNDL", Wlimits );
 
529
    ydel[0] = Wlimits[YMIN];
 
530
 
 
531
    AG_SSET( DU_TEXT );
 
532
    AG_SSET( DEF_CHAR );
 
533
    AG_SSET( HIGH_FONT );
 
534
    AG_SSET( BLUE_COLOR );
 
535
    for ( i = 0; i < NumLine; i++ ) {
 
536
        if ( Wave[i] != Rnull ) {
 
537
            sprintf( text, "%.1f", Wave[i] );
 
538
            AG_GTXT( X[i], Wlimits[YMAX], text, 17 );
 
539
            AG_VUPD();
 
540
        }
 
541
        if ( point_deleted(i) || Ident[i] != Rnull ) {
 
542
            AG_SSET( (point_deleted(i) ? RED_COLOR : GREEN_COLOR) );
 
543
            xdel[0] = xdel[1] = X[i];
 
544
            ydel[1] = Peak[i];
 
545
            AG_GPLL( xdel, ydel, 2 );
 
546
            AG_VUPD();
 
547
            AG_SSET( BLUE_COLOR );
 
548
        }
 
549
    }
 
550
    AG_SSET( LR_TEXT );
 
551
    AG_SSET( DEF_FONT );
 
552
    AG_SSET( DEF_COLOR );
 
553
}
 
554
 
 
555
void init_ident()
 
556
{
 
557
    int i;
 
558
    float xmin = MAXWAVE, xmax = MINWAVE;
 
559
    float ymin = MAXY, ymax = MINY;
 
560
    float x, y;
 
561
    char options[MAXOPTS];
 
562
 
 
563
    /*** get max. and min. values to plot ***/
 
564
 
 
565
    for ( i = 0; i < NumLine; i++ )
 
566
        if ( Ident[i] != Rnull ) {
 
567
            x = Wavec[i];
 
568
            y = Deltac[i] + Wavec[i] - Ident[i];
 
569
            if ( x < xmin ) xmin = x;
 
570
            if ( x > xmax ) xmax = x;
 
571
            if ( y < ymin ) ymin = y;
 
572
            if ( y > ymax ) ymax = y;
 
573
        }
 
574
 
 
575
    /* let's give a border for a nicer plot */
 
576
    xmin -= fabs( (double) (xmax - xmin) / DBORDER );
 
577
    xmax += fabs( (double) (xmax - xmin) / DBORDER );
 
578
    ymin -= fabs( (double) (ymax - ymin) / DBORDER );
 
579
    ymax += fabs( (double) (ymax - ymin) / DBORDER );
 
580
 
 
581
    strcpy( options, "LABY=Delta(Ident);LABX=Wavelength" );
 
582
    AG_AXES( xmin, xmax, ymin, ymax, options );
 
583
    AG_VUPD();
 
584
}
 
585
 
 
586
void plot_ident()
 
587
{
 
588
    int i, n = 0;
 
589
    float xmin = MAXWAVE, xmax = MINWAVE;
 
590
    float ymin = MAXY, ymax = MINY;
 
591
    float *x, *y;
 
592
    char options[MAXOPTS];
 
593
 
 
594
    /*** get the values to plot ***/
 
595
 
 
596
    x = fvector( 0, NumLine-1 );
 
597
    y = fvector( 0, NumLine-1 );
 
598
    for ( i = 0; i < NumLine; i++ )
 
599
        if ( Ident[i] != Rnull ) {
 
600
            x[n] = Wavec[i];
 
601
            y[n] = Deltac[i] + Wavec[i] - Ident[i];
 
602
            if ( x[n] < xmin ) xmin = x[n];
 
603
            if ( x[n] > xmax ) xmax = x[n];
 
604
            if ( y[n] < ymin ) ymin = y[n];
 
605
            if ( y[n] > ymax ) ymax = y[n];
 
606
            n++;
 
607
        }
 
608
 
 
609
    /*** plot... ***/
 
610
 
 
611
    /* let's give a border for a nicer plot */
 
612
    xmin -= fabs( (double) (xmax - xmin) / DBORDER );
 
613
    xmax += fabs( (double) (xmax - xmin) / DBORDER );
 
614
    ymin -= fabs( (double) (ymax - ymin) / DBORDER );
 
615
    ymax += fabs( (double) (ymax - ymin) / DBORDER );
 
616
 
 
617
    strcpy( options, "LABY=Delta(Ident);LABX=Wavelength" );
 
618
    AG_AXES( xmin, xmax, ymin, ymax, options );
 
619
    
 
620
    /* Plot the points */
 
621
    AG_SSET( BLUE_COLOR );
 
622
    AG_GPLM( x, y, n, BOX_MARKER );
 
623
    AG_VUPD();
 
624
    AG_SSET( DEF_COLOR );
 
625
 
 
626
    /* Plot horizontal line at y = 0 */
 
627
    AG_SSET( DASH_LINE );
 
628
    x[0] = xmin; x[1] = xmax;
 
629
    y[0] = y[1] = 0.0;
 
630
    AG_GPLL( x, y, 2 );
 
631
    AG_VUPD();
 
632
    AG_SSET( SOLID_LINE );
 
633
 
 
634
    free_fvector( x, 0, NumLine-1 );
 
635
    free_fvector( y, 0, NumLine-1 );
 
636
 
 
637
    save_limits( xmin, xmax, ymin, ymax ); /* for init_viewport() */
 
638
}
 
639
 
 
640
void getcur_wave()
 
641
{
 
642
    float cpx, cpy;
 
643
    float wlimits[4];   /* window limits */
 
644
    int key, valpix;
 
645
    char msg[30];
 
646
    double wave;
 
647
 
 
648
    AG_RGET( "WNDL", wlimits );
 
649
    cpx = wlimits[XMIN];
 
650
    cpy = wlimits[YMIN];
 
651
 
 
652
    SCTPUT( " " );
 
653
    SCTPUT( "   X         Wave" );
 
654
    SCTPUT( "--------------------" );
 
655
    while ( 1 ) { /* forever */
 
656
        AG_VLOC( &cpx, &cpy, &key, &valpix );
 
657
        if ( key == MID_BUT )
 
658
            break;
 
659
        wave = eval_dpoly( cpx, Coef-1, Fitd+1 );
 
660
        sprintf( msg, "%7.2f    %9.2f", cpx, wave );
 
661
        SCTPUT( msg );
 
662
    }
 
663
    SCTPUT( " " );
 
664
}
 
665
 
 
666
/***********************************************************
 
667
 * read_line_table(): reads all columns of 'LINTAB',
 
668
 *              for values :Y == Ystart.
 
669
 * OUT: 
 
670
 *      
 
671
 ***********************************************************/
 
672
void read_line_table()
 
673
{
 
674
    char erased, erast[3];
 
675
    int  nulval, sortcol, aw, ar, ncols, nrows;
 
676
    int i, id, selected; 
 
677
    int n = 0;
 
678
 
 
679
    NumLine = 0;
 
680
 
 
681
    if ( !file_exists( Lintab, ".tbl" ) ) {
 
682
        SCTPUT( "*** Lines have not been searched ***" );
 
683
        end_midas();
 
684
    }
 
685
    TCTOPN( Lintab, F_IO_MODE, &id );
 
686
    read_column_ids( id );
 
687
 
 
688
    TCIGET( id, &ncols, &nrows, &sortcol, &aw, &ar );
 
689
    
 
690
    for ( i = 1; i <= nrows; i++ ) {
 
691
        TCSGET(id, i, &selected);
 
692
        if ( selected )
 
693
            n++;
 
694
    }
 
695
 
 
696
    init_arrays_data( n );
 
697
 
 
698
    for ( i = 1; i <= nrows; i++ ) {
 
699
        TCSGET(id, i, &selected);
 
700
        if ( selected ) {
 
701
            TCERDR( id, i, Col.x, &X[NumLine], &nulval );
 
702
            TCERDR( id, i, Col.ident, &Ident[NumLine], &nulval );
 
703
            TCERDR( id, i, Col.peak, &Peak[NumLine], &nulval );
 
704
            TCERDR( id, i, Col.wave, &Wave[NumLine], &nulval );
 
705
            TCERDR( id, i, Col.wavec, &Wavec[NumLine], &nulval );
 
706
            TCERDR( id, i, Col.delta, &Delta[NumLine], &nulval );
 
707
            TCERDR( id, i, Col.deltac, &Deltac[NumLine], &nulval );
 
708
            erased = VAL_OK;
 
709
            TCERDC(id, i, Col.erased, erast, &nulval );
 
710
            erased = erast[0];
 
711
            if ( erased == VAL_ERASED )
 
712
                DelList[NumDel++] = NumLine;
 
713
            Row[NumLine] = i;
 
714
            NumLine++;
 
715
        }
 
716
    }
 
717
    TCTCLO(id);
 
718
 
 
719
    if ( NumDel > 0 )
 
720
        compute_calib_values();
 
721
}
 
722
 
 
723
void init_arrays_data( n )
 
724
int n;
 
725
{
 
726
    int i;
 
727
    float *fvector();
 
728
    int *ivector();
 
729
 
 
730
    /* allocate space for data */
 
731
    Row      = ivector( 0, n-1 );
 
732
    X        = fvector( 0, n-1 );
 
733
    Ident    = fvector( 0, n-1 );
 
734
    Peak     = fvector( 0, n-1 );
 
735
    Wave     = fvector( 0, n-1 );
 
736
    Wavec    = fvector( 0, n-1 );
 
737
    Delta    = fvector( 0, n-1 );
 
738
    Deltac   = fvector( 0, n-1 );
 
739
 
 
740
    /* initialize data */
 
741
    for ( i = 0; i < n; i++ )
 
742
        X[i] = Ident[i] = Peak[i] = Wave[i] = Wavec[i] = 
 
743
               Delta[i] = Deltac[i] = Rnull;
 
744
}
 
745
 
 
746
void free_data( n )
 
747
int  n;
 
748
 
 
749
{
 
750
    free_catalog_table(Lc);
 
751
    free_ivector( Row, 0, n-1 );
 
752
    free_fvector( X, 0, n-1 );
 
753
    free_fvector( Ident, 0, n-1 );
 
754
    free_fvector( Peak, 0, n-1 );
 
755
    free_fvector( Wave, 0, n-1 );
 
756
    free_fvector( Wavec, 0, n-1 );
 
757
    free_fvector( Delta, 0, n-1 );
 
758
    free_fvector( Deltac, 0, n-1 );
 
759
}
 
760
 
 
761
void read_column_ids( id )
 
762
int id;
 
763
{
 
764
    TCCSER( id, ":X",       &Col.x );
 
765
    TCCSER( id, ":Y",       &Col.y );
 
766
    TCCSER( id, ":PEAK",    &Col.peak );
 
767
    TCCSER( id, ":IDENT",   &Col.ident );
 
768
    TCCSER( id, ":WAVE",    &Col.wave );
 
769
    TCCSER( id, ":WAVEC",   &Col.wavec );
 
770
    TCCSER( id, ":DELTA",   &Col.delta );
 
771
    TCCSER( id, ":DELTAC",  &Col.deltac );
 
772
    TCCSER( id, ":ERASED",  &Col.erased );
 
773
 
 
774
    if ( Col.x == -1 || Col.y == -1 || Col.peak == -1 || Col.ident == -1 ||
 
775
         Col.wave == -1 || Col.wavec == -1 || 
 
776
         Col.delta == -1 || Col.deltac == -1 ) {
 
777
        SCTPUT( "*** Starting line has not been calibrated ***" );
 
778
        end_midas();
 
779
    }
 
780
 
 
781
    if ( Col.erased == -1 )
 
782
        TCCINI( id, D_C_FORMAT, 1, "A1", " ", "ERASED", &Col.erased );
 
783
}
 
784
 
 
785
 
 
786
void init_midas()
 
787
{
 
788
    SCSPRO( "SPPLOT" );
 
789
    TCMNUL( &Inull, &Rnull, &Dnull ); /* obtain NULL values */
 
790
}
 
791
 
 
792
void end_midas()
 
793
{
 
794
    SCSEPI();
 
795
}
 
796
 
 
797
void read_parameters()
 
798
{
 
799
    int actval;         /* actual values returned */
 
800
    int  unit;          /* useless */
 
801
    int nulval;         /* useless */
 
802
    int pltkey[2];
 
803
 
 
804
    SCKRDI( "SPPLT", 1, 2, &actval, pltkey, &unit, &nulval );
 
805
    SCKRDI( "DCX", 1, 1, &actval, &Fitd, &unit, &nulval );
 
806
    SCKRDI( "YSTART", 1, 1, &actval, &Ystart, &unit, &nulval );
 
807
    SCKRDR( "IMIN", 1, 1, &actval, &Imin, &unit, &nulval );
 
808
    SCKRDI( "WRANG", 1, 2, &actval, Wrang, &unit, &nulval );
 
809
    SCKGETC( "WLC", 1, 20, &actval, Wlc );
 
810
    SCKGETC( "LINTAB", 1, 20, &actval, Lintab );
 
811
    SCKGETC( "LINCAT", 1, 20, &actval, Lincat );
 
812
    SCKGETC( "COERBR", 1, 20, &actval, Coerbr );
 
813
    SCKGETC( "INPUTC", 1, 20, &actval, PlotFile );
 
814
    PlotAction = pltkey[0];
 
815
    PlotType   = pltkey[1];
 
816
}
 
817
 
 
818
/******************************************************
 
819
 * save_limits() : Save current viewport coordinates.
 
820
 ******************************************************/
 
821
void save_limits( xmin, xmax, ymin, ymax )
 
822
float xmin, xmax, ymin, ymax;
 
823
{
 
824
    int  unit;
 
825
    float lims[4];
 
826
 
 
827
    lims[0] = xmin;
 
828
    lims[1] = xmax;
 
829
    lims[2] = ymin;
 
830
    lims[3] = ymax;
 
831
    SCKWRR( "AGLIMS", lims, 1, 4, &unit );
 
832
}
 
833
 
 
834
/********************************************************
 
835
 * read_limits() : Read last saved viewport coordinates.
 
836
 ********************************************************/
 
837
void read_limits( xmin, xmax, ymin, ymax )
 
838
float *xmin, *xmax, *ymin, *ymax;
 
839
{
 
840
    int  unit;
 
841
    float lims[4];
 
842
    int  actval, nulval;        /* garbage */
 
843
 
 
844
    SCKRDR( "AGLIMS", 1, 4, &actval, lims, &unit, &nulval );
 
845
    *xmin = lims[0];
 
846
    *xmax = lims[1];
 
847
    *ymin = lims[2];
 
848
    *ymax = lims[3];
 
849
}
 
850
 
 
851
/*************************************************************
 
852
 * init_viewport() : Init viewport according to the last one.
 
853
 *************************************************************/
 
854
void init_viewport()
 
855
{
 
856
    float xmin, xmax, ymin, ymax;
 
857
 
 
858
    read_limits( &xmin, &xmax, &ymin, &ymax ); 
 
859
    AG_AXES( xmin, xmax, ymin, ymax, " " );
 
860
    AG_VUPD();
 
861
}
 
862
 
 
863
 
 
864
#define FMT_TITLE       "File: %s  Line: %d  Image: %s"
 
865
#define FMT_OPTIONS     "TITLE=%s;LABX=Position;LABY=Pixel value"
 
866
 
 
867
#define ORDMAX          9999999.0
 
868
#define ORDMIN          -9999999.0
 
869
 
 
870
/********************************************************
 
871
 * read_image(): plots the 'row'nth line of 'image'.
 
872
 ********************************************************/
 
873
void read_image( row, image )
 
874
int row;
 
875
char *image;
 
876
{
 
877
    int      i, framid;
 
878
    int ncols, nrows, npix[2];
 
879
    double start, step;
 
880
    float cuts[4];
 
881
    int  unit;
 
882
    int  nulval, retval; /* garbage */
 
883
    char options[512], title[512];
 
884
    float x[MAXDATA], y[MAXDATA];
 
885
    char ident[MAXIDENT+1];
 
886
    float xmin, xmax, ymin, ymax;
 
887
 
 
888
    SCFOPN( image, D_R4_FORMAT, 0, F_IMA_TYPE, &framid );
 
889
    SCDRDI( framid, "NPIX", 1, 2, &retval, npix, &unit, &nulval);
 
890
    ncols = npix[0];
 
891
    nrows = npix[1];
 
892
 
 
893
    SCDRDR( framid, "LHCUTS", 1, 4, &retval, cuts, &unit, &nulval );
 
894
    SCDRDD( framid, "START", 1, 1, &retval, &start, &unit, &nulval );
 
895
    SCDRDD( framid, "STEP", 1, 1, &retval,  &step, &unit, &nulval );
 
896
    SCDGETC( framid, "IDENT", 1, MAXIDENT, &retval, ident );
 
897
 
 
898
    /* mapping 'row' row (Y axis) */
 
899
    SCFGET( framid, (row-1)*ncols+1, ncols, &retval, (char *)y );
 
900
 
 
901
    for ( i = 0; i < ncols; i++ )
 
902
        x[i] = start + step*i;
 
903
 
 
904
    xmin = start;
 
905
    xmax = start + step*(ncols - 1);
 
906
 
 
907
    if ( cuts[1] != 0.0 ) {
 
908
        ymin = cuts[0]; 
 
909
        ymax = cuts[1];
 
910
    }
 
911
    else if ( cuts[3] != 0.0 ) {
 
912
        ymin = cuts[2]; 
 
913
        ymax = cuts[3];
 
914
    }
 
915
    else {      /* Max. and min. cuts are not assigned */
 
916
        ymin = ORDMAX;
 
917
        ymax = ORDMIN;
 
918
        for ( i = 0; i < ncols; i++ ) {
 
919
            if ( y[i] > ymax )
 
920
                ymax = y[i];
 
921
            if ( y[i] < ymin )
 
922
                ymin = y[i];
 
923
        }
 
924
        cuts[0] = ymin;
 
925
        cuts[1] = ymax;
 
926
        SCDWRR( framid, "LHCUTS", cuts, 3, 2, &unit );
 
927
    }
 
928
    SCFCLO(framid);
 
929
 
 
930
    sprintf( title, FMT_TITLE, image, row, ident );
 
931
    sprintf( options, FMT_OPTIONS, title );
 
932
 
 
933
    AG_VERS();
 
934
    AG_AXES( xmin, xmax, ymin, ymax, options );
 
935
 
 
936
    AG_GPLL( x, y, ncols );
 
937
    AG_VUPD();
 
938
 
 
939
    save_limits( xmin, xmax, ymin, ymax ); 
 
940
}
 
941
 
 
942
void init_graphic( devtype )
 
943
int devtype;
 
944
{
 
945
    if ( !graphwin_exists() ) {
 
946
        SCTPUT( "*** Please create the graphic window ***" );
 
947
        end_midas();
 
948
    }
 
949
 
 
950
    if ( InitGraphic ) 
 
951
        return;
 
952
 
 
953
    InitGraphic = TRUE;
 
954
    get_agldev();
 
955
    switch ( devtype ) {
 
956
        case DEV_ERASE:
 
957
                AG_VDEF( DevErase, 0.05, 1.0, 0.0, 1.0, 0.0, 0.0 );
 
958
                break;
 
959
        case DEV_NO_ERASE:
 
960
                AG_VDEF( DevNoErase, 0.05, 1.0, 0.0, 1.0, 0.0, 0.0 );
 
961
                break;
 
962
    }
 
963
}
 
964
 
 
965
void end_graphic()
 
966
{
 
967
    if ( InitGraphic && graphwin_exists() )
 
968
        AG_CLS();
 
969
 
 
970
    InitGraphic = FALSE;
 
971
}
 
972
 
 
973
/*********************************************************
 
974
 * get_agldev(): translate IDI device to devices erasable
 
975
 * and non-erasable suitable for AG_VDEF calls.
 
976
 *********************************************************/
 
977
void get_agldev()
 
978
{
 
979
    char device[21];    /* name of assoc. device as in pltdevices.dat */
 
980
 
 
981
    /* read & translate type of device to a specific device name */
 
982
    /* SCKGETC( "MID$PLOT", 1, 20, &actval, devkeyw );
 
983
       get_dev( devkeyw, device ); Not by now... */
 
984
    strcpy( device, "GRAPH_WND0" );
 
985
 
 
986
    /* now make the AGL device names */
 
987
    strcpy( DevErase, device ); 
 
988
    strcat( DevErase, ":" ); /* see AG_VDEF definition */
 
989
 
 
990
    strcpy( DevNoErase, device ); 
 
991
    strcat( DevNoErase, "/n:" ); /* see AG_VDEF definition */
 
992
}
 
993
 
 
994
void dpoly( x, p, np )  /* (p[i] <-- x**(i-1)), i = 1,np) */
 
995
double x, p[];
 
996
int np;
 
997
{
 
998
    int j;
 
999
 
 
1000
    p[1] = 1.0;
 
1001
    for ( j = 2; j <= np; j++ )
 
1002
        p[j] = p[j-1] * x;
 
1003
}
 
1004
 
 
1005
void read_rebin_parameters()
 
1006
{
 
1007
    int diff, min_diff = 32767;  /* :-) */
 
1008
    int i, id;
 
1009
    int col_y, col_rms, col_pixel;
 
1010
    int  nulval;        /* useless */
 
1011
    double y, pixel, rms;
 
1012
    int  sortcol, aw, ar, ncols, nrows;
 
1013
    
 
1014
 
 
1015
    if ( ! file_exists(Coerbr, ".tbl") ) {
 
1016
        SCTPUT( "Coefficients table couldn't be opened. Stop." );
 
1017
        end_midas();
 
1018
    }
 
1019
    TCTOPN(Coerbr, F_IO_MODE, &id);
 
1020
    TCIGET( id, &ncols, &nrows, &sortcol, &aw, &ar );
 
1021
 
 
1022
    if ( nrows == 0 ) {
 
1023
        SCTPUT( "Error: coefficients table is empty." );
 
1024
        end_midas();
 
1025
    }
 
1026
 
 
1027
    TCCSER( id, ":ROW", &col_y );
 
1028
    TCCSER( id, ":RMS", &col_rms );
 
1029
    TCCSER( id, ":PIXEL", &col_pixel );
 
1030
    if ( col_y == -1 || col_rms == -1 || col_pixel == -1 ) {
 
1031
        SCTPUT( "Calibration process has not been performed. Stop." );
 
1032
        end_midas();
 
1033
    }
 
1034
 
 
1035
    for ( i = 1; i <= nrows; i++ ) {
 
1036
        TCERDD( id, i, col_y, &y, &nulval );
 
1037
        TCERDD( id, i, col_pixel, &pixel, &nulval );
 
1038
        TCERDD( id, i, col_rms, &rms, &nulval );
 
1039
        diff = fabs(Ystart - y);
 
1040
        if ( diff < min_diff ) {
 
1041
            min_diff = diff;
 
1042
            Ycoerbr = i;
 
1043
            Rms = rms;
 
1044
            Dispersion = pixel;
 
1045
        }
 
1046
    }
 
1047
    TCTCLO( id );
 
1048
}
 
1049
 
 
1050
void read_coefs_ystart()
 
1051
{
 
1052
    int  nulval;        /* useless */
 
1053
    double value;
 
1054
    int i, id;
 
1055
 
 
1056
    TCTOPN(Coerbr, F_IO_MODE, &id);
 
1057
    for ( i = 3; i <= Fitd + 3; i++ ) {
 
1058
        TCERDD( id, Ycoerbr, i, &value, &nulval );
 
1059
        Coef[i-3] = value;
 
1060
    }
 
1061
    TCTCLO(id);
 
1062
}
 
1063
 
 
1064
void compute_calib_values()
 
1065
{
 
1066
    double *a; /* coefficients for fit = 1 */
 
1067
    double val;
 
1068
    double *xfit, *yfit;
 
1069
    double delta, delta_min;
 
1070
    int id, i, j, k, n = 1;
 
1071
 
 
1072
    a = dvector( 1, Fitd+1 );
 
1073
    xfit = dvector( 1, NumLine );
 
1074
    yfit = dvector( 1, NumLine );
 
1075
 
 
1076
    for ( k = 0; k < NumDel; k++ ) {
 
1077
        i = DelList[k];
 
1078
        Wavec[i] = eval_dpoly( X[i], Coef-1, Fitd+1 );
 
1079
        delta_min = MAXWAVE;
 
1080
        for ( j = 0; j < Lc->nrows; j++ )
 
1081
            if ( (delta = fabs(Lc->wave[j] - Wavec[i])) < delta_min &&
 
1082
                 ! Lc->sel[j] ) {
 
1083
                Wave[i] = Lc->wave[j];
 
1084
                delta_min = delta;
 
1085
            }
 
1086
    }
 
1087
 
 
1088
    for ( j = 0; j < NumLine; j++ )
 
1089
        if ( Wave[j] != Rnull ) {
 
1090
            xfit[n] = X[j];
 
1091
            yfit[n] = Wave[j];
 
1092
            n++;
 
1093
        }
 
1094
 
 
1095
    /* get the coefficients of the relation w=a0+a1x,
 
1096
       calculate Delta and Deltac */
 
1097
    lfit( xfit, yfit, n - 1, a, 2, dpoly ); /* fitting */
 
1098
    for ( k = 0; k < NumDel; k++ ) {
 
1099
        i = DelList[k];
 
1100
        val = a[1] + a[2] * X[i];
 
1101
        Deltac[i] = val - Wavec[i];
 
1102
        Delta[i] = val - Wave[i];
 
1103
    }
 
1104
    
 
1105
    /* now save the wave values in the table (for lnerase.exe) */
 
1106
    TCTOPN( Lintab, F_IO_MODE, &id );
 
1107
    for ( i = 0; i < NumDel; i++ )
 
1108
        TCEWRR( id, Row[DelList[i]], Col.wave, &Wave[DelList[i]] );
 
1109
    TCTCLO(id);
 
1110
 
 
1111
    free_dvector( xfit, 1, NumLine );
 
1112
    free_dvector( yfit, 1, NumLine );
 
1113
}
 
1114
 
 
1115
/*****************************************************************
 
1116
 * read_lincat_table(): reads line catalog named by 'Lincat' and
 
1117
 * stores the values in the 'Lc' structure, under the
 
1118
 * restrictions Wrang & Imin.
 
1119
 *****************************************************************/
 
1120
int read_lincat_table()
 
1121
{
 
1122
    if ( ! file_exists( Lincat, ".tbl" ) ) {
 
1123
        SCTPUT( "*** Line catalogue doesn't exist ***" );
 
1124
        return(FALSE);
 
1125
    }
 
1126
 
 
1127
    if ( Lc != NULL )
 
1128
        free_catalog_table(Lc);
 
1129
    Lc  = (LCTAB *)osmmget( sizeof(LCTAB) );
 
1130
 
 
1131
    /* call to a library routine */
 
1132
    if( !read_catalog_table( Lc, Lincat, Wrang, Imin ) ) {
 
1133
        Lc = NULL;
 
1134
        return FALSE;
 
1135
    }
 
1136
    return TRUE;
 
1137
}