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

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.COPYRIGHT   (c) 1993 European Southern Observatory
.IDENTifer   module PLCON
.AUTHOR      R.M. van Hees IPG-ESO Garching
.KEYWORDS    plot software, graphics, bulk data fame, contour plotting
.LANGUAGE    C
.PURPOSE     Contour plotting routine for two dimensional data
.COMMENTS    holds PLCON and PLCONI
.ENVIRONment MIDAS and AGL
             #include <agl.h>           Prototypes for AGL application programs
             #include <midas_def.h>     Prototypes for MIDAS interfaces
             #include <plot_def.h>      Symbols used by the PLT interfaces

.VERSION     1.3     14-Jul-1994   add option for colors RvH
             1.2     25-Nov-1993   rewritten             RvH
             1.1     10-Sep-1993   FORTRAN --> C         RvH
------------------------------------------------------------*/
/*
 * Define _POSIX_SOURCE to indicate
 * that this is a POSIX program
 */
#define  _POSIX_SOURCE 1

/*
 * definition of the used functions in this module
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <agl.h>
#include <midas_def.h>

/*
 * define some macros and constants
 */
#include <plot_def.h>

#define NEXT_S		1
#define DRAW_C		0

/*
 * A segment from a isophote is defined by:
 */
struct segmnt { int   row;                /* row in which the point is found */
                int   use;                      /* is it already used or not */
                double x1;                         /* x coord of first point */
                double y1;                         /* y coord of first point */
                double x2;                        /* x coord of second point */
                double y2;                        /* y coord of second point */
              };

/*++++++++++++++++++++
.IDENTifer   COLLECT_SEGMNT
.PURPOSE     collect segments which together form the isophotes
      input:   double *start        start value of the first pixel in w.c.
               double *step         distance between pixels in w.c.
               int    *npix         size of the frame in pixels
               float  *p_img        data of the frame
               float  levl          value of the isophote
     output:   int    nrdata        number of structs point
               struct segmnt *data  pointer first struct point
.COMMENTS    static function
--------------------*/
#ifdef __STDC__
      static void COLLECT_SEGMNT( double *start, double *step, int *npix,
                                  float *p_img, float levl, 
                                  int *nrdata, struct segmnt *data )
#else
      static void COLLECT_SEGMNT( start, step, npix, p_img, levl, 
                                                     nrdata, data )
      int    *npix, *nrdata;
      float  levl, *p_img;
      double *start, *step; 
      struct segmnt *data;
#endif
{
register int ix, iy;

int     icase;
double  v1, v2, v3, v4, x, y;

struct  segmnt *first;

first = data;

y = start[1];
for ( iy = 0; iy < npix[1]-1; iy++ )
    { x = start[0];
      for ( ix = 0; ix < npix[0]-1; ix++ )
          { v1 = *p_img;
            v3 = *(p_img + npix[0]);
            v2 = *++p_img;
            v4 = *(p_img + npix[0]);

            icase = 0;
            if ( levl > v1 ) icase += 1;
            if ( levl > v2 ) icase += 2;
            if ( levl > v3 ) icase += 4;
            if ( levl > v4 ) icase = 7 - icase;

            if ( icase > 0 )
               { if ( icase != 7 )
                    { if ( icase == 3 || icase == 4 )
                         { data->x1  = x;
                           data->y1  = y + step[1] * (levl-v1)/(v3-v1);;
                         }
                      else 
                         { data->x1  = x + step[0] * (levl-v1)/(v2-v1);
                           data->y1  = y;
                         }

                      if ( icase == 1 || icase == 6 )
                         { data->x2  = x;
                           data->y2  = y + step[1] * (levl-v1)/(v3-v1);
                           if ( icase == 6 && (levl != v2 || levl != v3) )
                              { data->row = iy + 1;
                                data++->use = FALSE;
                                data->x1 = x + step[0] * (levl-v3)/(v4-v3);
                                data->y1 = y + step[1];
                                data->x2 = x + step[0];
                                data->y2 = y + step[1] * (levl-v2)/(v4-v2);
                              }
                         }
                      else if ( icase == 2 || icase == 3 )
                         { data->x2  = x + step[0];
                           data->y2  = y + step[1] * (levl-v2)/(v4-v2);
                         }
                      else 
                         { data->x2  = x + step[0] * (levl-v3)/(v4-v3);
                           data->y2  = y + step[1];
                         }
                      data->row = iy + 1;
                      data++->use = FALSE;
                    }
                 else
                    { data->x1 = x + step[0] * (levl-v3)/(v4-v3);
                      data->y1 = y + step[1];
                      data->x2 = x + step[0];
                      data->y2 = y + step[1] * (levl-v2)/(v4-v2);
                      data->row = iy + 1;
                      data++->use = FALSE;
                    }               
               }
            x += step[0];
          }
      y += step[1];
      p_img++;
    }
*nrdata = data - first;
return;
}

