/* $Header: /afs/athena.mit.edu/astaff/project/atdev/src/plotter/RCS/ContourPlot.c,v 3.6 95/05/18 12:53:10 dot Exp Locker: dot $ */

/*******************************************************************
  Copyright (C) 1990 by the Massachusetts Institute of Technology

   Export of this software from the United States of America is assumed
   to require a specific license from the United States Government.
   It is the responsibility of any person or organization contemplating
   export to obtain such a license before exporting.

WITHIN THAT CONSTRAINT, permission to use, copy, modify, and
distribute this software and its documentation for any purpose and
without fee is hereby granted, provided that the above copyright
notice appear in all copies and that both that copyright notice and
this permission notice appear in supporting documentation, and that
the name of M.I.T. not be used in advertising or publicity pertaining
to distribution of the software without specific, written prior
permission.  M.I.T. makes no representations about the suitability of
this software for any purpose.  It is provided "as is" without express
or implied warranty.

***************************************************************** */

#include <stdio.h>
#include <X11/Intrinsic.h>
#include <X11/StringDefs.h>

#ifdef _AtDevelopment_
#include "Scale.h"    
#include "ContourPlotP.h"
#include "AtConverters.h"
#else
#include <At/Scale.h>
#include <At/ContourPlotP.h>
#include <At/AtConverters.h>
#endif

#ifdef SYSV
#define bcopy(s1, s2, l) memcpy (s2, s1, l)
#endif /* SYSV */

#define offset(field) XtOffset(AtContourPlotWidget, field)
static XtResource resources[] = {
  { XtNcontours, XtCContours, XtRInt, sizeof(int),
      offset(contourplot.contours), XtRImmediate, (caddr_t) 10 },
  { XtNxNumPoints, XtCXNumPoints, XtRInt, sizeof(int),
      offset(contourplot.xNumPoints), XtRImmediate, (caddr_t) 0 },
  { XtNyNumPoints, XtCYNumPoints, XtRInt, sizeof(int),
      offset(contourplot.yNumPoints), XtRImmediate, (caddr_t) 0 },
  { XtNxMin, XtCXMin, XtRDouble, sizeof(double),
      offset(contourplot.xMin), XtRString, "0.0" },
  { XtNxMax, XtCXMax, XtRDouble, sizeof(double),
      offset(contourplot.xMax), XtRString, "1.0" },
  { XtNyMin, XtCYMin, XtRDouble, sizeof(double),
      offset(contourplot.yMin), XtRString, "0.0" },
  { XtNyMax, XtCYMax, XtRDouble, sizeof(double),
      offset(contourplot.yMax), XtRString, "1.0" },
  { XtNxPoints, XtCXPoints, XtRDoubleArray, sizeof(double *),
      offset(contourplot.xPoints), XtRImmediate, NULL },
  { XtNyPoints, XtCYPoints, XtRDoubleArray, sizeof(double *),
      offset(contourplot.yPoints), XtRImmediate, NULL },
  { XtNzPoints, XtCZPoints, XtRDoubleArray, sizeof(double *),
      offset(contourplot.zPoints), XtRImmediate, NULL },
};
#undef offset

static void Initialize(AtContourPlotWidget, AtContourPlotWidget);
static void Destroy(AtContourPlotWidget);
static Boolean SetValues(AtContourPlotWidget, AtContourPlotWidget, 
			 AtContourPlotWidget);
static void Recompute(AtPlotWidget, AtScale *, AtScale *);
static void Draw(Display *, Window, AtPlotWidget, AtScale *, AtScale *,
		 Region);
static void DrawPS(FILE *, AtPlotWidget, AtScale * , AtScale *);
static void DrawIcon(Display *, Window, AtPlotWidget, int, int, int, int);
static void DrawIconPS(FILE *, AtPlotWidget, int, int, int, int);
static void ComputeContourLines(AtContourPlotWidget);
static void Triangle(AtContourPlotWidget, int, int, MidPoint *);
static double Intersect(double, double, double);
static void AddLine(AtContourPlotWidget, RealSegment *);
static void minmax(double *, int, double *, double *);
static double GetValue(AtContourPlotWidget, int, int);

