/* @(#)lnplot.c	17.1.1.1 (ES0-DMD) 01/25/02 17:54:00 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/* @(#)lnplot.c	17.1.1.1  (ESO)  01/25/02  17:54:00 */
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* .COPYRIGHT   (C) 1993 European Southern Observatory     */
/* .IDENT       lnplot.c                                   */
/* .AUTHORS     Pascal Ballester (ESO/Garching)            */
/*              Cristian Levin   (ESO/La Silla)            */
/* .KEYWORDS    Spectroscopy, Long-Slit                    */
/* .PURPOSE                                                */
/* .VERSION     1.0  Package Creation  17-MAR-1993         */
/* ------------------------------------------------------- */

#include <stdio.h>
#include <math.h>

#include <midas_def.h>

#include <gl_defs.h>
#include <main_defs.h>
#include <spec_defs.h>
#include <agl_defs.h>

#define MINWAVE		-99999.0
#define MAXWAVE		99999.0
#define MINY		-99999.0
#define MAXY		99999.0
#define DBORDER		10.0
#define NPOINTS		500
#define EPS             0.001

/* Actions */
#define PLOT		0
#define DELETE		1
#define GETCUR		2
#define PRINT		3

#define MAXDEL		80		/* max. temp. deleted points */

/* Keywords */
char 	Lincat[MAXLINE];
char	Wlc[21];
int	Wrang[2];
float 	Imin;
int	Ystart;
int	Fitd;
int 	PlotType;	/* type of current plot */
int 	PlotAction;	/* Action to do */

/* NULL values */
int    	Inull;
float  	Rnull;
double 	Dnull;

/* structures for allocating column values of 'LINTAB' */
float *X, *Peak, *Ident, *Wave, *Wavec, *Delta, *Deltac;
int *Row;	/* row number associated to each value */

int  	NumLine;		/* Number of rows read from 'LINTAB' */
float 	Wlimits[4];		/* window limits */

int	DelList[MAXDEL];	/* List of temporarily deleted points */
int	NumDel = 0;		/* DelList[] cardinality */

char 	CurrFile[81];		/* Global variable used for plot spectra */
char	DevErase[80], DevNoErase[80];

char 	PlotFile[MAXLINE], Lintab[MAXLINE], Coerbr[MAXLINE];
int 	Ycoerbr; 	/* row nearest to :ystart in 'COERBR' */
double	Rms, Dispersion;
double	Coef[20];

LCTAB * Lc = NULL;

int 	InitGraphic = FALSE;

struct column_ids {
	int x, y, peak, ident;
	int wave, wavec, delta, deltac;
	int erased;
} Col;

void compute_calib_values();
void read_coefs_ystart();
void read_rebin_parameters();
void read_column_ids();
void free_data();
void init_arrays_data();
void read_line_table();
void read_parameters();
void init_midas();
void end_midas();
void get_agldev();
void end_graphic();
void init_graphic();
void read_image();
void init_viewport();
void read_limits();
void save_limits();

void getcur_wave();
void plot_wave();
void plot_ident();
void plot_splabel();
void init_ident();
void init_splabel();
void plot_curve();
void plot_delete();
void save_limits();
void read_limits();

void del_ident();
void del_point();
void undel_point();

float *fvector();
double *dvector();
double eval_dpoly();
void dpoly();

main()
{
    init_midas();

    read_parameters();
    read_lincat_table();
    read_rebin_parameters();
    read_coefs_ystart();
    read_line_table();

    if ( PlotAction == PLOT ) {
        init_graphic(DEV_ERASE);
	AG_MOPN( PlotFile );
        AG_SSET("FONT=1");     /* nice font */
    	switch( PlotType ) {
	    case PLOT_IDENT: 	plot_ident(); break;
	    case PLOT_WAVE: 	plot_wave(); break;
	    case PLOT_SPLABEL: 	plot_splabel(); break;
    	}
	AG_MCLS();
    }
    else if ( PlotAction == DELETE ) {
        init_graphic(DEV_NO_ERASE);
        AG_SSET("FONT=1");
    	switch( PlotType ) {
	    case PLOT_IDENT: 	
	    case PLOT_WAVE: 	
	    case PLOT_SPLABEL: 	init_viewport(); break;
    	}
    }
    else if ( PlotAction == GETCUR ) {
        init_graphic(DEV_NO_ERASE);
        AG_SSET("FONT=1");
	init_viewport();
	getcur_wave();
    }

    if ( PlotAction == DELETE && PlotType != PLOT_NONE )
	plot_delete();

    free_data();
    end_graphic();
    end_midas();
}

