/* @(#)spgrap.c	17.1.1.1 (ES0-DMD) 01/25/02 17:29:51 */
/*===========================================================================
  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
===========================================================================*/

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.IDENT        amain.c
.MODULE       subroutines -- gmidas.exe
.LANGUAGE     C
.AUTHOR       Cristian Levin - ESO La Silla
.PURPOSE      The functions of this module implement the
              applications of the Alice package and the 
              spectrum manipulation functions.
.KEYWORDS     spectrum plotting, line integration.
.VERSION 1.0  1-Mar-1991   Implementation
.ENVIRONMENT  UNIX
------------------------------------------------------------*/

#include <stdio.h>

#include <gl_defs.h> 	/* local headers */
#include <agl_defs.h> 	/* local headers */
#include <main_defs.h>
#include <spec_defs.h>

#include <fileexts.h> 	/* midas headers */

#include <agl.h> 	/* AGL headers */
#include <aglerror.h>

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

#define ORDMAX		9999999.0
#define ORDMIN		-9999999.0

#define NEXT_SPEC	0
#define PREV_SPEC	1
#define FORM_SPEC	2

#define MIN_DIF_CUTS	0.01

/* Definitions for printing a plot */
#define CURRFILE	"current.plt"

void init_gmidas();
void end_gmidas();
void end_graphic();
void get_agldev();

int  Areadim();
void Agetcur();
int  Acutx();
int  Acuty();
void Aunzoom();
void Ahelp();

void open_plotfile();
void close_plotfile();

void spec_copy();
void get_image_name();

void redraw_spectrum();

void define_viewport();
void set_viewport();
static float vwx_min, vwy_min, vwx_max, vwy_max;
static float cwx_min, cwy_min, cwx_max, cwy_max;

float *fvector(), **fmatrix();
double *dvector(), **dmatrix();
int *ivector();
double eval_dpoly();
double exp(), sqrt();

static SPEC *	Sp[MAXSTACK]; 			/* stack of spectra */
static SPEC *	Spcur;				/* current spectrum */
static int 	Top = 0;    			/* index to top of stack */
static char 	DevErase[80], DevNoErase[80];	/* devices for AG_VDEF */
static int	InitGraphic = FALSE;
static int      XLimDefined = FALSE;
static double   XLimits[2];

/******************************************************
 * init_gmidas(): initialize Midas & Alice structures.
 ******************************************************/
void init_gmidas( progname )
char *progname;
{
    int i;

    SCSPRO( progname );

    for ( i = 0; i < MAXSTACK; i++ ) /* initialize stack */
    	Sp[i] = (SPEC *)osmmget( sizeof(SPEC) );
    Spcur = (SPEC *)osmmget( sizeof(SPEC) );

    TCMNUL( &Inull, &Rnull, &Dnull ); /* NULL values */
}

int exists_graphic()
{
    if ( !graphwin_exists() ) {
SCTPUT( "*** Please create the graphic window ***" );
	return(FALSE);
    }
    return TRUE;
}

int init_graphic( devtype )
int devtype;
{
    if ( InitGraphic )
	return(TRUE);

    if ( !graphwin_exists() ) {
	SCTPUT( "*** Please create the graphic window ***" );
	return(FALSE);
    }

    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;
    }
    AG_SSET("FONT=1");

    return(TRUE);
}

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

    InitGraphic = FALSE;
}

void open_plotfile()
{
    AG_MOPN( CURRFILE );
}

void close_plotfile()
{
    AG_MCLS();
}

void clear_graphic()
{
    if ( InitGraphic && graphwin_exists() ) {
        AG_VERS();
        AG_VUPD();
    }
}

/******************************************************
 * end_gmidas(): end Midas & Alice. 
 ******************************************************/
void end_gmidas()
{
    int i;

    for ( i = 0; i < MAXSTACK; i++ ) /* deallocate stack */
    	osmmfree(Sp[i]);
    osmmfree(Spcur);

    SCSEPI();
}