AtContourPlotClassRec atContourPlotClassRec = {
{ /* core part */
    /* superclass         */    (WidgetClass)&atPlotClassRec,
    /* class_name         */    "AtContourPlot",
    /* widget_size        */    sizeof(AtContourPlotRec),
    /* class_initialize   */    NULL,
    /* class_part_initialize*/  NULL,
    /* class_inited       */    FALSE,
    /* initialize         */    (XtInitProc) Initialize,
    /* initialize_hook    */    NULL,
    /* pad                */    NULL,
    /* pad                */    NULL,
    /* pad                */    0,
    /* resources          */    resources,
    /* num_resources      */    XtNumber(resources),
    /* xrm_class          */    NULLQUARK,
    /* pad                */    FALSE,
    /* pad                */    FALSE,
    /* pad                */    FALSE,
    /* pad                */    FALSE,
    /* destroy            */    (XtWidgetProc) Destroy,
    /* pad                */    NULL,
    /* pad                */    NULL,
    /* set_values         */    (XtSetValuesFunc) SetValues,
    /* set_values_hook    */    NULL,
    /* pad                */    NULL,
    /* get_values_hook    */    NULL,
    /* pad                */    NULL,
    /* version            */    XtVersion,
    /* callback_offsets   */    NULL,
    /* pad                */    NULL,
    /* pad                */    NULL,
    /* pad                */    NULL,
    /* extension            */  NULL
},
/* AtPlotClassPart initialization */
{
    /* draw()       */   Draw,
    /* drawIcon()   */   DrawIcon,
    /* resize()     */   Recompute,
    /* rescale()    */   Recompute,
    /* drawPS()     */	 DrawPS,
    /* drawIconPS() */   DrawIconPS,
    /* checkhit()   */   NULL,
}
};

WidgetClass atContourPlotWidgetClass = (WidgetClass) &atContourPlotClassRec;

#define AtWarning(w,msg) XtWarning(msg)
#define MINSEGS 50
#define XVALUE 1
#define YVALUE 2

/**********************************************************************/

static void SetBoundingBox(AtContourPlotWidget w)
{
  Arg arg;
  BoundingBox b;

  if (w->contourplot.xPoints)
    minmax(w->contourplot.xPoints, w->contourplot.xNumPoints,
	   &b.xmin, &b.xmax);
  else  {
    b.xmin = w->contourplot.xMin;
    b.xmax = w->contourplot.xMax;
  }
  if (w->contourplot.yPoints)
    minmax(w->contourplot.yPoints, w->contourplot.yNumPoints,
	   &b.ymin, &b.ymax);
  else  {
    b.ymin = w->contourplot.yMin;
    b.ymax = w->contourplot.yMax;
  }
  
  XtSetArg(arg, XtNboundingBox, &b);
  XtSetValues(w, &arg, 1);
}

/**********************************************************************/

static void Initialize(AtContourPlotWidget request, AtContourPlotWidget new)
{
  /*** XXX Do error checking here ***/
  
  if (new->contourplot.xPoints != NULL)  {
    new->contourplot.xpts =
      (double *)XtMalloc(sizeof(double) * new->contourplot.xNumPoints);
    bcopy((char *) new->contourplot.xPoints, (char *)new->contourplot.xpts,
	  sizeof(double) * new->contourplot.xNumPoints);
  }
  else
    new->contourplot.xpts = NULL;

  if (new->contourplot.yPoints != NULL)  {
    new->contourplot.ypts =
      (double *)XtMalloc(sizeof(double) * new->contourplot.yNumPoints);
    bcopy((char *)new->contourplot.yPoints, (char *)new->contourplot.ypts,
	  sizeof(double) * new->contourplot.yNumPoints);
  }
  else
    new->contourplot.ypts = NULL;
  
  new->contourplot.zpts = 
    (double *) XtMalloc(sizeof(double) * new->contourplot.xNumPoints *
			new->contourplot.yNumPoints);
  bcopy((char *) new->contourplot.zPoints, (char *) new->contourplot.zpts,
	sizeof(double) * new->contourplot.xNumPoints * 
	new->contourplot.yNumPoints);

  new->contourplot.realsegs = NULL;
  new->contourplot.segments = NULL;
  new->contourplot.nsegs = 0;
  new->contourplot.seglistsize = 0;
  new->contourplot.valid = False;
  ComputeContourLines(new);
  new->contourplot.contour_valid = True;

  /* figure out our bounding box */
  SetBoundingBox(new);
}