void plot_delete()
{
    float cpx = 0.0, cpy = 0.0;
    int key, valpix;
    int i = 0, min = 0;
    float min_deltax, min_deltay;
    float x[2], y[2], xi, yi;
    char str[80];

    while (1) {
    	min_deltax =  min_deltay = MAXWAVE;

    	AG_VLOC( &cpx, &cpy, &key, &valpix );
    	if ( key == MID_BUT )
	    break;

        AG_SSET( RED_COLOR );
    	switch( PlotType ) {
	   case PLOT_IDENT:
    	        for ( i = 0; i < NumLine; i++ ) {
		    xi = Wavec[i]; 
	            yi = Deltac[i] + Wavec[i] - Ident[i];
	            if ( Wave[i] != Rnull &&
		         fabs((double) cpx - xi) < min_deltax &&
		         fabs((double) cpy - yi) < min_deltay ) {
	        	min_deltax = fabs((double) cpx - xi);
	        	min_deltay = fabs((double) cpx - yi);
			min = i;
		    }
		}
		x[0] = Wavec[min]; 
	        y[0] = Deltac[min] + Wavec[min] - Ident[min];
		del_ident(min);
		sprintf(str, "point deleted : %10.3f %10.3f", x[0], y[0]);
		SCTPUT(str);
    		AG_GPLM( x, y, 1, BOX_MARKER );
		break;
	   case PLOT_WAVE:
    	        for ( i = 0; i < NumLine; i++ ) {
		    xi = Wave[i]; 
	            yi = Delta[i];
	            if ( Wave[i] != Rnull &&
		         fabs((double) cpx - xi) < min_deltax &&
		         fabs((double) cpy - yi) < min_deltay ) {
	        	min_deltax = fabs((double) cpx - xi);
	        	min_deltay = fabs((double) cpx - yi);
			min = i;
		    }
		}
		x[0] = Wave[min]; 
	        y[0] = Delta[min];
		if ( point_deleted(min) ) {
        	    AG_SSET( BLUE_COLOR );
		    undel_point(min);
		    sprintf(str, "point added   : %10.3f %10.3f", x[0], y[0]);
		}
		else {
		    del_point(min);
		    sprintf(str, "point deleted : %10.3f %10.3f", x[0], y[0]);
		}
		SCTPUT(str);
    		AG_GPLM( x, y, 1, X_MARKER );
		break;
	   case PLOT_SPLABEL:
    	        for ( i = 0; i < NumLine; i++ ) {
	            if ( Wave[i] != Rnull &&
		         fabs((double) cpx - X[i]) < min_deltax ) {
	        	min_deltax = fabs((double) cpx - X[i]);
			min = i;
		    }
		}
    		x[0] = x[1] = X[min];
    		AG_RGET( "WNDL", Wlimits );
    		y[0] = Wlimits[YMIN];
    		y[1] = Peak[min];
		if ( point_deleted(min) ) {
        	    AG_SSET( ((Ident[min] != Rnull) ? GREEN_COLOR:BLUE_COLOR) );
		    undel_point(min);
		    sprintf(str, "point added   : %10.3f", x[0]);
		}
		else {
		    del_point(min);
		    sprintf(str, "point deleted : %10.3f", x[0]);
		}
		SCTPUT(str);
    		AG_GPLL( x, y, 2 );
		break;
	}
    	AG_VUPD();
        AG_SSET( DEF_COLOR );
    }
    NumDel = 0; /* no temporarily deleted points */
}

void del_ident(n)
int n;
{
    int id;

    TCTOPN( Lintab, F_IO_MODE, &id );
    TCEWRR( id, Row[n], Col.ident, &Rnull );
    TCTCLO(id);
}

void del_point(n)
int n;
{
    int id;
    char del_str[2];

    sprintf(del_str, "%c", VAL_ERASED);
    DelList[NumDel++] = n;

    TCTOPN( Lintab, F_IO_MODE, &id );
    TCEWRC( id, Row[n], Col.erased, del_str );
    TCTCLO(id);
}

void undel_point(n)
int n;
{
    int  select = TRUE;
    int i, id, min = -1;
    char undel_str[2];
    float min_delta = MAXWAVE;

    sprintf(undel_str, "%c", VAL_OK);

    TCTOPN( Lintab, F_IO_MODE, &id );
    for ( i = 0; i < NumDel; i++ )
	if ( n == DelList[i] ) {
	    DelList[i] = DelList[NumDel-1];
    	    TCEWRC( id, Row[n], Col.erased, undel_str );
	    NumDel--;
	    break;
	}
    TCTCLO(id);

    TCTOPN( Lincat, F_IO_MODE, &id );
    for ( i = 0; i < Lc->nrows; i++ )
        if ( fabs(Lc->wave[i] - Wave[n]) < min_delta && ! Lc->sel[i] ) {
	    min = i;
	    min_delta = fabs(Lc->wave[i] - Wave[n]);
	}
    if ( min != -1 )
    	TCSPUT(id, Lc->row[min], &select);
    TCTCLO(id);
}