/********************************************************************
 * Areadim(): reads the 'row'nth line plus 'ahead' lines ahead
 *            of 'image_name' and put the average onto the stack. 
 ********************************************************************/
int Areadim( image_name, row, ahead, save_lims )
char *image_name;
int row;	/* row number to read */
int ahead;	/* ahead rows to read */
int save_lims;	/* TRUE if should save limits */
{
    int i, j;
    char msg[30];
    char ident[MAXIDENT+1];
    int ncols, nrows, npix[2], naxis;
    double start, step;
    float cuts[4];
    float *group, sum;
    int framid;
    int status;
    int unit;
    int nulval, retval; /* garbage */
    double fabs();

    if ( ! file_exists(image_name, ".bdf" ) )
	return(FALSE);

    SCFOPN( image_name, D_R4_FORMAT, 0, F_IMA_TYPE, &framid );

    SCDRDI( framid, "NPIX", 1, 2, &retval, npix, &unit, &nulval );
    SCDRDI( framid, "NAXIS", 1, 1, &retval, &naxis, &unit, &nulval );
    ncols = npix[0];
    nrows = (naxis == 1) ? 1 : npix[1];

    if ( row + ahead > nrows || row + ahead < 1 ) {
	row = 1;
	ahead = 0;
    }

    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 and averaging of "row + ahead" rows */
    group = fvector( 0, (ahead+1)*ncols - 1);
    SCFGET( framid, (row-1)*ncols+1, (ahead+1)*ncols, &retval, (char *)group );

    for ( i = 0; i < ncols; i++ ) {
	sum = 0;
    	for ( j = 0; j <= ahead; j++ )
	    sum += group[i+j*ncols]; 
	Sp[Top]->y[i] = sum/(ahead+1);
    }
    free_fvector(group, 0, (ahead+1)*ncols - 1 );

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


    if (XLimDefined == FALSE) {
      Sp[Top]->xmin = start;
      Sp[Top]->xmax = start + step*(ncols - 1);
      XLimits[0]    = start;
      XLimits[1]    = start + step*(ncols - 1);
    }
    else {
      Sp[Top]->xmin = Spcur->xmin;
      Sp[Top]->xmax = Spcur->xmax;
    }

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

    strcpy( Sp[Top]->name, image_name );
    strcpy( Sp[Top]->ident, ident );
    get_image_name( Sp[Top]->name );

    Sp[Top]->ncols = ncols;
    Sp[Top]->nrows = nrows;
    Sp[Top]->start = start;
    Sp[Top]->step  = step;
    Sp[Top]->row   = row;
    Sp[Top]->ahead = ahead;

    spec_copy( Spcur, Sp[Top] );

    init_graphic(DEV_ERASE);
    redraw_spectrum();
    end_graphic();

    return(TRUE);
}

void Alabelx( xvals, yvals, num_xvals, color )
float *xvals, *yvals;
int num_xvals;
int color;
{
    float zeros[MAXDATA];
    float wlimits[4];		/* window limits */
    float x[2], y[2];
    int i;

    init_graphic(DEV_NO_ERASE);
    set_viewport();

    AG_SSET( LARGE_SYM );
    AG_SSET( (char *)color );
    AG_RGET( "WNDL", wlimits );

    for ( i = 0; i < num_xvals; i++ )
	zeros[i] = wlimits[YMIN];

    /* first, plot a small tick, enough to see the identified line */
    AG_GPLM( xvals, zeros, num_xvals, VBAR_MARKER );

    /* now, plot a line as high as the intensity value is, for each
       identified line */
    for ( i = 0; i < num_xvals; i++ ) {
	x[0] = x[1] = xvals[i];
	y[0] = zeros[i];
	y[1] = yvals[i];
        AG_GPLL( x, y, 2 );
    }
    AG_SSET( DEF_SYM );
    AG_SSET( DEF_COLOR );
 
    end_graphic();
}