/**********************************************************************/

static void Destroy(AtContourPlotWidget w)
{
    XtFree(w->contourplot.xpts);
    XtFree(w->contourplot.ypts);
    XtFree(w->contourplot.zpts);
}

/**********************************************************************/

static Boolean SetValues(AtContourPlotWidget current,
			 AtContourPlotWidget request,
			 AtContourPlotWidget new)
{
#define Changed(field) (new->contourplot.field != current->contourplot.field)
  Boolean xchanged = False, ychanged = False;
  int i;
  Arg a;

  if (Changed(contours) || Changed(xMin) || Changed(xMax) || 
      Changed(yMin) || Changed(yMax))  {
    new->contourplot.valid = False;
    new->contourplot.contour_valid = False;
    new->plot.redisplay = True;
  }

  if (Changed(xPoints))
    xchanged = True;
  else
    for (i = 0; i < new->contourplot.xNumPoints; i++) 
      if ((current->contourplot.xpts && new->contourplot.xPoints) &&
          (current->contourplot.xpts[i] != new->contourplot.xPoints[i]))  {
	xchanged = True;
	break;
      }
  if (Changed(yPoints))
    ychanged = True;
  else
    for (i = 0; i < new->contourplot.yNumPoints; i++)
      if ((current->contourplot.ypts && new->contourplot.yPoints) &&
          (current->contourplot.ypts[i] != new->contourplot.yPoints[i]))  {
	ychanged = True;
	break;
      }

  if (Changed(xNumPoints))  {
    new->contourplot.valid = False;
    xchanged = True;
    if (new->contourplot.xPoints != NULL)
      new->contourplot.xpts = (double *) XtRealloc(new->contourplot.xpts, 
		new->contourplot.xNumPoints * sizeof(double));
    else  {
      XtFree(new->contourplot.xpts);
      new->contourplot.xpts = NULL;
    }
  }
  if (Changed(yNumPoints))  {
    new->contourplot.valid = False;
    ychanged = True;
    if (new->contourplot.yPoints != NULL)
      new->contourplot.ypts = (double *) XtRealloc(new->contourplot.ypts, 
		new->contourplot.yNumPoints * sizeof(double));
    else  {
      XtFree(new->contourplot.ypts);
      new->contourplot.ypts = NULL;
    }
  }

  if (xchanged)  {
    new->plot.redisplay = True;
    new->contourplot.contour_valid = False;
    bcopy((char *) new->contourplot.xPoints, (char *) new->contourplot.xpts,
	  sizeof(double) * new->contourplot.xNumPoints);
  }
  if (ychanged)  {
    new->plot.redisplay = True;
    new->contourplot.contour_valid = False;
    bcopy((char *) new->contourplot.yPoints, (char *) new->contourplot.ypts,
	  sizeof(double) * new->contourplot.yNumPoints);
  }

  if (Changed(xNumPoints) || Changed(yNumPoints))  {
    new->contourplot.zpts = (double *) XtRealloc(new->contourplot.zpts, 
              new->contourplot.xNumPoints * new->contourplot.yNumPoints * 
	      sizeof(double));
    bcopy((char *) new->contourplot.zPoints, new->contourplot.zpts,
	  new->contourplot.xNumPoints * sizeof(double) * 
	  new->contourplot.yNumPoints);
  }

  if ((new->plot.redisplay) &&
      (new->object.widget_class == atContourPlotWidgetClass)) {
    SetBoundingBox(new);
    new->plot.redisplay = False;
    XtSetArg(a, XtNneedsRedraw, True);
    XtSetValues(new, &a, 1);        /* set a constraint resource */
  }

    return False;
#undef Changed    
}