int point_deleted(n)
int n;
{
    int i;

    for ( i = 0; i < NumDel; i++ )
	if ( n == DelList[i] )
	    return( TRUE );
    return( FALSE );
}

void plot_wave()
{
    int i, n1 = 0, n2 = 0;
    int id, ncoef;
    int ndel = 0;
    float *xdel, *ydel;
    float xmin = MAXWAVE, xmax = MINWAVE;
    float ymin = MAXY, ymax = MINY;
    float *x1, *x2, *y1, *y2;
    char options[MAXOPTS];
    int nulval, actval; /* useless */
    int  unit;		/* useless */
    char text[80];

    /*** get the values to plot ***/

    x1 = fvector( 0, NumLine-1 );
    y1 = fvector( 0, NumLine-1 );
    x2 = fvector( 0, NumLine-1 );
    y2 = fvector( 0, NumLine-1 );
    xdel = fvector( 0, NumLine-1 );
    ydel = fvector( 0, NumLine-1 );

    for ( i = 0; i < NumLine; i++ ) {
	if ( Ident[i] != Rnull ) { 	/* get :IDENT points */
	    x1[n1] = Wavec[i];
	    y1[n1] = Deltac[i] + Wavec[i] - Ident[i];
	    if ( x1[n1] < xmin ) xmin = x1[n1];
	    if ( x1[n1] > xmax ) xmax = x1[n1];
	    if ( y1[n1] < ymin ) ymin = y1[n1];
	    if ( y1[n1] > ymax ) ymax = y1[n1];
	    n1++;
	}
	if ( point_deleted(i) ) {	/* get deleted points & continue */
	    xdel[ndel] = Wave[i];
	    ydel[ndel] = Delta[i];
	    if ( xdel[ndel] < xmin ) xmin = xdel[ndel];
	    if ( xdel[ndel] > xmax ) xmax = xdel[ndel];
	    if ( ydel[ndel] < ymin ) ymin = ydel[ndel];
	    if ( ydel[ndel] > ymax ) ymax = ydel[ndel];
	    ndel++;
	    continue;
        }
	if ( Wave[i] != Rnull ) {	/* get :WAVE points not deleted */
	    x2[n2] = Wave[i];
	    y2[n2] = Delta[i];
	    if ( x2[n2] < xmin ) xmin = x2[n2];
	    if ( x2[n2] > xmax ) xmax = x2[n2];
	    if ( y2[n2] < ymin ) ymin = y2[n2];
	    if ( y2[n2] > ymax ) ymax = y2[n2];
	    n2++;
	}
    }

    /*** plot... ***/

    /* let's give a border for a nicer plot */
    xmin -= fabs( (double) (xmax - xmin) / DBORDER );
    xmax += fabs( (double) (xmax - xmin) / DBORDER );
    ymin -= fabs( (double) (ymax - ymin) / DBORDER );
    ymax += fabs( (double) (ymax - ymin) / DBORDER );

    strcpy( options, "LABY=Delta(Wave);LABX=Wavelength" );
    AG_AXES( xmin, xmax, ymin, ymax, options );
    
    /* Plot the :IDENT related points */
    if (n1 > 0) {
        AG_SSET( GREEN_COLOR );
        AG_GPLM( x1, y1, n1, BOX_MARKER );
        AG_VUPD();
        AG_SSET( DEF_COLOR );
      }

    /* Plot the :WAVEC related points */
    AG_SSET( BLUE_COLOR );
    AG_GPLM( x2, y2, n2, X_MARKER );
    AG_VUPD();
    AG_SSET( DEF_COLOR );

    /* Plot the deleted points */
    if ( ndel > 0 ) {
    	AG_SSET( RED_COLOR );
    	AG_GPLM( xdel, ydel, ndel, X_MARKER );
    	AG_VUPD();
    	AG_SSET( DEF_COLOR );
    }

    /* Plot horizontal line at y = 0 */
    AG_SSET( DASH_LINE );
    x1[0] = xmin; x1[1] = xmax;
    y1[0] = y1[1] = 0.0;
    AG_GPLL( x1, y1, 2 );
    AG_VUPD();
    AG_SSET( SOLID_LINE );

    ncoef = (n2 > Fitd) ? Fitd+1 : n2; /* coeffs. to obtain */
    plot_curve( x2, y2, n2, xmin, xmax, ncoef );

    AG_RGET( "WNDL", Wlimits );
    AG_SSET( LR_TEXT );
    sprintf( text, "RMS=%.3f   DISP=%.3f", Rms, Dispersion );
    AG_GTXT( Wlimits[XMAX], Wlimits[YMAX], text, 18 );
    AG_VUPD();

    free_fvector( x1, 0, NumLine-1 );
    free_fvector( y1, 0, NumLine-1 );
    free_fvector( x2, 0, NumLine-1 );
    free_fvector( y2, 0, NumLine-1 );
    free_fvector( xdel, 0, NumLine-1 );
    free_fvector( ydel, 0, NumLine-1 );

    save_limits( xmin, xmax, ymin, ymax ); /* for init_viewport() */
}