/*++++++++++++++++++++
.IDENTifer   JOIN_SEGMNT
.PURPOSE     Draw isophotes given a set of points
      input:   int    nrdata        number of structs point
               struct segmnt *data  pointer first struct point
.COMMENTS    static function
--------------------*/
#ifdef __STDC__
      static void JOIN_SEGMNT( int nrdata, struct segmnt *data )
#else
      static void JOIN_SEGMNT( nrdata, data)
      int    nrdata;
      struct segmnt *data;
#endif
{
register int ii, indx;

int    bgn, line, nr, stat, str_indx;
float  *xdata, *ydata;
double *xx, *yy;

/*
 * if no points are found --> return
 */
if ( nrdata < 1 ) return;

xx = (double *) osmmget( (nrdata+1) * sizeof( double ));
yy = (double *) osmmget( (nrdata+1) * sizeof( double ));

/*
 * initialisation
 */
nr = 1;
indx = str_indx = 0;
line  = data->row;
*xx   = data->x1;
*yy   = data->y1;
*++xx = data->x2;
*++yy = data->y2;
data->use = TRUE;
stat = NEXT_S;

while ( indx < nrdata )
/*
 * search for next segment
 */
   { bgn = indx;
     while ( ++indx < nrdata && (data+indx)->row <= line + 1 )
        { if ( (data+indx)->use == FALSE &&
               ( (*xx == (data+indx)->x1 && *yy == (data+indx)->y1)
              || (*xx == (data+indx)->x2 && *yy == (data+indx)->y2)
               )
             )
             break;
        }
     if ( indx == nrdata || (data+indx)->row > line + 1 )
        { indx = bgn;
          while ( --indx >= 0 && (data+indx)->row >= line - 1 )
             { if ( (data+indx)->use == FALSE &&
                    ( (*xx == (data+indx)->x1 && *yy == (data+indx)->y1)
                   || (*xx == (data+indx)->x2 && *yy == (data+indx)->y2)
                    )
                  )
                  break;
             }
          if ( indx < 0 || (data+indx)->row < line - 1 ) stat = DRAW_C;
        }
/*
 * do we have to draw a contour?
 */
     if ( stat == DRAW_C )
        { xx -= nr;
          yy -= nr;
          if ( nr > 1 )
             { xdata = (float *) osmmget( ++nr * sizeof( float ));
               ydata = (float *) osmmget(   nr * sizeof( float ));
               for ( ii = 0; ii < nr; ii++ )
                   { xdata[ii] = (float) xx[ii];
                     ydata[ii] = (float) yy[ii];
                   }
               AG_GPLL( xdata, ydata, nr );
               (void) osmmfree( (char *) xdata );
               (void) osmmfree( (char *) ydata );
              }
/*
 * go to next starting point
 */
          indx = str_indx;
          while ( indx < nrdata && (data+indx)->use == TRUE ) indx++;
          str_indx = indx;
          line  = (data+indx)->row;
          *xx   = (data+indx)->x1;
          *yy   = (data+indx)->y1;
          *++xx = (data+indx)->x2;
          *++yy = (data+indx)->y2;
          (data+indx)->use = TRUE;
          nr = 1;
          stat = NEXT_S;
        }
     else
        { nr++;
          if ( *xx == (data+indx)->x1 && *yy == (data+indx)->y1 )
             { *++xx = (data+indx)->x2;
               *++yy = (data+indx)->y2;
             }
          else
             { *++xx = (data+indx)->x1;
               *++yy = (data+indx)->y1;
             }
          line = (data+indx)->row;
          (data+indx)->use = TRUE;
        }
   }
xx -= nr;
yy -= nr;
if ( nr > 1 ) 
   { xdata = (float *) osmmget( ++nr * sizeof( float ));
     ydata = (float *) osmmget(   nr * sizeof( float ));
     for ( ii = 0; ii < nr; ii++ )
         { xdata[ii] = (float) xx[ii];
           ydata[ii] = (float) yy[ii];
         }
     AG_GPLL( xdata, ydata, nr );
     (void) osmmfree( (char *) xdata );
     (void) osmmfree( (char *) ydata );
   }
(void) osmmfree( (char *) xx );
(void) osmmfree( (char *) yy );
return;
}