/**********************************************************************/
/**********************************************************************/

static void Recompute(AtPlotWidget plot, AtScale *xscale, AtScale *yscale)
{
  AtContourPlotWidget cplot = (AtContourPlotWidget) plot;
  int i;

  if(!cplot->contourplot.contour_valid)  {
    ComputeContourLines(cplot);
    cplot->contourplot.contour_valid = True;
  }

  if (cplot->contourplot.segments == NULL)
    cplot->contourplot.segments = (XSegment *) 
      XtMalloc(cplot->contourplot.nsegs * sizeof(XSegment));

  for(i = 0; i < cplot->contourplot.nsegs; i++)  {
    cplot->contourplot.segments[i].x1 = 
      AtScaleUserToPixel(xscale, cplot->contourplot.realsegs[i].x1);
    cplot->contourplot.segments[i].x2 = 
      AtScaleUserToPixel(xscale, cplot->contourplot.realsegs[i].x2);
    cplot->contourplot.segments[i].y1 = 
      AtScaleUserToPixel(yscale, cplot->contourplot.realsegs[i].y1);
    cplot->contourplot.segments[i].y2 = 
      AtScaleUserToPixel(yscale, cplot->contourplot.realsegs[i].y2);
  }
  cplot->contourplot.valid = True;
}

/**********************************************************************/

static void Draw(Display *dpy, Window win, AtPlotWidget plot,
		 AtScale *xscale, AtScale *yscale, Region region)
{
  AtContourPlotWidget cplot = (AtContourPlotWidget) plot;
  int clipx, clipy;
  static GC tmpGc = NULL;
  Pixmap clip_mask;
  XRectangle clip_rect;
  XSegment *base;
  int size;
  XGCValues values;
  unsigned long valueMask;

  if (cplot->contourplot.valid == False)
    Recompute(plot, xscale, yscale);

  if (tmpGc == (GC)NULL) {
    tmpGc = XCreateGC (dpy, win, 0, NULL);
  }
  valueMask = GCClipXOrigin | GCClipYOrigin | GCClipMask;
  XCopyGC (dpy, cplot->plot.gc, valueMask, tmpGc);

  if (region == NULL)  {
    clip_rect.x = AtScaleGetLowPix(xscale);
    clip_rect.y = AtScaleGetHighPix(yscale);
    clip_rect.width = AtScaleGetHighPix(xscale) - clip_rect.x;
    clip_rect.height = AtScaleGetLowPix(yscale) - clip_rect.y;
    XSetClipRectangles(dpy, cplot->plot.gc, 0,0,&clip_rect, 1, Unsorted);
  }
  else
    XSetRegion(dpy, cplot->plot.gc, region);

  /* do the drawing */
  base = cplot->contourplot.segments;
  size = cplot->contourplot.nsegs;
  while(size > 1000)  {
    XDrawSegments(dpy, win, cplot->plot.gc, base, 1000);
    base += 1000;
    size -= 1000;
  }
  if (size) XDrawSegments(dpy, win, cplot->plot.gc, base, size);

  /* now restore the clip region */
  /* We depend on valueMask being unchanged */
  XCopyGC (dpy, tmpGc, valueMask, cplot->plot.gc);
}

/**********************************************************************/