void plot_curve( x, y, n, xmin, xmax, ncoef )
float *x, *y;
float xmin, xmax;
int n, ncoef;
{
    float *px, *py, dx, step;
    double *a; /* coefficients for fitting */
    double *p, *xfit, *yfit;
    int np, intervals;
    int i, j;

    a = dvector( 1, ncoef );
    p = dvector( 1, ncoef );
    xfit = dvector( 1, n );
    yfit = dvector( 1, n );

    intervals = NPOINTS + 2;
    step = (xmax - xmin) / NPOINTS;
    px = fvector( 0, intervals-1 );
    py = fvector( 0, intervals-1 );

    for ( i = 1; i <= n; i++ ) {
	xfit[i] = x[i-1];
	yfit[i] = y[i-1];
    }
    lfit( xfit, yfit, n, a, ncoef, dpoly ); /* fitting */
   
    np = 0;
    for ( dx = xmin; dx <= xmax; dx += step ) {
	px[np] = dx;
        py[np] = eval_dpoly(dx, a, ncoef);
	np++;
    }
    AG_GPLL( px, py, np );
    AG_VUPD();

    free_fvector( px, 0, intervals-1 );
    free_fvector( py, 0, intervals-1 );
    free_dvector( a, 1, Fitd+1 );
    free_dvector( p, 1, Fitd+1 );
    free_dvector( xfit, 1, n );
    free_dvector( yfit, 1, n );
}

void plot_splabel()
{
    char text[80];
    int i;
    float xdel[2], ydel[2];

    if ( !file_exists( Wlc, ".bdf" ) ) {
	SCTPUT( "*** Calibration image doesn't exist ***" );
	end_midas();
    }
    read_image( Ystart, Wlc );

    AG_RGET( "WNDL", Wlimits );
    ydel[0] = Wlimits[YMIN];

    AG_SSET( DU_TEXT );
    AG_SSET( DEF_CHAR );
    AG_SSET( HIGH_FONT );
    AG_SSET( BLUE_COLOR );
    for ( i = 0; i < NumLine; i++ ) {
	if ( Wave[i] != Rnull ) {
	    sprintf( text, "%.1f", Wave[i] );
	    AG_GTXT( X[i], Wlimits[YMAX], text, 17 );
	    AG_VUPD();
	}
	if ( point_deleted(i) || Ident[i] != Rnull ) {
    	    AG_SSET( (point_deleted(i) ? RED_COLOR : GREEN_COLOR) );
	    xdel[0] = xdel[1] = X[i];
	    ydel[1] = Peak[i];
    	    AG_GPLL( xdel, ydel, 2 );
	    AG_VUPD();
    	    AG_SSET( BLUE_COLOR );
        }
    }
    AG_SSET( LR_TEXT );
    AG_SSET( DEF_FONT );
    AG_SSET( DEF_COLOR );
}

void init_ident()
{
    int i;
    float xmin = MAXWAVE, xmax = MINWAVE;
    float ymin = MAXY, ymax = MINY;
    float x, y;
    char options[MAXOPTS];

    /*** get max. and min. values to plot ***/

    for ( i = 0; i < NumLine; i++ )
	if ( Ident[i] != Rnull ) {
	    x = Wavec[i];
	    y = Deltac[i] + Wavec[i] - Ident[i];
	    if ( x < xmin ) xmin = x;
	    if ( x > xmax ) xmax = x;
	    if ( y < ymin ) ymin = y;
	    if ( y > ymax ) ymax = y;
	}

    /* let's give a border for a nicer plot */
    xmin -= fabs( (double) (xmax - xmin) / DBORDER );
    xmax += fabs( (double) (xmax - xmin) / DBORDER );
    ymin -= fabs( (double) (ymax - ymin) / DBORDER );
    ymax += fabs( (double) (ymax - ymin) / DBORDER );

    strcpy( options, "LABY=Delta(Ident);LABX=Wavelength" );
    AG_AXES( xmin, xmax, ymin, ymax, options );
    AG_VUPD();
}