void Agetcur()
{
    float cpx, cpy;
    int key, valpix;
    char msg[30];

    init_graphic(DEV_NO_ERASE);
    set_viewport();

    cpx = Spcur->xmin;
    cpy = Spcur->ymin;

    SCTPUT(" ");
    SCTPUT("   X-axis         Y-axis");
    SCTPUT("--------------------------");
    while ( 1 ) { /* forever */
    	AG_VLOC( &cpx, &cpy, &key, &valpix );
	if ( key == MID_BUT ) 
	    break;
	sprintf( msg, "%10.2f    %10.2f", cpx, cpy );
    	SCTPUT( msg );
    }
    end_graphic();
}

int Acutx()
{
    float cpx[2], cpy[2];
    int i, key, valpix;

    init_graphic(DEV_NO_ERASE);
    set_viewport();

    cpx[0] = Spcur->xmin;
    cpy[0] = Spcur->ymin;

    /* get cursor position */
    for ( i = 0; i < 2; i++ ) {
    	AG_VLOC( cpx+i, cpy+i, &key, &valpix );
	if ( key == MID_BUT ) {
            end_graphic();
	    return( MID_BUT );
	}
    }

    if ( cpx[0] < cpx[1] ) {
    	Spcur->xmin = cpx[0]; 
	Spcur->xmax = cpx[1];
    }
    else {
    	Spcur->xmin = cpx[1]; 
	Spcur->xmax = cpx[0];
    }
    redraw_spectrum();
    end_graphic();

    XLimDefined = TRUE;

    return( key );
}

int Acuty()
{
    float cpx[2], cpy[2], cuts[2];
    int i, key, valpix, id;
    int unit;

    init_graphic(DEV_NO_ERASE);
    set_viewport();

    cpx[0] = Spcur->xmin;
    cpy[0] = Spcur->ymin;

    /* get cursor position */
    for ( i = 0; i < 2; i++ ) {
    	AG_VLOC( cpx+i, cpy+i, &key, &valpix );
	if ( key == MID_BUT ) {
            end_graphic();
	    return( MID_BUT );
	}
    }

    if ( cpy[0] < cpy[1] ) {
    	Spcur->ymin = cpy[0]; 
	Spcur->ymax = cpy[1];
    }
    else {
    	Spcur->ymin = cpy[1]; 
	Spcur->ymax = cpy[0];
    }
    redraw_spectrum();
    end_graphic();

    cuts[0] = Spcur->ymin;
    cuts[1] = Spcur->ymax;
    SCFOPN( Spcur->name, D_R4_FORMAT, 0, F_IMA_TYPE, &id );
    SCDWRR( id, "LHCUTS", cuts, 1, 2, &unit );
    SCFCLO(id);

    return( key );
}

int Ashift()
{
    float cpx, cpy;
    int key, valpix, middle;

    init_graphic(DEV_NO_ERASE);
    set_viewport();

    cpx = Spcur->xmin;
    cpy = Spcur->ymin;

    /* get cursor position */
    AG_VLOC( &cpx, &cpy, &key, &valpix );
    if ( key == MID_BUT ) {
        end_graphic();
	return( MID_BUT );
    }

    middle = (Spcur->xmax - Spcur->xmin) / 2.0;
    if ( cpx + middle > Sp[Top]->xmax ) {
    	Spcur->xmin = Spcur->xmin + Sp[Top]->xmax - Spcur->xmax;
	Spcur->xmax = Sp[Top]->xmax;
    }
    else if ( cpx - middle < Sp[Top]->xmin ) {
	Spcur->xmin = Sp[Top]->xmin;
    	Spcur->xmax = Spcur->xmax + Sp[Top]->xmin - Spcur->xmin;
    }
    else {
	Spcur->xmin = cpx - middle;
    	Spcur->xmax = cpx + middle;
    }
    redraw_spectrum();
    end_graphic();

    return( key );
}