static void DrawPS(FILE *f, AtPlotWidget plot, AtScale *xs, AtScale *ys)
{
  AtContourPlotWidget cplot = (AtContourPlotWidget) plot;
  int i;

  if(!cplot->contourplot.contour_valid)  {
    ComputeContourLines(cplot);
    cplot->contourplot.contour_valid = True;
  }

  /* adjust scales so I get greater accuracy */
  AtScaleResize(xs, AtScaleGetLowPix(xs)*100, AtScaleGetHighPix(xs)*100);
  AtScaleResize(ys, AtScaleGetLowPix(ys)*100, AtScaleGetHighPix(ys)*100);

  fprintf(f, "%%%% BeginObject: AtContourPlot\ngsave\n");
  fprintf(f, "gsave\n");
  fprintf(f, "%d setlinewidth\n",
	  (cplot->plot.lineWidth == 0) ? 1 : cplot->plot.lineWidth);
  if (cplot->plot.lineStyle != LineSolid)
    fprintf(f, "[%d] 0 setdash\n", cplot->plot.dashLength);
  
  for (i = 0; i < cplot->contourplot.nsegs; i++)  {
    fprintf(f, "%g %g M\n",
	    ((double) AtScaleUserToPixel(xs, cplot->contourplot.realsegs[i].x1) 
	     / 100.0),
	    ((double) AtScaleUserToPixel(ys, cplot->contourplot.realsegs[i].y1)
	     / 100.0));
    fprintf(f, "%g %g L\n", 
	    ((double) AtScaleUserToPixel(xs, cplot->contourplot.realsegs[i].x2)
	     / 100.0),
	    ((double) AtScaleUserToPixel(ys, cplot->contourplot.realsegs[i].y2)
	     / 100.0));
    if ((i % 25) == 0)
      fprintf(f, "stroke\n");
  }
  fprintf(f, "stroke\n");
  fprintf(f, "grestore\n%%%% EndObject\n");

  /* adjust scales back to original value */
  AtScaleResize(xs, AtScaleGetLowPix(xs)/100, AtScaleGetHighPix(xs)/100);
  AtScaleResize(ys, AtScaleGetLowPix(ys)/100, AtScaleGetHighPix(ys)/100);
}

/**********************************************************************/

static void DrawIcon(Display *dpy, Window win, AtPlotWidget pw, int x, 
		     int y, int width, int height)

{
  XSegment segs[2];
  int fourth;

  fourth = width / 4;

  segs[0].y1 = segs[1].y1 = y + height - 1;
  segs[0].y2 = segs[1].y2 = y + 1;

  segs[0].x1 = x + fourth;
  segs[0].x2 = x + (width / 2);
  segs[1].x1 = segs[0].x1 + fourth;
  segs[1].x2 = segs[0].x2 + fourth;

  XDrawSegments(dpy, win, pw->plot.gc, segs, 2);
}

/**********************************************************************/

static void DrawIconPS(FILE *f, AtPlotWidget pw, int x, int y, 
		       int width, int height)

{
  int fourth;

  fourth = width / 4;

  fprintf(f, "%%%% BeginObject: AtContourPlot icon\ngsave\n");
  fprintf(f, "gsave\n");
  if (pw->plot.lineStyle != LineSolid)
    fprintf(f, "[%d] 0 setdash\n", pw->plot.dashLength);

  fprintf(f, "%d %d M\n", x + fourth, y - height + 1);
  fprintf(f, "%d %d L\n", x + (width / 2), y - 1);
  fprintf(f, "%d %d M\n", x + 2*fourth, y - height + 1);
  fprintf(f, "%d %d L\n", x + (width / 2) + fourth, y - 1);
  fprintf(f, "stroke\n");
  fprintf(f, "grestore\n%%%% EndObject\n");
}

/**********************************************************************/