void plot_ident()
{
    int i, n = 0;
    float xmin = MAXWAVE, xmax = MINWAVE;
    float ymin = MAXY, ymax = MINY;
    float *x, *y;
    char options[MAXOPTS];

    /*** get the values to plot ***/

    x = fvector( 0, NumLine-1 );
    y = fvector( 0, NumLine-1 );
    for ( i = 0; i < NumLine; i++ )
	if ( Ident[i] != Rnull ) {
	    x[n] = Wavec[i];
	    y[n] = Deltac[i] + Wavec[i] - Ident[i];
	    if ( x[n] < xmin ) xmin = x[n];
	    if ( x[n] > xmax ) xmax = x[n];
	    if ( y[n] < ymin ) ymin = y[n];
	    if ( y[n] > ymax ) ymax = y[n];
	    n++;
	}

    /*** plot... ***/

    /* let's give a border for a nicer plot */
    xmin -= fabs( (double) (xmax - xmin) / DBORDER );
    xmax += fabs( (double) (xmax - xmin) / DBORDER );
    ymin -= fabs( (double) (ymax - ymin) / DBORDER );
    ymax += fabs( (double) (ymax - ymin) / DBORDER );

    strcpy( options, "LABY=Delta(Ident);LABX=Wavelength" );
    AG_AXES( xmin, xmax, ymin, ymax, options );
    
    /* Plot the points */
    AG_SSET( BLUE_COLOR );
    AG_GPLM( x, y, n, BOX_MARKER );
    AG_VUPD();
    AG_SSET( DEF_COLOR );

    /* Plot horizontal line at y = 0 */
    AG_SSET( DASH_LINE );
    x[0] = xmin; x[1] = xmax;
    y[0] = y[1] = 0.0;
    AG_GPLL( x, y, 2 );
    AG_VUPD();
    AG_SSET( SOLID_LINE );

    free_fvector( x, 0, NumLine-1 );
    free_fvector( y, 0, NumLine-1 );

    save_limits( xmin, xmax, ymin, ymax ); /* for init_viewport() */
}

void getcur_wave()
{
    int  unit;		/* useless */
    int actval, nulval;	/* useless */
    float cpx, cpy;
    float wlimits[4];	/* window limits */
    int i, id, key, valpix;
    char msg[30];
    double x, wave;

    AG_RGET( "WNDL", wlimits );
    cpx = wlimits[XMIN];
    cpy = wlimits[YMIN];

    SCTPUT( " " );
    SCTPUT( "   X         Wave" );
    SCTPUT( "--------------------" );
    while ( 1 ) { /* forever */
    	AG_VLOC( &cpx, &cpy, &key, &valpix );
	if ( key == MID_BUT )
	    break;
	wave = eval_dpoly( cpx, Coef-1, Fitd+1 );
	sprintf( msg, "%7.2f    %9.2f", cpx, wave );
    	SCTPUT( msg );
    }
    SCTPUT( " " );
}

/***********************************************************
 * read_line_table(): reads all columns of 'LINTAB',
 *		for values :Y == Ystart.
 * OUT: 
 *      
 ***********************************************************/
void read_line_table()
{
    char erased, erast[3];
    int  nulval, sortcol, aw, ar, ncols, nrows;
    int i, id, selected; 
    int n = 0;

    NumLine = 0;

    if ( !file_exists( Lintab, ".tbl" ) ) {
	SCTPUT( "*** Lines have not been searched ***" );
	end_midas();
    }
    TCTOPN( Lintab, F_IO_MODE, &id );
    read_column_ids( id );

    TCIGET( id, &ncols, &nrows, &sortcol, &aw, &ar );
    
    for ( i = 1; i <= nrows; i++ ) {
        TCSGET(id, i, &selected);
        if ( selected )
            n++;
    }

    init_arrays_data( n );

    for ( i = 1; i <= nrows; i++ ) {
        TCSGET(id, i, &selected);
        if ( selected ) {
	    TCERDR( id, i, Col.x, &X[NumLine], &nulval );
	    TCERDR( id, i, Col.ident, &Ident[NumLine], &nulval );
	    TCERDR( id, i, Col.peak, &Peak[NumLine], &nulval );
	    TCERDR( id, i, Col.wave, &Wave[NumLine], &nulval );
	    TCERDR( id, i, Col.wavec, &Wavec[NumLine], &nulval );
	    TCERDR( id, i, Col.delta, &Delta[NumLine], &nulval );
	    TCERDR( id, i, Col.deltac, &Deltac[NumLine], &nulval );
	    erased = VAL_OK;
	    TCERDC(id, i, Col.erased, erast, &nulval );
            erased = erast[0];
	    if ( erased == VAL_ERASED )
    		DelList[NumDel++] = NumLine;
	    Row[NumLine] = i;
	    NumLine++;
	}
    }
    TCTCLO(id);

    if ( NumDel > 0 )
    	compute_calib_values();
}