void Aunzoom()
{
    int unit;
    int retval, nulval, id;
    float cuts[2];

    init_graphic(DEV_NO_ERASE);
    set_viewport();

    spec_copy( Spcur, Sp[Top] );

    SCFOPN( Spcur->name, D_R4_FORMAT, 0, F_IMA_TYPE, &id );
    SCDRDR( id, "LHCUTS", 3, 2, &retval, cuts, &unit, &nulval );
    SCDWRR( id, "LHCUTS", cuts, 1, 2, &unit );
    SCFCLO(id);

    Spcur->xmin = XLimits[0];
    Spcur->xmax = XLimits[1];

    Spcur->ymin = cuts[0];
    Spcur->ymax = cuts[1];

    redraw_spectrum();
    end_graphic();
}

void Ahelp()
{
}

/******************************************************
 * spec_copy()
 ******************************************************/
void spec_copy( sp_to, sp_from )
SPEC *sp_to, *sp_from;
{
    int i;

    sp_to->ncols = sp_from->ncols;
    sp_to->nrows = sp_from->nrows;
    sp_to->start = sp_from->start;
    sp_to->step  = sp_from->step;
 
    sp_to->xmin = sp_from->xmin;
    sp_to->xmax = sp_from->xmax;
    sp_to->ymin = sp_from->ymin;
    sp_to->ymax = sp_from->ymax;

    sp_to->row   = sp_from->row;
    sp_to->ahead = sp_from->ahead;

    for ( i = 0; i < sp_from->ncols; i++ ) {
	sp_to->x[i] = sp_from->x[i];
	sp_to->y[i] = sp_from->y[i];
    }
    strcpy( sp_to->name, sp_from->name ); 
    strcpy( sp_to->ident, sp_from->ident ); 
}

/*************************************************************
 * get_image_name()
 *************************************************************/
void get_image_name( c )
char c[];
{
    int i;

    /* strip root of the filename */
    for ( i = strlen(c); i > 0; i-- )
	if ( c[i] == '/' ) {
	    strcpy( c, c+i+1 );
	    break;
	}

    /* strip suffix */
    for ( i = strlen(c); i > 0; i-- )
	if ( c[i] == '.' ) {
	    c[i] = '\0';
	    return;
	}	
}

/*********************************************************
 * get_agldev(): translate IDI device to devices erasable
 * and non-erasable suitable for AG_VDEF calls.
 *********************************************************/
void get_agldev()
{
    int actval; 	/* actual values returned */
    char *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 define_viewport()
{
    float vals[4];

    AG_RGET("WNDL", vals);
    vwx_min = vals[0];
    vwx_max = vals[1];
    vwy_min = vals[2];
    vwy_max = vals[3];
    AG_RGET("CLPL", vals);
    cwx_min = vals[0];
    cwx_max = vals[1];
    cwy_min = vals[2];
    cwy_max = vals[3];
}

void set_viewport()
{
    AG_CDEF(cwx_min, cwx_max, cwy_min, cwy_max);
    AG_WDEF(vwx_min, vwx_max, vwy_min, vwy_max);
}

/***************************************************************
 * redraw_spectrum(): redraw the current spectrum of the stack.
 ***************************************************************/
void redraw_spectrum()
{
    char options[512], title[512], interval[512];
    char small_ident[11];

    if ( Spcur->ahead == 0 )
        sprintf( interval, "%d", Spcur->row );
    else
        sprintf( interval, "%1d-%1d", Spcur->row, Spcur->row + Spcur->ahead );

    strncpy(small_ident, Spcur->ident, 10);
    small_ident[10] = '\0';
    sprintf( title, FMT_TITLE, Spcur->name, interval, small_ident );
    sprintf( options, FMT_OPTIONS, title );

    clear_graphic();
    open_plotfile();

    AG_AXES( Spcur->xmin, Spcur->xmax, Spcur->ymin, Spcur->ymax, options );
    define_viewport();

    AG_GPLL( Spcur->x, Spcur->y, Spcur->ncols );

    close_plotfile();
}