static void ComputeContourLines(AtContourPlotWidget cplot)
{
  MidPoint midpoint;
  int i,j, pt1, pt2, pt3, pt4;

  cplot->contourplot.xinterval = 
    (cplot->contourplot.xMax - cplot->contourplot.xMin) / 
      cplot->contourplot.xNumPoints;
  cplot->contourplot.yinterval = 
    (cplot->contourplot.yMax - cplot->contourplot.yMin) / 
      cplot->contourplot.yNumPoints;

  minmax(cplot->contourplot.zPoints, cplot->contourplot.xNumPoints *
	 cplot->contourplot.yNumPoints, &cplot->contourplot.minz, 
	 &cplot->contourplot.maxz);
  cplot->contourplot.interval = (cplot->contourplot.maxz -
				 cplot->contourplot.minz) / 
				   cplot->contourplot.contours;
  
  for(i = 0; i < cplot->contourplot.yNumPoints - 1; i++)
    for(j = 0; j < cplot->contourplot.xNumPoints - 1; j++)  {
      pt1 = i * cplot->contourplot.xNumPoints + j;
      pt2 = pt1 + 1;
      pt3 = pt1 + cplot->contourplot.xNumPoints;
      pt4 = pt3 + 1;
      midpoint.x = .25 * (GetValue(cplot, XVALUE, pt1) + 
			  GetValue(cplot, XVALUE, pt2) + 
			  GetValue(cplot, XVALUE, pt3) +
			  GetValue(cplot, XVALUE, pt4));
      midpoint.y = .25 * (GetValue(cplot, YVALUE, pt1) + 
			  GetValue(cplot, YVALUE, pt2) +
			  GetValue(cplot, YVALUE, pt3) +
			  GetValue(cplot, YVALUE, pt4));
      midpoint.z = .25 * (cplot->contourplot.zPoints[pt1] + 
			  cplot->contourplot.zPoints[pt2] + 
			  cplot->contourplot.zPoints[pt3] + 
			  cplot->contourplot.zPoints[pt4]);
	
      Triangle(cplot, pt1, pt2, &midpoint);
      Triangle(cplot, pt2, pt4, &midpoint);
      Triangle(cplot, pt1, pt3, &midpoint);
      Triangle(cplot, pt3, pt4, &midpoint);
    }

/*
  target->realsegs =
    (RealSegment*)XtRealloc(target->realsegs,
			    target->nsegs * sizeof(RealSegment));
  fprintf(stderr, "%d segments generated.\n", target->nsegs);
*/
}

/**********************************************************************/

static void Triangle(AtContourPlotWidget cplot, int pt1, int pt2,
		     MidPoint *mp)
{
  RealSegment segment;
  double contourvalue, *smallestz, *biggestz, p1, p2, p3;
  double xp1, yp1, xp2, yp2;

  smallestz = biggestz = &(mp->z);

  if (cplot->contourplot.zPoints[pt1] < *smallestz) 
    smallestz = &(cplot->contourplot.zPoints[pt1]);
  else if (cplot->contourplot.zPoints[pt1] > *biggestz)
    biggestz = &(cplot->contourplot.zPoints[pt1]);

  if (cplot->contourplot.zPoints[pt2] < *smallestz)
    smallestz = &(cplot->contourplot.zPoints[pt2]);
  else if (cplot->contourplot.zPoints[pt2] > *biggestz)
    biggestz = &(cplot->contourplot.zPoints[pt2]);

  if (*smallestz == *biggestz) 
    return; /* if triangle is flat, no
	       contours will be drawn on it */

  xp1 = GetValue(cplot, XVALUE, pt1); xp2 = GetValue(cplot, XVALUE, pt2);
  yp1 = GetValue(cplot, YVALUE, pt1); yp2 = GetValue(cplot, YVALUE, pt2);

  for(contourvalue = cplot->contourplot.minz + (.5*cplot->contourplot.interval); 
      contourvalue < cplot->contourplot.maxz; 
      contourvalue += cplot->contourplot.interval)  {
    if ((*smallestz <= contourvalue) &&
	(*biggestz >= contourvalue))  {
      p1 = Intersect(cplot->contourplot.zPoints[pt1], 
		     cplot->contourplot.zPoints[pt2], contourvalue);
      p2 = Intersect(cplot->contourplot.zPoints[pt2], mp->z, contourvalue);
      p3 = Intersect(mp->z, cplot->contourplot.zPoints[pt1], contourvalue);

      if (p1 == 2)  {
	if (p2 == 2) continue; /* flat triangle */
	if (p2 == -1)  {       /* draw a line from pt1 to pt2 */
	  segment.x1 = xp1; segment.x2 = xp2;
	  segment.y1 = yp1; segment.y2 = yp2;
	  AddLine(cplot, &segment);
	  continue;
	}
      }
      if (p2 == 2)  {          /* we know no other p? can be two */
	segment.x1 = xp2; segment.x2 = mp->x;
	segment.y1 = yp2; segment.y2 = mp->y;
	AddLine(cplot, &segment);
	continue;
      }

      /* at this point, the contour does not fall on an edge */
      if (p1 >= 0)  {
	segment.x1 = ((1-p1)*xp1 + p1 * xp2);
	segment.y1 = ((1-p1)*yp1 + p1 * yp2);
	if (p2 >= 0)  {
	  segment.x2 = ((1-p2)* xp2 + p2 * mp->x);
	  segment.y2 = ((1-p2)*yp2+ p2 * mp->y);
	  AddLine(cplot, &segment);
	  continue;
	}
	else if (p3 >= 0)  {
	  segment.x2 = ((1-p3)* mp->x + p3 * xp1);
	  segment.y2 = ((1-p3)* mp->y + p3 * yp1);
	  AddLine(cplot, &segment);
	  continue;
	}
      }
      if (p2 >= 0)  {
	segment.x1 = ((1-p2)*xp2+ p2 * mp->x);
	segment.y1 = ((1-p2)*yp2+ p2 * mp->y);
	segment.x2 = ((1-p3)* mp->x + p3 * xp1);
	segment.y2 = ((1-p3)* mp->y + p3 * yp1);
	AddLine(cplot, &segment);
	continue;
      }
    }
  }
}