void init_arrays_data( n )
int n;
{
    int i;
    float *fvector();
    int *ivector();

    /* allocate space for data */
    Row      = ivector( 0, n-1 );
    X        = fvector( 0, n-1 );
    Ident    = fvector( 0, n-1 );
    Peak     = fvector( 0, n-1 );
    Wave     = fvector( 0, n-1 );
    Wavec    = fvector( 0, n-1 );
    Delta    = fvector( 0, n-1 );
    Deltac   = fvector( 0, n-1 );

    /* initialize data */
    for ( i = 0; i < n; i++ )
	X[i] = Ident[i] = Peak[i] = Wave[i] = Wavec[i] = 
	       Delta[i] = Deltac[i] = Rnull;
}

void free_data( n )
{
    free_catalog_table(Lc);
    free_ivector( Row, 0, n-1 );
    free_fvector( X, 0, n-1 );
    free_fvector( Ident, 0, n-1 );
    free_fvector( Peak, 0, n-1 );
    free_fvector( Wave, 0, n-1 );
    free_fvector( Wavec, 0, n-1 );
    free_fvector( Delta, 0, n-1 );
    free_fvector( Deltac, 0, n-1 );
}

void read_column_ids( id )
int id;
{
    TCCSER( id, ":X",       &Col.x );
    TCCSER( id, ":Y",       &Col.y );
    TCCSER( id, ":PEAK",    &Col.peak );
    TCCSER( id, ":IDENT",   &Col.ident );
    TCCSER( id, ":WAVE",    &Col.wave );
    TCCSER( id, ":WAVEC",   &Col.wavec );
    TCCSER( id, ":DELTA",   &Col.delta );
    TCCSER( id, ":DELTAC",  &Col.deltac );
    TCCSER( id, ":ERASED",  &Col.erased );

    if ( Col.x == -1 || Col.y == -1 || Col.peak == -1 || Col.ident == -1 ||
         Col.wave == -1 || Col.wavec == -1 || 
         Col.delta == -1 || Col.deltac == -1 ) {
	SCTPUT( "*** Starting line has not been calibrated ***" );
	end_midas();
    }

    if ( Col.erased == -1 )
	TCCINI( id, D_C_FORMAT, 1, "A1", " ", "ERASED", &Col.erased );
}


void init_midas()
{
    SCSPRO( "SPPLOT" );
    TCMNUL( &Inull, &Rnull, &Dnull ); /* obtain NULL values */
}

void end_midas()
{
    SCSEPI();
}

void read_parameters()
{
    int actval; 	/* actual values returned */
    int  unit;		/* useless */
    int nulval;		/* useless */
    int pltkey[2];

    SCKRDI( "SPPLT", 1, 2, &actval, pltkey, &unit, &nulval );
    SCKRDI( "DCX", 1, 1, &actval, &Fitd, &unit, &nulval );
    SCKRDI( "YSTART", 1, 1, &actval, &Ystart, &unit, &nulval );
    SCKRDR( "IMIN", 1, 1, &actval, &Imin, &unit, &nulval );
    SCKRDI( "WRANG", 1, 2, &actval, Wrang, &unit, &nulval );
    SCKGETC( "WLC", 1, 20, &actval, Wlc );
    SCKGETC( "LINTAB", 1, 20, &actval, Lintab );
    SCKGETC( "LINCAT", 1, 20, &actval, Lincat );
    SCKGETC( "COERBR", 1, 20, &actval, Coerbr );
    SCKGETC( "INPUTC", 1, 20, &actval, PlotFile );
    PlotAction = pltkey[0];
    PlotType   = pltkey[1];
}

/******************************************************
 * save_limits() : Save current viewport coordinates.
 ******************************************************/
void save_limits( xmin, xmax, ymin, ymax )
float xmin, xmax, ymin, ymax;
{
    int i;
    int  unit;
    float lims[4];

    lims[0] = xmin;
    lims[1] = xmax;
    lims[2] = ymin;
    lims[3] = ymax;
    SCKWRR( "AGLIMS", lims, 1, 4, &unit );
}

/********************************************************
 * read_limits() : Read last saved viewport coordinates.
 ********************************************************/
void read_limits( xmin, xmax, ymin, ymax )
float *xmin, *xmax, *ymin, *ymax;
{
    int i;
    int  unit;
    float lims[4];
    int  actval, nulval;	/* garbage */

    SCKRDR( "AGLIMS", 1, 4, &actval, lims, &unit, &nulval );
    *xmin = lims[0];
    *xmax = lims[1];
    *ymin = lims[2];
    *ymax = lims[3];
}

/*************************************************************
 * init_viewport() : Init viewport according to the last one.
 *************************************************************/
void init_viewport()
{
    float xmin, xmax, ymin, ymax;

    read_limits( &xmin, &xmax, &ymin, &ymax ); 
    AG_AXES( xmin, xmax, ymin, ymax, " " );
    AG_VUPD();
}