/*++++++++++++++++++++++++++++++
.IDENTifer   PLCON
.PURPOSE     Contour plotting routine for two dimensional data
     input:   float  *p_img   pointer to the data of the frame
              float  *image   image size in pixel coordinates
              float  *area    image size in world coordinates
              double *step    distance between pixels in world units
              int    nlevl    number of contours
              float  *clevl   level of the contours
              int    *ctype   ctype >= 0: array with line-types of contours
                              ctype < 0 : array with colors of contours

.COMMENTS    if image[0] > image[1] or image[2] > image[3] then
             the values are put in the right order
--------------------------------*/
void PLCON( p_img, image, area, step, nlevl, clevl, ctype )
int    nlevl, *ctype;
float  *p_img, *image, *area, *clevl;
double *step;

{
register int ic;

char   text[21];
int    actvals, color, ctagl, ltype, ltagl, nrdata, npix[2];
double start[2];

struct segmnt *data;

/*
 * return if no contours have to be drawn
 */
if ( nlevl == 0 ) return;

/*
 * find start values in x and y
 */
if ( step[0] < 0.0 )
   start[0] = MYMAX( area[0], area[1] );
else
   start[0] = MYMIN( area[0], area[1] );

if ( step[1] < 0.0 )
   start[1] = MYMAX( area[2], area[3] );
else
   start[1] = MYMIN( area[2], area[3] );

/*
 * get line type from PLISTAT
 */
PCKRDI( "LTYPE", 1, &actvals, &ltype );
PCKRDI( "COLOR", 1, &actvals, &color );
ltagl = ltype = MYMAX( ltype-1, 0 );
ctagl = -color;
(void) sprintf( text, "lstyl=%1d", ltype );
AG_SSET( text );

/*
 * get the points along the isophotes
 */
npix[0] = (int) fabs( (double) image[1] - image[0] ) + 1;
npix[1] = (int) fabs( (double) image[3] - image[2] ) + 1; 
data = (struct segmnt *) osmmget( npix[0]* npix[1]* sizeof( struct segmnt));

/*
 * collect the data isophote by isophote
 */
for ( ic = 0; ic < nlevl; ic++ )
    { if ( ctype[ic] != ltagl && ctype[ic] != ctagl )
         { if ( ctype[ic] < 0 )                              /* change color */
              { ctagl = ctype[ic];
                (void) sprintf( text, "color=%1d", -ctagl );
              }
           else                                              /* change ltype */
              { ltagl = ctype[ic];
                (void) sprintf( text, "lstyl=%1d", ltagl );
              }
           AG_SSET( text );
         }
/*
 * collect the segments
 */
      COLLECT_SEGMNT( start, step, npix, p_img, clevl[ic], &nrdata, data );
/*
 * join the segments and draw them in the current viewport
 */
      JOIN_SEGMNT( nrdata, data );
    }
/*
 * reset line style
 */
if ( ltype != ltagl ) 
   { (void) sprintf( text, "lstyl=%1d", ltype );
     AG_SSET( text );
   }
else if ( color != -ctagl ) 
   { (void) sprintf( text, "color=%1d", color );
     AG_SSET( text );
   }
(void) osmmfree( (char *) data );
}

