/*  pdnmesh - a 2D finite element solver
    Copyright (C) 2001-2005 Sarod Yatawatta <sarod@users.sf.net>  
  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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  $Id: contour.c,v 1.36 2005/02/18 18:27:10 sarod Exp $
*/

/* Werner Hoch provided code to produce color in EPS files
 */


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "types.h"


#define LEGEND_WIDTH 40 
#define LEGEND_HEIGHT DEFAULT_HEIGHT
#define LEGEND_BORDER 5

int contour_levels; /* no of contour levels */
int current_plotting_contour; /* number in z array */
/* current max. value of gradient */
MY_DOUBLE max_gradient=0.01; /* initially a low one */


/* contour plotting */
static int 
plot_contour(triangle *t, double potential, Mesh MM)
{
  /* variables for contour plotting */
  int H,L,E;
  
  MY_INT p1,p2,p3;
  MY_DOUBLE x1,x2,y1,y2;

#ifdef DEBUG
  printf("contour for triangle %d\n",t->n);
#endif
		/* initialize the following to ignore compiler warnings */
 p1=p2=p3=0;
	x1=x2=y1=y2=0;
  
  /* now determine the contours if any */
  /* first determine point potentials */
  E=H=L=0;
  if ( Mzz(t->p0,current_plotting_contour,MM) == potential )
    E++;
  else if ( Mzz(t->p0,current_plotting_contour,MM) < potential )
    L++;
  else
    H++;
  if ( Mzz(t->p1,current_plotting_contour,MM) == potential )
    E++;
  else if ( Mzz(t->p1,current_plotting_contour,MM) < potential )
    L++;
  else
    H++;
  if ( Mzz(t->p2,current_plotting_contour,MM) == potential )
    E++;
  else if ( Mzz(t->p2,current_plotting_contour,MM) < potential )
    L++;
  else
    H++;
  
  /* now according to the values of E,H and L */
  /* fit the general case */
  if ( E == 3 ) { /* 3 contours */
    /* we do not plot this */
#ifdef DEBUG
    glBegin(GL_LINE_LOOP);
    /* get first point */
    glVertex2f( (GLfloat)Mx(t->p0,MM), (GLfloat)My(t->p0,MM));
    /* get second point */
    glVertex2f( (GLfloat)Mx(t->p1,MM), (GLfloat)My(t->p1,MM));  
    /* get third point */
    glVertex2f( (GLfloat)Mx(t->p2,MM), (GLfloat)My(t->p2,MM));
    glEnd();
#endif
#ifdef DEBUG
    printf("3 lines\n");
#endif
    return(3);
  }
  if ( E == 2 ) { /* 1 line */
    /* generalize to p1,p2 */
    if ( Mzz(t->p0,current_plotting_contour,MM) != potential ) {
      p1 = t->p1; p2 = t->p2;
    } 
    else if ( Mzz(t->p1,current_plotting_contour,MM) != potential ) {
      p1 = t->p2; p2 = t->p0;
    } 
    else if ( Mzz(t->p2,current_plotting_contour,MM) != potential ) {
      p1 = t->p0; p2 = t->p1;
    } 
    else
      return(-1); /*error*/
    /* now plot the general line */
     glBegin(GL_LINE_LOOP);
     /* get first point */
     glVertex2f( (GLfloat)Mx(p1,MM), (GLfloat)My(p1,MM));
    /* get second point */
    glVertex2f( (GLfloat)Mx(p2,MM), (GLfloat)My(p2,MM));  
    glEnd();
#ifdef DEBUG
    printf("1 line boundary\n");
#endif
    return(1);
  }
  if ( (E == 1) && ( H == 1 ) && ( L == 1 ) ) {
    /* againe one line generalize */
    if ( Mzz(t->p0,current_plotting_contour,MM) == potential ) {
      p1 = t->p1; p2 = t->p2; p3 = t->p0;
    } 
    else if ( Mzz(t->p1,current_plotting_contour,MM) == potential ) {
      p1 = t->p2; p2 = t->p0; p3 = t->p1;
    } 
    else if ( Mzz(t->p2,current_plotting_contour,MM) == potential ) {
      p1 = t->p0; p2 = t->p1; p3 = t->p2;
    } 
    else 
      return(-1);
    /* now interpolate p1,p2 to find the correct point */
    if ( Mzz(p1,current_plotting_contour,MM) == Mzz(p2,current_plotting_contour,MM) )
      return(-1); /* error */
    x1 = Mx(p2,MM) 
      +( Mx(p1,MM)-Mx(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
    y1 = My(p2,MM)
      +( My(p1,MM)-My(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
    
    /* now draw line between p3,(x1,y1) */
    glBegin(GL_LINE_LOOP);
    /* get first point */
    glVertex2f( (GLfloat)Mx(p3,MM), (GLfloat)My(p3,MM));
    glVertex2f( (GLfloat)x1,(GLfloat) y1 );
    glEnd();
#ifdef DEBUG
    printf("1 line equal\n");
#endif
    return(1);
  }
  
  if (( ( H == 2 ) && ( L ==1 ) ) || ( ( H == 1 ) && ( L == 2 ))) {
    /* two lines, need to generalize */
    
    if ( ( H == 2 ) && ( L ==1 ) ) {
      /* find the lowest point, assign it to p1 */
      if ( Mzz(t->p0,current_plotting_contour,MM) < potential ) {
	p1 = t->p0; p2 = t->p1; p3 = t->p2;
      } else  
      if ( Mzz(t->p1,current_plotting_contour,MM) < potential ) {
	p1 = t->p1; p2 = t->p2; p3 = t->p0;
      } else
      if ( Mzz(t->p2,current_plotting_contour,MM) < potential ) {
	p1 = t->p2; p2 = t->p0; p3 = t->p1;
      } else 
	return(-1); /* error */
#ifdef DEBUG
      printf("1 line case1\n");
#endif
    }
    if ( ( H == 1 ) && ( L == 2 ) ) {
      /* find the highest point, assign it to p1 */
      if ( Mzz(t->p0,current_plotting_contour,MM) > potential ) {
	p1 = t->p0; p2 = t->p1; p3 = t->p2;
      } else
      if ( Mzz(t->p1,current_plotting_contour,MM) > potential ) {
	p1 = t->p1; p2 = t->p2; p3 = t->p0;
      } else
      if ( Mzz(t->p2,current_plotting_contour,MM) > potential ) {
	p1 = t->p2; p2 = t->p0; p3 = t->p1;
      } else
	return(-1);
#ifdef DEBUG
      printf("1 line case2\n");
#endif
    }
    /* now interpolate between (p1,p2) and (p1,p3) */
  
    /* first p1,p2 interpolation */
    if ( Mzz(p1,current_plotting_contour,MM) != Mzz(p2,current_plotting_contour,MM) ) {
      x1 = Mx(p2,MM) 
	+ (Mx(p1,MM)-Mx(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
      y1 = My(p2,MM)
	+ (My(p1,MM)-My(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
    }
    /* now p1,p3 interpolation */
    if ( Mzz(p1,current_plotting_contour,MM) != Mzz(p3,current_plotting_contour,MM) ) {
    x2 = Mx(p3,MM) 
      + (Mx(p1,MM)-Mx(p3,MM))*(potential-Mzz(p3,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p3,current_plotting_contour,MM));
    y2 = My(p3,MM)
      + (My(p1,MM)-My(p3,MM))*(potential-Mzz(p3,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p3,current_plotting_contour,MM));
  
    }
    /* now draw line (x1,y1) (x2,y2) */
    glBegin(GL_LINE_LOOP);
    /* get first point */
    glVertex2f( (GLfloat)x1, (GLfloat) y1 );
    glVertex2f( (GLfloat)x2, (GLfloat) y2 );
    glEnd();
    return(1); 

  }
  
  /* if we reach here, no contour */
  return(0);
}

/* 3d plot with potential */
static int 
plot_contour_in_3d(triangle *t, double potential,Mesh MM)
{
  /* variables for contour plotting */
  int H,L,E;
  
  unsigned long int p1,p2,p3;
  long double x1,x2,y1,y2;
  double elevation;/*=potential/g_maxpot;*/
	if ( ABS(g_maxpot)<ABS(g_minpot) ) {
     elevation=potential/ABS(1.1*g_minpot);
	} else {
     elevation=potential/ABS(1.1*g_maxpot);
	}
#ifdef DEBUG
  printf("contour for triangle %d\n",t->n);
#endif
  p1=p2=p3=0;
	 x1=x2=y1=y2=0;

  
  /* now determine the contours if any */
  /* first determine point potentials */
  E=H=L=0;
  if ( Mzz(t->p0,current_plotting_contour,MM) == potential )
    E++;
  else if ( Mzz(t->p0,current_plotting_contour,MM) < potential )
    L++;
  else
    H++;
  if ( Mzz(t->p1,current_plotting_contour,MM) == potential )
    E++;
  else if ( Mzz(t->p1,current_plotting_contour,MM) < potential )
    L++;
  else
    H++;
  if ( Mzz(t->p2,current_plotting_contour,MM) == potential )
    E++;
  else if ( Mzz(t->p2,current_plotting_contour,MM) < potential )
    L++;
  else
    H++;
  
  /* now according to the values of E,H and L */
  /* fit the general case */
  if ( E == 3 ) { /* 3 contours */
    /* we do not plot this */
#ifdef DEBUG
    glBegin(GL_LINE_LOOP);
    /* get first point */
    glVertex3f( (GLfloat)Mx(t->p0,MM), (GLfloat)My(t->p0,MM),(GLfloat)elevation);
    /* get second point */
    glVertex3f( (GLfloat)Mx(t->p1,MM), (GLfloat)My(t->p1,MM),(GLfloat)elevation);  
    /* get third point */
    glVertex3f( (GLfloat)Mx(t->p2,MM), (GLfloat)My(t->p2,MM),(GLfloat)elevation);
    glEnd();
#endif
#ifdef DEBUG
    printf("3 lines\n");
#endif
    return(3);
  }
  if ( E == 2 ) { /* 1 line */
    /* generalize to p1,p2 */
    if ( Mzz(t->p0,current_plotting_contour,MM) != potential ) {
      p1 = t->p1; p2 = t->p2;
    } 
    else if ( Mzz(t->p1,current_plotting_contour,MM) != potential ) {
      p1 = t->p2; p2 = t->p0;
    } 
    else if ( Mzz(t->p2,current_plotting_contour,MM) != potential ) {
      p1 = t->p0; p2 = t->p1;
    } 
    else
      return(-1); /*error*/
    /* now plot the general line */
     glBegin(GL_LINE_LOOP);
     /* get first point */
     glVertex3f( (GLfloat)Mx(p1,MM), (GLfloat)My(p1,MM),(GLfloat)elevation);
    /* get second point */
    glVertex3f( (GLfloat)Mx(p2,MM), (GLfloat)My(p2,MM),(GLfloat)elevation);  
    glEnd();
#ifdef DEBUG
    printf("1 line boundary\n");
#endif
    return(1);
  }
  if ( (E == 1) && ( H == 1 ) && ( L == 1 ) ) {
    /* againe one line generalize */
    if ( Mzz(t->p0,current_plotting_contour,MM) == potential ) {
      p1 = t->p1; p2 = t->p2; p3 = t->p0;
    } 
    else if ( Mzz(t->p1,current_plotting_contour,MM) == potential ) {
      p1 = t->p2; p2 = t->p0; p3 = t->p1;
    } 
    else if ( Mzz(t->p2,current_plotting_contour,MM) == potential ) {
      p1 = t->p0; p2 = t->p1; p3 = t->p2;
    } 
    else 
      return(-1);
    /* now interpolate p1,p2 to find the correct point */
    if ( Mzz(p1,current_plotting_contour,MM) == Mzz(p2,current_plotting_contour,MM) )
      return(-1); /* error */
    x1 = Mx(p2,MM) 
      +( Mx(p1,MM)-Mx(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
    y1 = My(p2,MM)
      +( My(p1,MM)-My(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
    
    /* now draw line between p3,(x1,y1) */
    glBegin(GL_LINE_LOOP);
    /* get first point */
    glVertex3f( (GLfloat)Mx(p3,MM), (GLfloat)My(p3,MM),(GLfloat)elevation);
    glVertex3f( (GLfloat)x1,(GLfloat) y1,(GLfloat)elevation );
    glEnd();
#ifdef DEBUG
    printf("1 line equal\n");
#endif
    return(1);
  }
  
  if (( ( H == 2 ) && ( L ==1 ) ) || ( ( H == 1 ) && ( L == 2 ))) {
    /* two lines, need to generalize */
    
    if ( ( H == 2 ) && ( L ==1 ) ) {
      /* find the lowest point, assign it to p1 */
      if ( Mzz(t->p0,current_plotting_contour,MM) < potential ) {
	p1 = t->p0; p2 = t->p1; p3 = t->p2;
      } else  
      if ( Mzz(t->p1,current_plotting_contour,MM) < potential ) {
	p1 = t->p1; p2 = t->p2; p3 = t->p0;
      } else
      if ( Mzz(t->p2,current_plotting_contour,MM) < potential ) {
	p1 = t->p2; p2 = t->p0; p3 = t->p1;
      } else 
	return(-1); /* error */
#ifdef DEBUG
      printf("1 line case1\n");
#endif
    }
    if ( ( H == 1 ) && ( L == 2 ) ) {
      /* find the highest point, assign it to p1 */
      if ( Mzz(t->p0,current_plotting_contour,MM) > potential ) {
	p1 = t->p0; p2 = t->p1; p3 = t->p2;
      } else
      if ( Mzz(t->p1,current_plotting_contour,MM) > potential ) {
	p1 = t->p1; p2 = t->p2; p3 = t->p0;
      } else
      if ( Mzz(t->p2,current_plotting_contour,MM) > potential ) {
	p1 = t->p2; p2 = t->p0; p3 = t->p1;
      } else
	return(-1);
#ifdef DEBUG
      printf("1 line case2\n");
#endif
    }
    /* now interpolate between (p1,p2) and (p1,p3) */
  
    /* first p1,p2 interpolation */
    if ( Mzz(p1,current_plotting_contour,MM) != Mzz(p2,current_plotting_contour,MM) ) {
      x1 = Mx(p2,MM) 
	+ (Mx(p1,MM)-Mx(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
      y1 = My(p2,MM)
	+ (My(p1,MM)-My(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
    }
    /* now p1,p3 interpolation */
    if ( Mzz(p1,current_plotting_contour,MM) != Mzz(p3,current_plotting_contour,MM) ) {
    x2 = Mx(p3,MM) 
      + (Mx(p1,MM)-Mx(p3,MM))*(potential-Mzz(p3,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p3,current_plotting_contour,MM));
    y2 = My(p3,MM)
      + (My(p1,MM)-My(p3,MM))*(potential-Mzz(p3,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p3,current_plotting_contour,MM));
  
    }
    /* now draw line (x1,y1) (x2,y2) */
    glBegin(GL_LINE_LOOP);
    /* get first point */
    glVertex3f( (GLfloat)x1, (GLfloat) y1,(GLfloat)elevation );
    glVertex3f( (GLfloat)x2, (GLfloat) y2,(GLfloat)elevation );
    glEnd();
    return(1); 

  }
  
  /* if we reach here, no contour */
  return(0);
}

/* to get a color to plot */
int 
get_colour(float* redp, float* greenp, float *bluep, int level, int max_levels)
{
       float cl  = (float)level/(float)max_levels;
/*
		 if ( cl < 0.166666667 ) {
            *redp = 0.0f;
            *greenp = 1.0f;
            *bluep = cl*6.0;
		 } else if ( cl < 0.333333333 ) {
             *redp = 0.0f;
            *greenp = 2.0f - 6.0*cl;
            *bluep = 1.0f;
		 } else if ( cl < 0.5 ) {
              *redp = cl*6.0-3.0;
            *greenp = 0.0f;
            *bluep = 1.0f;
		 } else if ( cl < 0.6666666667 ) {
               *redp = 1.0f;
            *greenp = 0.0f;
            *bluep = 4.0f-cl*6.0;
		 } else if ( cl < 0.833333333 ) {
            *redp = 1.0f;
            *greenp = cl*6.0 -4.0f;
            *bluep = 0.0f;
		 } else {
            *redp = 1.0f;
            *greenp = 1.0f;
            *bluep = cl*6.0-5.0f;
		 }
*/
		 if ( cl < 0.25 ) {
            *redp = 0.0f;
            *greenp = cl*4.0f;
            *bluep = 1.0;
		 } else if ( cl < 0.5 ) {
             *redp = 0.0f;
            *greenp = 1.0f;
            *bluep = 2.0f -cl*4.0f;
		 } else if ( cl < 0.75 ) {
              *redp = cl*4.0f-2.0f;
            *greenp = 1.0f;
            *bluep = 0.0f;
		 } else  {
               *redp = 1.0f;
            *greenp = 4.0f - cl*4.0f;
            *bluep = 0.0f;
		 }	

     /* blue - green */
		  /*
   if ( cl < 1.5 ) { 
   *redp = 0.5;
   *greenp = cl*0.5;
   *bluep = 1.0 - *greenp;
	*/
/* green - red */
		  /*
  } else  { 
   *redp = 0.25 + (cl-1.5)*0.5;
   *bluep = 0.0;
   *greenp = 1.0 - *redp;
  }
  */
  
		 /*
  *redp = (float)level/(float)max_levels;
  *greenp = 0.7f;
  *bluep = 1.0f - *redp;
  */
  return(0);

}

/* function to plot the contour */
void 
plot_contour_all(Mesh MM)
{
  triangle *t;
  int j; /* for contour calculation */
  float r,g,b; /* colours */
  double step =  (g_maxpot-g_minpot)/(float)contour_levels; /* contour step */

  
  
#ifdef DEBUG
  printf("plotting from %lf to %lf at %lf intervals\n",g_maxpot,g_minpot,step);
#endif

		DAG_traverse_list_reset(&dt);
		t=DAG_traverse_prune_list(&dt);
		while(t) {
    if ( t && t->status==LIVE ) {
      for ( j = 0; j < contour_levels ; j++ ) {
	/* get the colour */
	get_colour(&r,&g,&b,j, contour_levels);
	/* set the colour */
	glColor3f((GLfloat)r, (GLfloat)g, (GLfloat)b);
	plot_contour(t,g_minpot+j*step,MM);
      }
         
    }
				
		t=DAG_traverse_prune_list(&dt);
  }

}

void 
plot_contour_all_in_3d(Mesh MM)
{
  triangle *t;
  int j; /* for contour calculation */
  float r,g,b; /* colours */
  double step =  (g_maxpot-g_minpot)/(float)contour_levels; /* contour step */

  
  
#ifdef DEBUG
  printf("plotting from %lf to %lf at %lf intervals\n",g_maxpot,g_minpot,step);
#endif

		DAG_traverse_list_reset(&dt);
		t=DAG_traverse_prune_list(&dt);
		while(t) {
    if ( t && t->status==LIVE ) {
      for ( j = 0; j < contour_levels ; j++ ) {
	/* get the colour */
/*	get_colour(&r,&g,&b,j, contour_levels); */
  r=g=b=0.0; 
	glColor3f((GLfloat)r, (GLfloat)g, (GLfloat)b); 
	plot_contour_in_3d(t,g_minpot+j*step,MM);
      }
         
    }
				
		t=DAG_traverse_prune_list(&dt);
  }

}


/* a function to print the contour of a triangle to a postscript file */
static int 
print_contour(FILE *plotf, triangle *t, MY_DOUBLE potential, MY_DOUBLE myscale,
								float r, float g, float b, int color_flag, Mesh MM)
{
/* variables for contour plotting */
  int H,L,E;
  
  MY_INT p1,p2,p3;
  MY_DOUBLE x1,x2,y1,y2;
  
  /* note: we use shorthand here as follows:
     /l {lineto} bind def 
     /m {moveto} bind def 
     /n {newpath} bind def 
     /s {stroke} bind def 
     /slw {setlinewidth} bind def 
  */


#ifdef DEBUG
  printf("contour for triangle %d\n",t->n);
#endif
		/* initialize the following to ignore compiler warnings */
 p1=p2=p3=0;
x1=x2=y1=y2=0;
  
  /* now determine the contours if any */
  /* first determine point potentials */
  E=H=L=0;
  if ( Mzz(t->p0,current_plotting_contour,MM) == potential )
    E++;
  else if ( Mzz(t->p0,current_plotting_contour,MM) < potential )
    L++;
  else
    H++;
  if ( Mzz(t->p1,current_plotting_contour,MM) == potential )
    E++;
  else if ( Mzz(t->p1,current_plotting_contour,MM) < potential )
    L++;
  else
    H++;
  if ( Mzz(t->p2,current_plotting_contour,MM) == potential )
    E++;
  else if ( Mzz(t->p2,current_plotting_contour,MM) < potential )
    L++;
  else
    H++;
  
  /* now according to the values of E,H and L */
  /* fit the general case */
  if ( E == 3 ) { /* 3 contours */
    /* we do not plot this */
#ifdef DEBUG
    printf("3 lines\n");
#endif
    return(3);
  }
  if ( E == 2 ) { /* 1 line */
    /* generalize to p1,p2 */
    if ( Mzz(t->p0,current_plotting_contour,MM) != potential ) {
      p1 = t->p1; p2 = t->p2;
    } 
    else if ( Mzz(t->p1,current_plotting_contour,MM) != potential ) {
      p1 = t->p2; p2 = t->p0;
    } 
    else if ( Mzz(t->p2,current_plotting_contour,MM) != potential ) {
      p1 = t->p0; p2 = t->p1;
    } 
    else
      return(-1); /*error*/
		/* set color */
		if ( color_flag )
			fprintf(plotf,"%f %f %f srgb\n",r,g,b);
   /* now plot the general line */
    fprintf(plotf,"n ");
    /* get first point */
    fprintf(plotf,"%lf %lf m ",(double) (g_xrange+Mx(p1,MM)*g_xscale)*myscale, (double)(g_yrange + My(p1,MM)*g_yscale)*myscale);
    /* get second point */
    fprintf(plotf,"%lf %lf l ",(double) (g_xrange+Mx(p2,MM)*g_xscale)*myscale, (double)(g_yrange + My(p2,MM)*g_yscale)*myscale);
    fprintf(plotf,"cp ");
    fprintf(plotf,"s\n");
#ifdef DEBUG
    printf("1 line boundary\n");
#endif
    return(1);
  }
  if ( (E == 1) && ( H == 1 ) && ( L == 1 ) ) {
    /* againe one line generalize */
    if ( Mzz(t->p0,current_plotting_contour,MM) == potential ) {
      p1 = t->p1; p2 = t->p2; p3 = t->p0;
    } 
    else if ( Mzz(t->p1,current_plotting_contour,MM) == potential ) {
      p1 = t->p2; p2 = t->p0; p3 = t->p1;
    } 
    else if ( Mzz(t->p2,current_plotting_contour,MM) == potential ) {
      p1 = t->p0; p2 = t->p1; p3 = t->p2;
    } 
    else 
      return(-1);
    /* now interpolate p1,p2 to find the correct point */
    if ( Mzz(p1,current_plotting_contour,MM) == Mzz(p2,current_plotting_contour,MM) )
      return(-1); /* error */
    x1 = Mx(p2,MM) 
      +( Mx(p1,MM)-Mx(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
    y1 = My(p2,MM)
      +( My(p1,MM)-My(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
 
			/* set color */
		if ( color_flag )
			fprintf(plotf,"%f %f %f srgb\n",r,g,b);
	
    /* now draw line between p3,(x1,y1) */
    fprintf(plotf,"n ");
    /* get first point */
    fprintf(plotf,"%lf %lf m ",(double) (g_xrange+Mx(p3,MM)*g_xscale)*myscale, (double)(g_yrange + My(p3,MM)*g_yscale)*myscale);
    /* get second point */
    fprintf(plotf,"%lf %lf l ",(double) (g_xrange+x1*g_xscale)*myscale, (double)(g_yrange + y1*g_yscale)*myscale);
    fprintf(plotf,"cp ");
    fprintf(plotf,"s\n");
   
#ifdef DEBUG
    printf("1 line equal\n");
#endif
    return(1);
  }
  
  if (( ( H == 2 ) && ( L ==1 ) ) || ( ( H == 1 ) && ( L == 2 ))) {
    /* two lines, need to generalize */
    
    if ( ( H == 2 ) && ( L ==1 ) ) {
      /* find the lowest point, assign it to p1 */
      if ( Mzz(t->p0,current_plotting_contour,MM) < potential ) {
	p1 = t->p0; p2 = t->p1; p3 = t->p2;
      } else  
      if ( Mzz(t->p1,current_plotting_contour,MM) < potential ) {
	p1 = t->p1; p2 = t->p2; p3 = t->p0;
      } else
      if ( Mzz(t->p2,current_plotting_contour,MM) < potential ) {
	p1 = t->p2; p2 = t->p0; p3 = t->p1;
      } else 
	return(-1); /* error */
#ifdef DEBUG
      printf("1 line case1\n");
#endif
    }
    if ( ( H == 1 ) && ( L == 2 ) ) {
      /* find the highest point, assign it to p1 */
      if ( Mzz(t->p0,current_plotting_contour,MM) > potential ) {
	p1 = t->p0; p2 = t->p1; p3 = t->p2;
      } else
      if ( Mzz(t->p1,current_plotting_contour,MM) > potential ) {
	p1 = t->p1; p2 = t->p2; p3 = t->p0;
      } else
      if ( Mzz(t->p2,current_plotting_contour,MM) > potential ) {
	p1 = t->p2; p2 = t->p0; p3 = t->p1;
      } else
	return(-1);
#ifdef DEBUG
      printf("1 line case2\n");
#endif
    }
    /* now interpolate between (p1,p2) and (p1,p3) */
  
    /* first p1,p2 interpolation */
    if ( Mzz(p1,current_plotting_contour,MM) != Mzz(p2,current_plotting_contour,MM) ) {
      x1 = Mx(p2,MM) 
	+ (Mx(p1,MM)-Mx(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
      y1 = My(p2,MM)
	+ (My(p1,MM)-My(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
    }
    /* now p1,p3 interpolation */
    if ( Mzz(p1,current_plotting_contour,MM) != Mzz(p3,current_plotting_contour,MM) ) {
    x2 = Mx(p3,MM) 
      + (Mx(p1,MM)-Mx(p3,MM))*(potential-Mzz(p3,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p3,current_plotting_contour,MM));
    y2 = My(p3,MM)
      + (My(p1,MM)-My(p3,MM))*(potential-Mzz(p3,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p3,current_plotting_contour,MM));
  
    }
		/* set color */
		if ( color_flag )
			fprintf(plotf,"%f %f %f srgb\n",r,g,b);

    /* now draw line (x1,y1) (x2,y2) */
    fprintf(plotf,"n ");
    /* get first point */
    fprintf(plotf,"%lf %lf m ",(double) (g_xrange+x1*g_xscale)*myscale, (double)(g_yrange + y1*g_yscale)*myscale);
    /* get second point */
    fprintf(plotf,"%lf %lf l ",(double) (g_xrange+x2*g_xscale)*myscale, (double)(g_yrange + y2*g_yscale)*myscale);
    fprintf(plotf,"cp ");
    fprintf(plotf,"s\n");

  }
  
  /* if we reach here, no contour */
  return(0);

}

static void
print_polygon(polygon *poly, double myscale, FILE *plotf, Mesh MM)
{
 
 if (poly&&poly->head) {
	int i=poly->nedges;
	poly_elist *tempe=poly->head;
  for (;i;i--) {
  fprintf(plotf,"n ");
  fprintf(plotf,"%lf %lf  m ",(double)(Mx(tempe->data.p1,MM)*g_xscale+g_xrange)*myscale, (double)(My(tempe->data.p1,MM)*g_yscale+g_yrange)*myscale);
  fprintf(plotf,"%lf %lf  l ",(double)(Mx(tempe->data.p2,MM)*g_xscale+g_xrange)*myscale, (double)(My(tempe->data.p2,MM)*g_yscale+g_yrange)*myscale);
   fprintf(plotf,"s\n");  
	 tempe=tempe->next;
  }
 }

}


/* printing contour plot */
void
print_contour_all(const char* filename,int color_flag, Mesh MM)
{
 /* color_flag: 0, BW contours, 1: color contours, 2: color contours+legend */
	FILE *plotf;
  MY_DOUBLE myscale,delta;
  int i;
	float r,g,b;
  triangle *t;
  double step =  (g_maxpot-g_minpot)/(float)contour_levels; /* contour step */

	if ((plotf=fopen(filename,"w"))==NULL){
			fprintf(stderr,"%s: %d: could not open file %s\n",__FILE__,__LINE__,filename);
	   return;	
	}


  /* determine the right bound of bounding box */
  if ( g_xrange > g_yrange )
    myscale = 400.0/g_xrange;
  else
    myscale = 400.0/g_yrange;

  /* print the header */
  fprintf(plotf,"%%!PS-Adobe-2.0 EPSF-3.0\n");
  fprintf(plotf,"%%%%Title: %s\n",filename);
  fprintf(plotf,"%%%%Generated By: "PACKAGE" "VERSION"  contour levels: %d, min= "MDF", max= "MDF"\n", contour_levels, g_minpot, g_maxpot);
  /* now bounding box -add border of 20 */
 /* add x margin of 60 for legend 
   the myscale is selected such that max dimension is 820+80(x) of 820 (y)
 */
  if ( color_flag==2 ) {
  fprintf(plotf,"%%%%BoundingBox: 0 0 %4.2lf %4.2lf\n",(double)(2*g_xrange*myscale+20+100), (double)(2*g_yrange*myscale+20) );
  } else {
  fprintf(plotf,"%%%%BoundingBox: 0 0 %4.2lf %4.2lf\n",(double)(2*g_xrange*myscale+20), (double)(2*g_yrange*myscale+20) );
 }
  fprintf(plotf,"%%%%Magnification: 1.0000\n");
  
  

  /* start plotting from 10 10 */
  fprintf(plotf,"10 10 translate\n");


  /* now print some shorter definitions */
  fprintf(plotf,"/l {lineto} bind def\n"); /* lineto */
  fprintf(plotf,"/f {fill} bind def\n"); /* lineto */
  fprintf(plotf,"/m {moveto} bind def\n"); /* moveto */
  fprintf(plotf,"/n {newpath} bind def\n"); /* newpath */
  fprintf(plotf,"/s {stroke} bind def\n"); /* stroke */
  fprintf(plotf,"/slw {setlinewidth} bind def\n"); /* setlinewidth */
  fprintf(plotf,"/cp {closepath} bind def\n"); /* closepath */
	if ( color_flag )
	 fprintf(plotf,"/srgb {setrgbcolor} bind def\n"); /* rgb color */
  
  /* now print contours */
  fprintf(plotf,"1 slw\n");
 DAG_traverse_list_reset(&dt);
		t=DAG_traverse_prune_list(&dt);
		while(t) {
    if ( t && t->status==LIVE ) {
      for ( i = 0; i < contour_levels ; i++ ) {
	/* get the colour */
	get_colour(&r,&g,&b,i, contour_levels);
	/* set the colour */
	     print_contour(plotf,t,g_minpot+i*step,myscale,r,g,b,color_flag,MM);
      }
         
    }
				
		t=DAG_traverse_prune_list(&dt);
  }
 
  /* then plot the boundaries */
  /* set the line thikness */
  fprintf(plotf,"2 slw\n");

	if ( color_flag )
   fprintf(plotf,"0 0 0 srgb\n");
  for ( i = 0; i < bound.npoly; i++ ) {
		print_polygon(&bound.poly_array[i],myscale,plotf,MM);
	}

if ( color_flag==2 ) {
/* draw legend */
  fprintf(plotf,"1 slw\n");
/* box */
  fprintf(plotf,"n %lf %lf m %lf %lf l %lf %lf l %lf %lf l %lf %lf l s\n",(double)(2*g_xrange*myscale+20),0.0, (double)(2*g_xrange*myscale+60),0.0, (double)(2*g_xrange*myscale+60),(double)(2*g_yrange*myscale), (double)(2*g_xrange*myscale+20), (double)(2*g_yrange*myscale), (double)(2*g_xrange*myscale+20),0.0);
/* tick marks */
  if ( contour_levels > 1 ) {
   delta=(double)(2*g_yrange*myscale)/(contour_levels-1);
   for (i=1; i< contour_levels-1; i++) {
    fprintf(plotf,"n %lf %lf m %lf %lf l %lf %lf l %lf %lf l cp\n",(double)(2*g_xrange*myscale+20),(i-1)*delta, (double)(2*g_xrange*myscale+60),(i-1)*delta,(double)(2*g_xrange*myscale+60),(i)*delta,(double)(2*g_xrange*myscale+20),(i)*delta);
	 get_colour(&r,&g,&b,i, contour_levels);
   fprintf(plotf,"%f %f %f srgb f\n",r,g,b);
   fprintf(plotf,"n %lf %lf m %lf %lf l s\n",(double)(2*g_xrange*myscale+20),i*delta,(double)(2*g_xrange*myscale+20+5),i*delta);
   }
    fprintf(plotf,"n %lf %lf m %lf %lf l %lf %lf l %lf %lf l cp\n",(double)(2*g_xrange*myscale+20),(i-1)*delta, (double)(2*g_xrange*myscale+60),(i-1)*delta,(double)(2*g_xrange*myscale+60),(i)*delta,(double)(2*g_xrange*myscale+20),(i)*delta);
	 get_colour(&r,&g,&b,i, contour_levels);
   fprintf(plotf,"%f %f %f srgb f\n",r,g,b);

   /* now text for potentials */
   fprintf(plotf,"/Times-Roman findfont\n");
   fprintf(plotf,"0.5 slw\n");
   fprintf(plotf,"9 scalefont setfont\n");
   fprintf(plotf,"0 0 0 srgb\n");
   for (i=0; i< contour_levels; i++) {
    fprintf(plotf,"n %lf %lf m (%5.2lf) show\n",(double)(2*g_xrange*myscale+62),i*delta,g_minpot+i*step);
   }
  }
 }
 
  fprintf(plotf,"%%%%EOF\n");
  fclose(plotf);
}

/* printing mesh as eps*/
void
print_mesh_eps(const char* filename,Mesh MM)
{

  FILE *plotf;
  MY_DOUBLE myscale;
  int i;
  triangle *t;

	if ((plotf=fopen(filename,"w"))==NULL){
			fprintf(stderr,"%s: %d: could not open file %s\n",__FILE__,__LINE__,filename);
	   return;	
	}


  /* determine the right bound of bounding box */
  if ( g_xrange > g_yrange )
    myscale = 400.0/g_xrange;
  else
    myscale = 400.0/g_yrange;

  /* print the header */
  fprintf(plotf,"%%!PS-Adobe-2.0 EPSF-3.0\n");
  fprintf(plotf,"%%%%Title: %s\n",filename);
  fprintf(plotf,"%%%%Generated By: "PACKAGE" "VERSION"  contour levels: %d, min= "MDF", max= "MDF"\n", contour_levels, g_minpot, g_maxpot);
  /* now bounding box -add border of 20 */
  fprintf(plotf,"%%%%BoundingBox: 0 0 %4.2lf %4.2lf\n",(double)(2*g_xrange*myscale+20), (double)(2*g_yrange*myscale+20) );
  fprintf(plotf,"%%%%Magnification: 1.0000\n");

  

  /* start plotting from 10 10 */
  fprintf(plotf,"10 10 translate\n");


  /* now print some shorter definitions */
  fprintf(plotf,"/l {lineto} bind def\n"); /* lineto */
  fprintf(plotf,"/m {moveto} bind def\n"); /* moveto */
  fprintf(plotf,"/n {newpath} bind def\n"); /* newpath */
  fprintf(plotf,"/s {stroke} bind def\n"); /* stroke */
  fprintf(plotf,"/slw {setlinewidth} bind def\n"); /* setlinewidth */
  fprintf(plotf,"/cp {closepath} bind def\n"); /* closepath */
  
  /* now print contours */
  fprintf(plotf,"1 slw\n");
 DAG_traverse_list_reset(&dt);
		t=DAG_traverse_prune_list(&dt);
		while(t) {
    if ( t && t->status==LIVE ) {
	/* get the colour */
	/*get_colour(&r,&g,&b,j, contour_levels); */
	/* set the colour */
  fprintf(plotf,"n ");
  fprintf(plotf,"%lf %lf  m ",(double)(Mx(t->p0,MM)*g_xscale+g_xrange)*myscale, (double)(My(t->p0,MM)*g_yscale+g_yrange)*myscale);
  fprintf(plotf,"%lf %lf  l ",(double)(Mx(t->p1,MM)*g_xscale+g_xrange)*myscale, (double)(My(t->p1,MM)*g_yscale+g_yrange)*myscale);
   fprintf(plotf,"s\n");  
  fprintf(plotf,"n ");
  fprintf(plotf,"%lf %lf  m ",(double)(Mx(t->p1,MM)*g_xscale+g_xrange)*myscale, (double)(My(t->p1,MM)*g_yscale+g_yrange)*myscale);
  fprintf(plotf,"%lf %lf  l ",(double)(Mx(t->p2,MM)*g_xscale+g_xrange)*myscale, (double)(My(t->p2,MM)*g_yscale+g_yrange)*myscale);
   fprintf(plotf,"s\n");  
  fprintf(plotf,"n ");
  fprintf(plotf,"%lf %lf  m ",(double)(Mx(t->p2,MM)*g_xscale+g_xrange)*myscale, (double)(My(t->p2,MM)*g_yscale+g_yrange)*myscale);
  fprintf(plotf,"%lf %lf  l ",(double)(Mx(t->p0,MM)*g_xscale+g_xrange)*myscale, (double)(My(t->p0,MM)*g_yscale+g_yrange)*myscale);
   fprintf(plotf,"s\n");  
         
    }
				
		t=DAG_traverse_prune_list(&dt);
  }
 
  /* then plot the boundaries */
  /* set the line thikness */
  fprintf(plotf,"2 slw\n");
  for ( i = 0; i < bound.npoly; i++ ) {
		print_polygon(&bound.poly_array[i],myscale,plotf,MM);
	}
    fprintf(plotf,"%%%%EOF\n");
  
  fclose(plotf);
}


/* printing mesh as ACII file*/
void
print_mesh_ascii(const char* filename,Mesh MM)
{
  FILE *outf;
  unsigned int i,j;
  triangle *t;
  unsigned int *renumb,*re_renumb;
  unsigned MY_INT total,unknowns,triangles;

   if ((outf=fopen(filename,"w"))==NULL){
    fprintf(stderr,"%s: %d: could not open file %s\n",__FILE__,__LINE__,filename);
    return;	
   }

  /* allocate memory for arrays */
  if ((renumb = (MY_INT *)calloc(M.count,sizeof(MY_INT)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

  if((re_renumb = (MY_INT *)calloc(M.count,sizeof(MY_INT)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

  total=unknowns=triangles=0;
  /* need to solve for only points referred by triangles */
  DAG_traverse_list_reset(&dt);
  t=DAG_traverse_prune_list(&dt);
  while ( t ) {
    if ( t && t->status==LIVE ) {
      /* mark the points in this triangle */
      /* use re_renumber[] array for this */
      re_renumb[t->p0] = 1;
      re_renumb[t->p1] = 1;
      re_renumb[t->p2] = 1;

      triangles++;
    }
    t=DAG_traverse_prune_list(&dt);
  }

  /* now find unknowns, or INSIDE points */
  j = 0;
  for ( i = 0; i < M.count; i++ )
    if (  Mval(i,MM) == VAR  && re_renumb[i] == 1) {
      renumb[i] = j++; 
      unknowns++;
    }
  /* now renumber the known points */
  total = unknowns;
  for ( i = 0; i < M.count; i++ )
    if ( Mval(i,MM) == FX  && re_renumb[i]==1 ) {
      renumb[i] = j++; 
      total++;
    }

   /* finally give numbers to points 0,1,2 */
   /* we will not use these points */
   renumb[0]=j++;
   renumb[1]=j++;
   renumb[2]=j++;

   /* now construct re-renmbering array */
    i = 0;
    /* do an insertion sort */
    while ( i < M.count ) {
    for ( j = 0; j < M.count; j++ ) {
       if ( renumb[j] == i )
       break;
    }
    re_renumb[i] = j;
    i++;
    }
 
  fprintf(outf,"## Generated by "PACKAGE" "VERSION"\n");
  fprintf(outf,"## Total points %d, (known points %d, unknown points %d) Total triangles %d of %d\n",total,unknowns,total-unknowns,triangles,ntriangles);

   /* now print the unknown point coords, with the number */
   fprintf(outf,"## unknown points: point number | x | y\n");

    for ( i = 0; i < unknowns; i++ )
       fprintf(outf,"%d "MDF" "MDF"\n",i, (Mx(re_renumb[i],MM))*g_xscale-g_xoff, (My(re_renumb[i],MM))*g_yscale-g_yoff );
    /* now the known coordinates with potentials */
    fprintf(outf,"## known points: point number | x | y | potential\n");
    for ( i = unknowns; i < total; i++ )
      fprintf(outf,"%d "MDF" "MDF" "MDF"\n",i,  (Mx(re_renumb[i],MM))*g_xscale-g_xoff, (My(re_renumb[i],MM))*g_yscale-g_yoff, Mz(re_renumb[i],MM));
    /* now the triangles data with rho and epsilon */
    fprintf(outf,"## triangles : triangle number | point 1 | point 2 | point 3 | rho| mu or epsilon\n");


   DAG_traverse_list_reset(&dt);
   i=0;
   t=DAG_traverse_prune_list(&dt);
   while(t) {
    if ( t && t->status==LIVE ) {
     fprintf(outf,MIF" "MIF" "MIF" "MIF" ",i,renumb[t->p0],renumb[t->p1],renumb[t->p2]);
     fprintf(outf,MDF" ",(bound.poly_array[t->boundary].rhoexp?
       (ex(bound.poly_array[t->boundary].rhoexp,t->p0)+ex(bound.poly_array[t->boundary].rhoexp,t->p0)+ex(bound.poly_array[t->boundary].rhoexp,t->p0))*ONE_OVER_THREE:
       bound.poly_array[t->boundary].rho));
     fprintf(outf,MDF"\n",bound.poly_array[t->boundary].mu);
     i++;
    }
    t=DAG_traverse_prune_list(&dt);
   } 

  fclose(outf);
  free(renumb);
  free(re_renumb);

}

/* printing potentials as ACII file*/
void
print_potential(const char* filename, Mesh MM)
{
/* file format
  first line: no. of points
  other lines: point number, x, y, potential
*/

  FILE *outf;
  unsigned int i;
  if ((outf=fopen(filename,"w"))==NULL){
    fprintf(stderr,"%s: %d: could not open file %s\n",__FILE__,__LINE__,filename);
    return;	
   }
  
  fprintf(outf,"%d\n",M.count);
  for ( i = 0; i < M.count; i++ ) {
       fprintf(outf,"%d "MDF" "MDF" "MDF"\n",i, (Mx(i,MM))*g_xscale-g_xoff, (My(i,MM))*g_yscale-g_yoff, Mzz(i,current_plotting_contour,MM));
  }
  
  fclose(outf);
}


/* a function to plot the gradient of potential */
void
plot_gradient(Mesh MM)
{
  triangle *t;
  MY_INT n1,n2,n3;
  float r,g,bl;
  MY_DOUBLE del, a,b, yy12, yy23, yy31, x0,x1,x2,y0,y1,y2;

  max_gradient=0.01;
		DAG_traverse_list_reset(&dt);
		t=DAG_traverse_prune_list(&dt);
		while(t) {
    if ( t && t->status==LIVE ) {
#ifdef DEBUG
      printf("plotting gradient for "MIF" contour %d\n",t->n,current_plotting_contour);
#endif
      /* calculate gradient */
      /* phi =a*x+b*y+c for each triangle */
      /* need to find a and b only */
      /* using cramer's rule for this */
      yy23 = My(t->p1,MM)-My(t->p2,MM);
      yy31 = My(t->p2,MM)-My(t->p0,MM);
      yy12 = My(t->p0,MM)-My(t->p1,MM);
      del = ABS(Mx(t->p0,MM)*yy23 + Mx(t->p1,MM)*yy31 + Mx(t->p2,MM)*yy12);
      if ( del < TOL ) {
		t=DAG_traverse_prune_list(&dt);
 continue; /* do not plot */
   }
      /* else */
      del = 0.001/del;
      a = del*(Mzz(t->p0,current_plotting_contour,MM)*yy23 + Mzz(t->p1,current_plotting_contour,MM)*yy31 + Mzz(t->p2,current_plotting_contour,MM)*yy12);
      b = del*(Mx(t->p0,MM)*(Mzz(t->p1,current_plotting_contour,MM)-Mzz(t->p2,current_plotting_contour,MM))
        +Mx(t->p1,MM)*(Mzz(t->p2,current_plotting_contour,MM)-Mzz(t->p0,current_plotting_contour,MM))
        +Mx(t->p2,MM)*(Mzz(t->p0,current_plotting_contour,MM)-Mzz(t->p1,current_plotting_contour,MM)));

      /* now find the centroid of triangle */
      /* x coord */
      x0 = (Mx(t->p0,MM)+Mx(t->p1,MM)+Mx(t->p2,MM))*ONE_OVER_THREE;
      /* y coord */
      y0 = (My(t->p0,MM)+My(t->p1,MM)+My(t->p2,MM))*ONE_OVER_THREE;

      /* now find the two intersection points of the gradient line */
      /*               p2                                   */
      /*               /\                                   */
      /*            p1/  \                                  */
      /*             /`   \                                 */
      /*            /  `   \                                */
      /*           /    `   \                               */
      /*       p0 -------`----p1                            */
      /*                 p2                                 */

      /* at each corner, evaluate line value */
      /* check to see if n1 on line */
      yy12= a*(My(t->p0,MM)-y0) - b*(Mx(t->p0,MM)-x0);
      /* check to see if n2 on line */
      yy23= a*(My(t->p1,MM)-y0) - b*(Mx(t->p1,MM)-x0);
      /* check to see if n2 on line */
      yy31= a*(My(t->p2,MM)-y0) - b*(Mx(t->p2,MM)-x0);

      /* now 2 of the above values will have opposite signes from the other one */
      /* so generalize */
      if ( ((yy12 >= TOL) && (yy23 < 0) && (yy31 < 0)) || ((yy12 <= TOL) && (yy23 > 0) && (yy31 > 0))) {
 n1 = t->p0;
 n2 = t->p1;
 n3 = t->p2;
      } else if ( ((yy23 >= TOL) && (yy12 < 0) && (yy31 < 0)) || ((yy23 <= TOL) && (yy12 > 0) && (yy31 > 0))) {
 n1 = t->p1;
 n2 = t->p2;
 n3 = t->p0;
      } else if ( ((yy31 >= TOL) && (yy12 < 0) && (yy23 < 0)) || ((yy31 <= TOL) && (yy12 > 0) && (yy23 > 0)) ) {
 n1 = t->p2;
 n2 = t->p0;
 n3 = t->p1;
      } else {
		t=DAG_traverse_prune_list(&dt);
 continue; /* error */
  }
      /* now check intersection between (n1,n2) and (n1,n3) */
      /* first (n1,n2) */
      yy12= a*(My(n1,MM)-y0) - b*(Mx(n1,MM)-x0);
      yy23= a*(My(n2,MM)-y0) - b*(Mx(n2,MM)-x0);
      if ( (yy23 < TOL) && (yy23 > -TOL) ) {
 x1 = Mx(n2,MM);
 y1 = My(n2,MM);
      } else {
 del = -yy12/yy23;
 yy23 = 1/(1+del);
 x1 = (Mx(n1,MM)+del*Mx(n2,MM))*yy23;
 y1 = (My(n1,MM)+del*My(n2,MM))*yy23;
      }
      /* now (n1,n3) */
      yy23= a*(My(n3,MM)-y0) - b*(Mx(n3,MM)-x0);
      if ( (yy23 < TOL) && (yy23 > -TOL) ) {
 x2 = Mx(n3,MM);
 y2 = My(n3,MM);
      } else {
 del = -yy12/yy23;
 yy12 = 1/(1+del);
 x2 = (Mx(n1,MM)+del*Mx(n3,MM))*yy12;
 y2 = (My(n1,MM)+del*My(n3,MM))*yy12;
      }

  /* now we know the two endpoints */
  /* we draw an arrowhead in the direction point (x1,y1) */
  /* need to determine that now */
  if ( a > 0 ) {
      if ( x1 > x2 )
        n1 = 1; /* correct point */
          else
        n1 =2; /* need to switch */
  } else if ( a < 0 ) {
              if ( x1 < x2 )
        n1 = 1; /* correct point */
          else
        n1 =2; /* need to switch */
      } else { /* a == 0 */
          if ( b > 0 ) {
      if ( y1 > y2 )
        n1 = 1; /* correct point */
          else
        n1 =2; /* need to switch */
      } else if ( b < 0 ) {
              if ( y1 < y2 )
        n1 = 1; /* correct point */
          else
        n1 =2; /* need to switch */
          }
  }

  /* need to switch points? */
  if ( n1 == 2 ) {
      x0 = x1;
      y0 = y1;
      x1 = x2;
      y1 = y2;
      x2 = x0;
      y2 = y0;
    }


  /* arrow will look like this */
  /*                              (?  ,    ?)         */
  /*                                |\                */
  /*                                | \               */
  /*        (x2,y2)-----------------   \ (x1,y1)      */
  /*                         (x0,y0)| /               */
  /*                                |/                */
  /*                              (?  ?)              */

  /* keep 5:1 ratio */
  del=5;
  x0 = (x1*del+x2)/(del+1);
  y0 = (y1*del+y2)/(del+1);

  /* calculate value of gradient */
  a = sqrt(a*a + b*b);
  /* if this is higher than the maximum, update it */
  if ( max_gradient < a )
      max_gradient = a;
  /* set colour */
    get_colour(&r,&g,&bl,(int)(100*a), (int)(100*max_gradient));
      glColor3f(r, g, bl);
  /* arrow size */
  a = 0.2*sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
  /* gradient of normal */
  /* temporary variables */
  if ( y2 != y1 ) {
    del = -(x2-x1)/(y2-y1);
    yy12 = sqrt(1+del*del);
    yy23 = a/yy12;
            yy12 = del*a/yy12;
  } else {
      yy23 = 0;
      yy12 = a;
  }


      /* now plot gradient line (x1,y1) (x2,y2) */
      glBegin(GL_LINE_LOOP);
      /* draw line */
      glVertex2d( (GLfloat)(x2), (GLfloat)(y2) );
      glVertex2d( (GLfloat)(x1), (GLfloat)(y1) );
      glEnd();
     /* now arrow head*/
     glBegin(GL_TRIANGLES);
      glVertex2d( (GLfloat)(x0+yy23), (GLfloat)(y0+yy12) );
      glVertex2d( (GLfloat)(x0-yy23), (GLfloat)(y0-yy12) );
      glVertex2d( (GLfloat)(x1), (GLfloat)(y1) );
      glEnd();
  }
				
		t=DAG_traverse_prune_list(&dt);
  }
}
/**********************************************/
/* legend plotting routines */
#ifndef WIN32
static void
realize(GtkWidget *widget, gpointer data)
{
	GdkGLContext *glcontext =gtk_widget_get_gl_context(widget);
	GdkGLDrawable *gldrawable =gtk_widget_get_gl_drawable(widget);
	if (!gdk_gl_drawable_gl_begin(gldrawable,glcontext))
					return;
	gdk_gl_drawable_gl_end(gldrawable);
}

static gboolean
configure_event(GtkWidget *widget, GdkEventConfigure *event, gpointer data)
{
		GdkGLContext *glcontext=gtk_widget_get_gl_context(widget);
		GdkGLDrawable *gldrawable=gtk_widget_get_gl_drawable(widget);
		GLfloat w=widget->allocation.width;
		GLfloat h=widget->allocation.height;
		if (!gdk_gl_drawable_gl_begin(gldrawable,glcontext))
						return(FALSE);
		glViewport(0,0,w,h);
		glMatrixMode(GL_PROJECTION);
		glLoadIdentity();
		gluOrtho2D(-1,1,-1,1);
		glMatrixMode(GL_MODELVIEW);
		gdk_gl_drawable_gl_end(gldrawable);
		return(TRUE);
}

static gboolean
expose_event(GtkWidget *widget,
								GdkEventExpose *event,
								gpointer data)
{
  int i;
	float r,g,b;
	GLfloat dy; /* y step */
	GLfloat y1,y2,x1,x2,scale;
	GdkGLDrawable *gldrawable =gtk_widget_get_gl_drawable(widget);
       GdkGLContext *glcontext=gtk_widget_get_gl_context(widget);

  GLfloat w=(widget->allocation.width);
	GLfloat h=(widget->allocation.height);

if (!gdk_gl_drawable_gl_begin(gldrawable,glcontext))
   return(FALSE);

	/* rescale */
 if ( h >w )
		scale=1/h;
 else 
	  scale=1/w;
		w=2;
		h=2;

/* top, bottom border =h/20 */
 /* left border = w/10, color block=9w/10, right border=w/10 */
	/* drawing */
	 dy=(h*0.9)/(GLfloat)contour_levels;

  y2=-0.45*h;
	y1=dy+y2;
  x1=-0.45*w;
	x2=-x1;
   glClear(GL_COLOR_BUFFER_BIT);
#ifdef DEBUG
	printf("expose_ev: legend: rectangle (%f,%f)--(%f,%f)\n",x1,y1,x2,y2);
#endif


	/* horizontal line */
	glColor3f(0.0f,0.0f,0.0f);
	glBegin(GL_LINES);
	glVertex2d(x1,y2);
	glVertex2d(x1+0.9*w,y2);
	glEnd();

  for (i=0; i<contour_levels; i++) {
				get_colour(&r,&g,&b,i,contour_levels);
				glColor3f((GLfloat)r,(GLfloat)g,(GLfloat)b);
				glBegin(GL_QUADS);
				glVertex2d(x1,y1);
				glVertex2d(x2,y1);
				glVertex2d(x2,y2);
				glVertex2d(x1,y2);
				glEnd();

	/* horizontal line */
	glColor3f(0.0f,0.7f,0.7f);
	glBegin(GL_LINES);
	glVertex2d(x1,y1);
	glVertex2d(x1+0.0*w,y1);
	glEnd();

	y1+=dy;
	y2+=dy;

	}
if ( gdk_gl_drawable_is_double_buffered(gldrawable))
		gdk_gl_drawable_swap_buffers(gldrawable);
else
		glFlush();
gdk_gl_drawable_gl_end(gldrawable);
 return(TRUE);
}

/* mouse motion */
static gboolean
motion_notify_event(GtkWidget *widget,
   GdkEventMotion *event,
   gpointer data) 
{
 GtkWidget *label;
 GLfloat h,y;
 int i;
 gchar buf[128];
 h=widget->allocation.height;

 y=-2*event->y/h+1.0;
 label=(GtkWidget*)data;
 /* determine exact contour we are in */
 i=(int)floor((double)(y+0.9)*(double)contour_levels/1.8);
 /* sanity check: i has to be in [0,contour_levels-1]*/
 if ( (i>= 0) && (i<contour_levels)) {
 g_snprintf(buf,128,"%4.2f",(float)i*(g_maxpot-g_minpot)/(float)contour_levels+g_minpot);
/* printf("%4.3f\n",(float)i*(g_maxpot-g_minpot)/(float)contour_levels+g_minpot); */
 gtk_label_set_markup(GTK_LABEL(label),buf);
 return(TRUE);
 }
 return(FALSE);
}

static gint
close_legend_window(GtkWidget *widget, GdkEvent *event, gpointer *data)
{
  return(FALSE);
}

/* plotting legend */
void
display_legend_window(void)
{
  GtkWidget *window;
  GtkWidget *vbox;
  GtkWidget *drawing_area;
  GtkWidget *legend_label;
	GdkGLConfig *glconfig;


  /* Try double-buffered visual */
  glconfig = gdk_gl_config_new_by_mode (GDK_GL_MODE_RGB |
					GDK_GL_MODE_DOUBLE);
  if (glconfig == NULL)
    {
      g_print ("\n*** Cannot find the double-buffered visual.\n");
      g_print ("\n*** Trying single-buffered visual.\n");

      /* Try single-buffered visual */
      glconfig = gdk_gl_config_new_by_mode (GDK_GL_MODE_RGB);
      if (glconfig == NULL)
	{
	  g_print ("*** No appropriate OpenGL-capable visual found.\n");
	  exit(1);
	}
    }


  /*
   * Top-level window.
   */

  window = gtk_window_new (GTK_WINDOW_TOPLEVEL);
  gtk_window_set_title (GTK_WINDOW (window), "Legend");

  /* Get automatically redrawn if any of their children changed allocation. */
  gtk_container_set_reallocate_redraws (GTK_CONTAINER (window), TRUE);

  /* Connect signal handlers to the window */
  g_signal_connect (G_OBJECT (window), "destroy",
		    G_CALLBACK (gtk_widget_destroy), (gpointer)window);
  gtk_signal_connect (GTK_OBJECT (window), "destroy",
              GTK_SIGNAL_FUNC (close_legend_window), NULL);

  /*
   * VBox.
   */

  vbox = gtk_vbox_new (FALSE, 0);
  gtk_container_add (GTK_CONTAINER (window), vbox);
  gtk_widget_show (vbox);

  /*
   * Drawing area to draw OpenGL scene.
   */

  drawing_area = gtk_drawing_area_new ();
  gtk_widget_set_size_request (drawing_area, LEGEND_WIDTH, LEGEND_HEIGHT);

  /* Set OpenGL-capability to the widget */
  gtk_widget_set_gl_capability (drawing_area,
				glconfig,
				NULL,
				TRUE,
				GDK_GL_RGBA_TYPE);

  gtk_widget_add_events (drawing_area,
			 GDK_POINTER_MOTION_MASK    |  /* any pointer motion */
			 GDK_VISIBILITY_NOTIFY_MASK);

  /* Connect signal handlers to the drawing area */
  g_signal_connect_after (G_OBJECT (drawing_area), "realize",
                          G_CALLBACK (realize), NULL);
  g_signal_connect (G_OBJECT (drawing_area), "configure_event",
		    G_CALLBACK (configure_event), NULL);
  g_signal_connect (G_OBJECT (drawing_area), "expose_event",
		    G_CALLBACK (expose_event), NULL);

  gtk_box_pack_start (GTK_BOX (vbox), drawing_area, TRUE, TRUE, 0);
  gtk_widget_show (drawing_area);

  /* label to show potential */
  legend_label=gtk_label_new("000");
  gtk_label_set_justify(GTK_LABEL(legend_label),GTK_JUSTIFY_CENTER);
  gtk_box_pack_start(GTK_BOX(vbox),legend_label,FALSE,TRUE,0);
  gtk_widget_show(legend_label);

  /* connect mouse motion */
  g_signal_connect (G_OBJECT(drawing_area),"motion_notify_event",
        G_CALLBACK(motion_notify_event), (gpointer)legend_label);


  gtk_widget_show(window);

}
#endif/*WIN32*/

/* a function to print the gradient of potential */
void
print_gradient(const char *filename, Mesh MM)
{
  triangle *t;
  MY_INT n1,n2,n3;
  MY_DOUBLE del, a,b, yy12, yy23, yy31, x0,x1,x2,y0,y1,y2;
	FILE *plotf;
  MY_DOUBLE myscale;
  int i;

	if ((plotf=fopen(filename,"w"))==NULL){
			fprintf(stderr,"%s: %d: could not open file %s\n",__FILE__,__LINE__,filename);
	   return;	
	}


  /* determine the right bound of bounding box */
  if ( g_xrange > g_yrange )
    myscale = 400.0/g_xrange;
  else
    myscale = 400.0/g_yrange;

  /* print the header */
  fprintf(plotf,"%%!PS-Adobe-2.0 EPSF-3.0\n");
  fprintf(plotf,"%%%%Title: %s\n",filename);
  fprintf(plotf,"%%%%Generated By: "PACKAGE" "VERSION"  contour levels: %d, min= "MDF", max= "MDF"\n", contour_levels, g_minpot, g_maxpot);
  fprintf(plotf,"%%%%BoundingBox: 0 0 %4.2lf %4.2lf\n",(double)(2*g_xrange*myscale+20), (double)(2*g_yrange*myscale+20) );
  fprintf(plotf,"%%%%Magnification: 1.0000\n");
  
  

  /* start plotting from 10 10 */
  fprintf(plotf,"10 10 translate\n");

  /* now print some shorter definitions */
  fprintf(plotf,"/l {lineto} bind def\n"); /* lineto */
  fprintf(plotf,"/f {fill} bind def\n"); /* lineto */
  fprintf(plotf,"/m {moveto} bind def\n"); /* moveto */
  fprintf(plotf,"/n {newpath} bind def\n"); /* newpath */
  fprintf(plotf,"/s {stroke} bind def\n"); /* stroke */
  fprintf(plotf,"/slw {setlinewidth} bind def\n"); /* setlinewidth */
  fprintf(plotf,"/cp {closepath} bind def\n"); /* closepath */

  fprintf(plotf,"1 slw\n");
		DAG_traverse_list_reset(&dt);
		t=DAG_traverse_prune_list(&dt);
		while(t) {
    if ( t && t->status==LIVE ) {
#ifdef DEBUG
      printf("plotting gradient for "MIF" contour %d\n",t->n,current_plotting_contour);
#endif
      /* calculate gradient */
      /* phi =a*x+b*y+c for each triangle */
      /* need to find a and b only */
      /* using cramer's rule for this */
      yy23 = My(t->p1,MM)-My(t->p2,MM);
      yy31 = My(t->p2,MM)-My(t->p0,MM);
      yy12 = My(t->p0,MM)-My(t->p1,MM);
      del = ABS(Mx(t->p0,MM)*yy23 + Mx(t->p1,MM)*yy31 + Mx(t->p2,MM)*yy12);
      if ( del < TOL ) {
		t=DAG_traverse_prune_list(&dt);
 continue; /* do not plot */
   }
      /* else */
      del = 0.001/del;
      a = del*(Mzz(t->p0,current_plotting_contour,MM)*yy23 + Mzz(t->p1,current_plotting_contour,MM)*yy31 + Mzz(t->p2,current_plotting_contour,MM)*yy12);
      b = del*(Mx(t->p0,MM)*(Mzz(t->p1,current_plotting_contour,MM)-Mzz(t->p2,current_plotting_contour,MM))
        +Mx(t->p1,MM)*(Mzz(t->p2,current_plotting_contour,MM)-Mzz(t->p0,current_plotting_contour,MM))
        +Mx(t->p2,MM)*(Mzz(t->p0,current_plotting_contour,MM)-Mzz(t->p1,current_plotting_contour,MM)));

      /* now find the centroid of triangle */
      /* x coord */
      x0 = (Mx(t->p0,MM)+Mx(t->p1,MM)+Mx(t->p2,MM))*ONE_OVER_THREE;
      /* y coord */
      y0 = (My(t->p0,MM)+My(t->p1,MM)+My(t->p2,MM))*ONE_OVER_THREE;

      /* now find the two intersection points of the gradient line */
      /*               p2                                   */
      /*               /\                                   */
      /*            p1/  \                                  */
      /*             /`   \                                 */
      /*            /  `   \                                */
      /*           /    `   \                               */
      /*       p0 -------`----p1                            */
      /*                 p2                                 */

      /* at each corner, evaluate line value */
      /* check to see if n1 on line */
      yy12= a*(My(t->p0,MM)-y0) - b*(Mx(t->p0,MM)-x0);
      /* check to see if n2 on line */
      yy23= a*(My(t->p1,MM)-y0) - b*(Mx(t->p1,MM)-x0);
      /* check to see if n2 on line */
      yy31= a*(My(t->p2,MM)-y0) - b*(Mx(t->p2,MM)-x0);

      /* now 2 of the above values will have opposite signes from the other one */
      /* so generalize */
      if ( ((yy12 >= TOL) && (yy23 < 0) && (yy31 < 0)) || ((yy12 <= TOL) && (yy23 > 0) && (yy31 > 0))) {
 n1 = t->p0;
 n2 = t->p1;
 n3 = t->p2;
      } else if ( ((yy23 >= TOL) && (yy12 < 0) && (yy31 < 0)) || ((yy23 <= TOL) && (yy12 > 0) && (yy31 > 0))) {
 n1 = t->p1;
 n2 = t->p2;
 n3 = t->p0;
      } else if ( ((yy31 >= TOL) && (yy12 < 0) && (yy23 < 0)) || ((yy31 <= TOL) && (yy12 > 0) && (yy23 > 0)) ) {
 n1 = t->p2;
 n2 = t->p0;
 n3 = t->p1;
      } else {
		t=DAG_traverse_prune_list(&dt);
 continue; /* error */
  }
      /* now check intersection between (n1,n2) and (n1,n3) */
      /* first (n1,n2) */
      yy12= a*(My(n1,MM)-y0) - b*(Mx(n1,MM)-x0);
      yy23= a*(My(n2,MM)-y0) - b*(Mx(n2,MM)-x0);
      if ( (yy23 < TOL) && (yy23 > -TOL) ) {
 x1 = Mx(n2,MM);
 y1 = My(n2,MM);
      } else {
 del = -yy12/yy23;
 yy23 = 1/(1+del);
 x1 = (Mx(n1,MM)+del*Mx(n2,MM))*yy23;
 y1 = (My(n1,MM)+del*My(n2,MM))*yy23;
      }
      /* now (n1,n3) */
      yy23= a*(My(n3,MM)-y0) - b*(Mx(n3,MM)-x0);
      if ( (yy23 < TOL) && (yy23 > -TOL) ) {
 x2 = Mx(n3,MM);
 y2 = My(n3,MM);
      } else {
 del = -yy12/yy23;
 yy12 = 1/(1+del);
 x2 = (Mx(n1,MM)+del*Mx(n3,MM))*yy12;
 y2 = (My(n1,MM)+del*My(n3,MM))*yy12;
      }

  /* now we know the two endpoints */
  /* we draw an arrowhead in the direction point (x1,y1) */
  /* need to determine that now */
  if ( a > 0 ) {
      if ( x1 > x2 )
        n1 = 1; /* correct point */
          else
        n1 =2; /* need to switch */
  } else if ( a < 0 ) {
              if ( x1 < x2 )
        n1 = 1; /* correct point */
          else
        n1 =2; /* need to switch */
      } else { /* a == 0 */
          if ( b > 0 ) {
      if ( y1 > y2 )
        n1 = 1; /* correct point */
          else
        n1 =2; /* need to switch */
      } else if ( b < 0 ) {
              if ( y1 < y2 )
        n1 = 1; /* correct point */
          else
        n1 =2; /* need to switch */
          }
  }

  /* need to switch points? */
  if ( n1 == 2 ) {
      x0 = x1;
      y0 = y1;
      x1 = x2;
      y1 = y2;
      x2 = x0;
      y2 = y0;
    }


  /* arrow will look like this */
  /*                              (?  ,    ?)         */
  /*                                |\                */
  /*                                | \               */
  /*        (x2,y2)-----------------   \ (x1,y1)      */
  /*                         (x0,y0)| /               */
  /*                                |/                */
  /*                              (?  ?)              */

  /* keep 5:1 ratio */
  del=5;
  x0 = (x1*del+x2)/(del+1);
  y0 = (y1*del+y2)/(del+1);

  /* calculate value of gradient */
  a = sqrtf((float)(a*a + b*b));
  /* arrow size */
  a = 0.2*sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
  /* gradient of normal */
  /* temporary variables */
  if ( y2 != y1 ) {
    del = -(x2-x1)/(y2-y1);
    yy12 = sqrt(1+del*del);
    yy23 = a/yy12;
            yy12 = del*a/yy12;
  } else {
      yy23 = 0;
      yy12 = a;
  }


      /* now plot gradient line (x1,y1) (x2,y2) */
    fprintf(plotf,"n ");
    fprintf(plotf,"%lf %lf m ",(double) (g_xrange+x2*g_xscale)*myscale, (double)(g_yrange + y2*g_yscale)*myscale);
    fprintf(plotf,"%lf %lf l ",(double) (g_xrange+x0*g_xscale)*myscale, (double)(g_yrange + y0*g_yscale)*myscale);
    fprintf(plotf,"cp s\n");
     /* now arrow head*/
    fprintf(plotf,"n ");
    fprintf(plotf,"%lf %lf m ",(double) (g_xrange+(x0+yy23)*g_xscale)*myscale, (double)(g_yrange +(y0+yy12)*g_yscale)*myscale);
    fprintf(plotf,"%lf %lf l ",(double) (g_xrange+(x0-yy23)*g_xscale)*myscale, (double)(g_yrange +(y0-yy12)*g_yscale)*myscale);
    fprintf(plotf,"%lf %lf l ",(double) (g_xrange+x1*g_xscale)*myscale, (double)(g_yrange + y1*g_yscale)*myscale);
    fprintf(plotf,"cp s\n");
  }
				
		t=DAG_traverse_prune_list(&dt);
  }

  fprintf(plotf,"2 slw\n");
/* print boundaries */
  for ( i = 0; i < bound.npoly; i++ ) {
		print_polygon(&bound.poly_array[i],myscale,plotf,MM);
	}

  fprintf(plotf,"%%%%EOF\n");
  fclose(plotf);

}