#define FMT_TITLE	"File: %s  Line: %d  Image: %s"
#define FMT_OPTIONS	"TITLE=%s;LABX=Position;LABY=Pixel value"

#define ORDMAX		9999999.0
#define ORDMIN		-9999999.0

/********************************************************
 * read_image(): plots the 'row'nth line of 'image'.
 ********************************************************/
void read_image( row, image )
int row;
char *image;
{
    int      i, j, framid;
    int ncols, nrows, npix[2];
    double start, step;
    float cuts[4];
    int  unit;
    int  nulval, retval; /* garbage */
    char options[512], title[512];
    float x[MAXDATA], y[MAXDATA];
    char ident[MAXIDENT+1];
    float xmin, xmax, ymin, ymax;

    SCFOPN( image, D_R4_FORMAT, 0, F_IMA_TYPE, &framid );
    SCDRDI( framid, "NPIX", 1, 2, &retval, npix, &unit, &nulval);
    ncols = npix[0];
    nrows = npix[1];

    SCDRDR( framid, "LHCUTS", 1, 4, &retval, cuts, &unit, &nulval );
    SCDRDD( framid, "START", 1, 1, &retval, &start, &unit, &nulval );
    SCDRDD( framid, "STEP", 1, 1, &retval,  &step, &unit, &nulval );
    SCDGETC( framid, "IDENT", 1, MAXIDENT, &retval, ident );

    /* mapping 'row' row (Y axis) */
    SCFGET( framid, (row-1)*ncols+1, ncols, &retval, (char *)y );

    for ( i = 0; i < ncols; i++ )
	x[i] = start + step*i;

    xmin = start;
    xmax = start + step*(ncols - 1);

    if ( cuts[1] != 0.0 ) {
    	ymin = cuts[0]; 
    	ymax = cuts[1];
    }
    else if ( cuts[3] != 0.0 ) {
    	ymin = cuts[2]; 
    	ymax = cuts[3];
    }
    else {	/* Max. and min. cuts are not assigned */
    	ymin = ORDMAX;
    	ymax = ORDMIN;
    	for ( i = 0; i < ncols; i++ ) {
	    if ( y[i] > ymax )
	    	ymax = y[i];
	    if ( y[i] < ymin )
	    	ymin = y[i];
    	}
        cuts[0] = ymin;
        cuts[1] = ymax;
    	SCDWRR( framid, "LHCUTS", cuts, 3, 2, &unit );
    }
    SCFCLO(framid);

    sprintf( title, FMT_TITLE, image, row, ident );
    sprintf( options, FMT_OPTIONS, title );

    AG_VERS();
    AG_AXES( xmin, xmax, ymin, ymax, options );

    AG_GPLL( x, y, ncols );
    AG_VUPD();

    save_limits( xmin, xmax, ymin, ymax ); 
}

void init_graphic( devtype )
int devtype;
{
    if ( !graphwin_exists() ) {
	SCTPUT( "*** Please create the graphic window ***" );
	end_midas();
    }

    if ( InitGraphic ) 
	return;

    InitGraphic = TRUE;
    get_agldev();
    switch ( devtype ) {
	case DEV_ERASE:
    		AG_VDEF( DevErase, 0.05, 1.0, 0.0, 1.0, 0.0, 0.0 );
		break;
	case DEV_NO_ERASE:
    		AG_VDEF( DevNoErase, 0.05, 1.0, 0.0, 1.0, 0.0, 0.0 );
		break;
    }
}

void end_graphic()
{
    if ( InitGraphic && graphwin_exists() )
    	AG_CLS();

    InitGraphic = FALSE;
}

/*********************************************************
 * get_agldev(): translate IDI device to devices erasable
 * and non-erasable suitable for AG_VDEF calls.
 *********************************************************/
void get_agldev()
{
    int actval; 	/* actual values returned */
    int  unit;		/* useless */
    int nulval;		/* useless */
    char devkeyw[21];	/* type of device as in MID$PLOT */
    char device[21];	/* name of assoc. device as in pltdevices.dat */

    /* read & translate type of device to a specific device name */
    /* SCKGETC( "MID$PLOT", 1, 20, &actval, devkeyw );
       get_dev( devkeyw, device ); Not by now... */
    strcpy( device, "GRAPH_WND0" );

    /* now make the AGL device names */
    strcpy( DevErase, device ); 
    strcat( DevErase, ":" ); /* see AG_VDEF definition */

    strcpy( DevNoErase, device ); 
    strcat( DevNoErase, "/n:" ); /* see AG_VDEF definition */
}