/*++++++++++++++++++++++++++++++
.IDENTifer   PLCONI
.PURPOSE     produce plot information for a two-dimensional contour plot
     input:   int   plmode    plot mode, see PMODE in PLISTAT
              char  *name     data frame name
              char  *ident    ascii identifier of image
              float *clevl    values of the contours
              int    *ctype   ctype >= 0: array with line-types of contours
                              ctype < 0 : array with colors of contours
              int   nlevl     number of contour levels
--------------------------------*/
void PLCONI( plmode, name, ident, clevl, ctype, nlevl )
int   plmode, nlevl, *ctype; 
char  *name, *ident;
float *clevl;
{
int    actvals, color, ii;
float  one, ssize, tsize, x1, x2, y1, y2, xt, yt, yh, xstr;
float  scale[2], mnmx[2], xl[3], yl[3], image[4], clpl[4], wndl[4];
char   text[81];

/*
 * only plot mode is 2
 */
if ( plmode != 2 ) return;

/*
 * get the symbol and character dimensions, from the MIDAS keywords
 */
PCKRDR( "SSIZE", 1, &actvals, &ssize );
PCKRDR( "TSIZE", 1, &actvals, &tsize );

/*
 * if the symbol size or the text size is not equal to 1, set it to 1
 * and call PCTSET to get the proper sizes for the MIDAS layout
 */
if ( ssize != 1.0 || tsize != 1.0 )
   { one = 1.0;
     PCKWRR( "SSIZE", 1, &one );
     PCKWRR( "TSIZE", 1, &one );
   }
PCTSET();

AG_SSET( "norm");
AG_SSET( "linx");
AG_SSET( "liny");
(void) AG_RGET( "clpl", clpl );
(void) AG_RGET( "wndl", wndl );
x1 = clpl[1] + 0.01;
x2 = 1.0;
y1 = 0.0;
y2 = clpl[3];
AG_CDEF( x1, x2, y1, y2);
AG_WDEF( 0.0, 1.0, 0.0, 1.0 );

/*
 * plot MIDAS logo
 */
PLLOGI( &xt, &yt );

/* 
 * set character height
 */
AG_SSET( "sydi=0.75;chdi=0.75,0.75" );
AG_TGET( "M", xl, yl );
yh = 2.0 * yl[1];

/*
 * plot user name
 */
PLUSER( text );
AG_GTXT( xt, yt, text, 1 );

/*
 * name of frame
 */
yt -= 2*yh;
if ( strlen( name ) > (size_t) 12 )
   { AG_GTXT( xt, yt, "Frame:", 1 );
     yt -= yh;
     AG_GTXT( xt, yt, name, 1 ); 
   }
else
   { (void) sprintf( text, "Frame: %s", name );
     AG_GTXT( xt, yt, text, 1 );
   }
/*
 * identification
 */
if ( strlen( ident ) > (size_t) 0 )
   { yt -= 2*yh;
     AG_GTXT( xt, yt, "Identification:", 1 );
     yt -= yh;
     AG_GTXT( xt, yt, ident, 1 );
   }
/*
 * area
 */
PCKRDR( "PIXEL", 4, &actvals, image );
yt -= 2*yh;
AG_GTXT( xt, yt, "Area:", 1 );
yt -= yh;
(void) sprintf( text, "X: %.0f to %.0f", image[0], image[1] );
AG_GTXT( xt, yt, text, 1 );
yt -= yh;
(void) sprintf( text, "Y: %.0f to %.0f", image[2], image[3] );
AG_GTXT( xt, yt, text, 1 );
/*
 * scales
 */
PCKRDR( "SCALES", 2, &actvals, scale );
yt -= 2*yh;
AG_GTXT( xt, yt, "Scales:", 1 );
yt -= yh;
(void) sprintf( text, "X: %-.6g", scale[0] );
AG_GTXT( xt, yt, text, 1 );
yt -= yh;
(void) sprintf( text, "Y: %-.6g", scale[1] );
AG_GTXT( xt, yt, text, 1 );
/*
 * minimum and maximum
 */
PCKRDR( "ZWNDL", 2, &actvals, mnmx );
yt -= 2*yh;
(void) sprintf( text, "Min: %-.3g ", mnmx[0] );
AG_GTXT( xt, yt, text, 1 );
yt -= yh;
(void) sprintf( text, "Max: %-.3g ", mnmx[1] );
AG_GTXT( xt, yt, text, 1 );

/*
 * contour levels
 */
yt -= 2*yh;
AG_GTXT( xt, yt, "Contour levels:", 1 );
color = 1;
for( ii = 0; ii < nlevl; ii++ )
   { xstr = xt + 11 * xl[1];
     if ( ii % 2 == 0 )
        { xstr = xt;
          yt  -= yh;
        }

     if ( ctype[ii] < 0 && color != ctype[ii] )
        { (void) sprintf( text, "color=%1d", -ctype[ii] );
          color = ctype[ii];
          AG_SSET( text );
        }

     (void) sprintf( text, "%8.5g", clevl[ii] );
     AG_GTXT( xstr, yt, text, 1 );
   }
if ( color != 1 ) AG_SSET( "color=1" );

yt -= 2.*yh;
PLDATI( &xt, &yt );

/*
 * back to standard sizes
 */
AG_CDEF( clpl[0], clpl[1], clpl[2], clpl[3] );
AG_WDEF( wndl[0], wndl[1], wndl[2], wndl[3] );

/*
 * reset if necessary the symbol and text size
 */
if ( ssize != 1.0 || tsize != 1.0 )
   { PCKWRR( "SSIZE", 1, &ssize );
     PCKWRR( "TSIZE", 1, &tsize );
     PCTSET();
   }
}