/**********************************************************************/

static double Intersect(double zvalue1, double zvalue2, double cvalue)
{
  double percent;
  if (zvalue1 == zvalue2)  {
    if (zvalue1 == cvalue) return(2.0);
    else return(-1.0);
  }
  percent = (cvalue - zvalue1) / (zvalue2 - zvalue1);
  if ((percent >= 0.0) && (percent <= 1.0)) return(percent);
  else return (-1.0);
}

/**********************************************************************/

static void AddLine(AtContourPlotWidget cplot, RealSegment *segment)
{
  if ((segment->x1 == segment->x2) && (segment->y1 == segment->y2))
    return;

  if (cplot->contourplot.realsegs == NULL)  {
    cplot->contourplot.realsegs =
      (RealSegment *) XtMalloc(MINSEGS * sizeof(RealSegment));
    cplot->contourplot.seglistsize = MINSEGS;
    cplot->contourplot.nsegs = 0;
  }

  if (cplot->contourplot.nsegs == cplot->contourplot.seglistsize)  {
    cplot->contourplot.seglistsize *= 1.5;
    cplot->contourplot.realsegs =
      (RealSegment *) XtRealloc(cplot->contourplot.realsegs,
				cplot->contourplot.seglistsize * 
				sizeof(RealSegment));
  }

  bcopy(segment, &(cplot->contourplot.realsegs[cplot->contourplot.nsegs]),
	sizeof(RealSegment));
  cplot->contourplot.nsegs++;
}

/**********************************************************************/

/***********************************************
 *********           minmax          ***********
 ***********************************************
  Convenience routine for plots to find the
  maximum and minimum values in a double precision
  array.  The algorithm should be 25% faster
  than the obvious algorithm.
 ***********************************************/
void minmax(double *data, int n, double *minret, double *maxret)
{
  register int i;

  *maxret = *minret = data[0];

  for(i = 1; i < n - 1; i += 2)  {
    if (data[i] < data[i+1])  {
      if (data[i] < *minret) *minret = data[i];
      if (data[i+1]> *maxret) *maxret = data[i+1];
    }
    else  {
      if (data[i] > *maxret) *maxret = data[i];
      if (data[i+1] < *minret) *minret = data[i+1];
    }
  }

  if (i == n - 1)  {
    if (data[i] > *maxret) *maxret = data[i];
    if (data[i] < *minret) *minret = data[i];
  }
  return;
}

/**********************************************************************/

static double GetValue(AtContourPlotWidget cplot, int xory, int index)
{
  if (xory == XVALUE)  {
    if (cplot->contourplot.xPoints != NULL) 
      return cplot->contourplot.xPoints[index];
    else return (double) (index % cplot->contourplot.xNumPoints);
  }
  else  {
    if (cplot->contourplot.yPoints != NULL) 
      return cplot->contourplot.yPoints[index];
    else return (double) (index / cplot->contourplot.xNumPoints);
  }
}