void dpoly( x, p, np )  /* (p[i] <-- x**(i-1)), i = 1,np) */
double x, p[];
int np;
{
    int j;

    p[1] = 1.0;
    for ( j = 2; j <= np; j++ )
	p[j] = p[j-1] * x;
}

void read_rebin_parameters()
{
    int diff, min_diff = 32767;  /* :-) */
    int i, id;
    int col_y, col_rms, col_pixel;
    int  actval, nulval;        /* useless */
    double y, pixel, rms;
    int  sortcol, aw, ar, ncols, nrows;
    int  unit;
    

    if ( ! file_exists(Coerbr, ".tbl") ) {
        SCTPUT( "Coefficients table couldn't be opened. Stop." );
        end_midas();
    }
    TCTOPN(Coerbr, F_IO_MODE, &id);
    TCIGET( id, &ncols, &nrows, &sortcol, &aw, &ar );

    if ( nrows == 0 ) {
        SCTPUT( "Error: coefficients table is empty." );
        end_midas();
    }

    TCCSER( id, ":ROW", &col_y );
    TCCSER( id, ":RMS", &col_rms );
    TCCSER( id, ":PIXEL", &col_pixel );
    if ( col_y == -1 || col_rms == -1 || col_pixel == -1 ) {
	SCTPUT( "Calibration process has not been performed. Stop." );
	end_midas();
    }

    for ( i = 1; i <= nrows; i++ ) {
	TCERDD( id, i, col_y, &y, &nulval );
	TCERDD( id, i, col_pixel, &pixel, &nulval );
	TCERDD( id, i, col_rms, &rms, &nulval );
	diff = fabs(Ystart - y);
	if ( diff < min_diff ) {
	    min_diff = diff;
	    Ycoerbr = i;
	    Rms = rms;
	    Dispersion = pixel;
	}
    }
    TCTCLO( id );
}

void read_coefs_ystart()
{
    int  nulval;        /* useless */
    double value;
    int i, id;

    TCTOPN(Coerbr, F_IO_MODE, &id);
    for ( i = 3; i <= Fitd + 3; i++ ) {
	TCERDD( id, Ycoerbr, i, &value, &nulval );
	Coef[i-3] = value;
    }
    TCTCLO(id);
}

void compute_calib_values()
{
    double *a; /* coefficients for fit = 1 */
    double val;
    double *xfit, *yfit;
    double delta, delta_min;
    int id, i, j, k, n = 1;

    a = dvector( 1, Fitd+1 );
    xfit = dvector( 1, NumLine );
    yfit = dvector( 1, NumLine );

    for ( k = 0; k < NumDel; k++ ) {
	i = DelList[k];
        Wavec[i] = eval_dpoly( X[i], Coef-1, Fitd+1 );
    	delta_min = MAXWAVE;
        for ( j = 0; j < Lc->nrows; j++ )
            if ( (delta = fabs(Lc->wave[j] - Wavec[i])) < delta_min &&
		 ! Lc->sel[j] ) {
                Wave[i] = Lc->wave[j];
                delta_min = delta;
            }
    }

    for ( j = 0; j < NumLine; j++ )
	if ( Wave[j] != Rnull ) {
	    xfit[n] = X[j];
	    yfit[n] = Wave[j];
	    n++;
    	}

    /* get the coefficients of the relation w=a0+a1x,
       calculate Delta and Deltac */
    lfit( xfit, yfit, n - 1, a, 2, dpoly ); /* fitting */
    for ( k = 0; k < NumDel; k++ ) {
	i = DelList[k];
        val = a[1] + a[2] * X[i];
        Deltac[i] = val - Wavec[i];
        Delta[i] = val - Wave[i];
    }
    
    /* now save the wave values in the table (for lnerase.exe) */
    TCTOPN( Lintab, F_IO_MODE, &id );
    for ( i = 0; i < NumDel; i++ )
        TCEWRR( id, Row[DelList[i]], Col.wave, &Wave[DelList[i]] );
    TCTCLO(id);

    free_dvector( xfit, 1, NumLine );
    free_dvector( yfit, 1, NumLine );
}

/*****************************************************************
 * read_lincat_table(): reads line catalog named by 'Lincat' and
 * stores the values in the 'Lc' structure, under the
 * restrictions Wrang & Imin.
 *****************************************************************/
int read_lincat_table()
{
    if ( ! file_exists( Lincat, ".tbl" ) ) {
        SCTPUT( "*** Line catalogue doesn't exist ***" );
        return(FALSE);
    }

    if ( Lc != NULL )
	free_catalog_table(Lc);
    Lc  = (LCTAB *)osmmget( sizeof(LCTAB) );

    /* call to a library routine */
    if( !read_catalog_table( Lc, Lincat, Wrang, Imin ) ) {
        Lc = NULL;
        return FALSE;
    }
    return TRUE;
}