Index: /tags/PIGC_V1/SophyaPI/PIGcont/alloc.h
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/alloc.h	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/alloc.h	(revision 1830)
@@ -0,0 +1,92 @@
+/*
+ * $Id: alloc.h,v 1.1.1.1 2001-12-19 18:01:18 perderos Exp $
+ */
+
+/* GNUPLOT - alloc.h */
+
+/*[
+ * Copyright 1986 - 1993, 1998   Thomas Williams, Colin Kelley
+ *
+ * Permission to use, copy, and distribute this software and its
+ * documentation for any purpose with or 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.
+ *
+ * Permission to modify the software is granted, but not the right to
+ * distribute the complete modified source code.  Modifications are to
+ * be distributed as patches to the released version.  Permission to
+ * distribute binaries produced by compiling modified sources is granted,
+ * provided you
+ *   1. distribute the corresponding source modifications from the
+ *    released version in the form of a patch file along with the binaries,
+ *   2. add special version identification to distinguish your version
+ *    in addition to the base release version number,
+ *   3. provide your name and address as the primary contact for the
+ *    support of your modified version, and
+ *   4. retain our contact information in regard to use of the base
+ *    software.
+ * Permission to distribute the released version of the source code along
+ * with corresponding source modifications in the form of a patch file is
+ * granted with same provisions 2 through 4 for binary distributions.
+ *
+ * This software is provided "as is" without express or implied warranty
+ * to the extent permitted by applicable law.
+]*/
+
+#ifndef GNUPLOT_ALLOC_H
+# define GNUPLOT_ALLOC_H
+
+#include <stdio.h>
+#include <stdlib.h>
+#include "graph3d.h"
+/****
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include "stdfn.h"
+****/
+/* prototypes from "alloc.c". This file figures out if the free hack is needed
+ * and redefines free if necessary.
+ */
+ 
+void *gp_alloc(size_t size, const char *message);
+void *gp_realloc(void *p, size_t size, const char *message);
+struct iso_curve *iso_alloc(int);
+
+void iso_free(struct iso_curve *);
+void gp_free(void *);
+void contour_free(struct gnuplot_contours *ptr);
+/* dont define CHECK_HEAP_USE on a FARALLOC machine ! */
+
+#ifdef CHECK_HEAP_USE
+
+/* all allocated blocks have guards at front and back.
+ * CHECK_POINTER checks guards on block, and checks that p is in range
+ * START_LEAK_CHECK and END_LEAK_CHECK allow assert that no net memory
+ * is allocated within enclosed block
+ */
+ 
+void checked_free(void *p);
+void check_pointer_in_block(void *block, void *p, int size, char *file, int line);
+void start_leak_check(char *file,int line);
+void end_leak_check(char *file,int line);
+/*# define free(x) checked_free(x)*/
+# define CHECK_POINTER(block, p) check_pointer_in_block(block, p, sizeof(*p), __FILE__, __LINE__)
+# define START_LEAK_CHECK() start_leak_check(__FILE__, __LINE__)
+# define END_LEAK_CHECK() end_leak_check(__FILE__, __LINE__)
+#else
+# define CHECK_POINTER(block, p) /*nowt*/
+# define START_LEAK_CHECK() /*nowt*/
+# define END_LEAK_CHECK() /*nowt*/
+#endif
+
+/***
+#if defined(MSDOS) && defined(__TURBOC__) && !defined(DOSX286) || defined(_Windows) && !defined(WIN32) 
+#define FARALLOC 
+void gpfree((generic *p)); 
+#define free gpfree 
+#endif 
+***/
+#endif /* GNUPLOT_ALLOC_H */
Index: /tags/PIGC_V1/SophyaPI/PIGcont/contour.c
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/contour.c	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/contour.c	(revision 1830)
@@ -0,0 +1,1534 @@
+#ifndef lint
+static char *RCSid() { return RCSid("$Id: contour.c,v 1.1.1.1 2001-12-19 18:01:18 perderos Exp $"); }
+#endif
+
+/* GNUPLOT - contour.c */
+
+/*[
+ * Copyright 1986 - 1993, 1998   Thomas Williams, Colin Kelley
+ *
+ * Permission to use, copy, and distribute this software and its
+ * documentation for any purpose with or 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.
+ *
+ * Permission to modify the software is granted, but not the right to
+ * distribute the complete modified source code.  Modifications are to
+ * be distributed as patches to the released version.  Permission to
+ * distribute binaries produced by compiling modified sources is granted,
+ * provided you
+ *   1. distribute the corresponding source modifications from the
+ *    released version in the form of a patch file along with the binaries,
+ *   2. add special version identification to distinguish your version
+ *    in addition to the base release version number,
+ *   3. provide your name and address as the primary contact for the
+ *    support of your modified version, and
+ *   4. retain our contact information in regard to use of the base
+ *    software.
+ * Permission to distribute the released version of the source code along
+ * with corresponding source modifications in the form of a patch file is
+ * granted with same provisions 2 through 4 for binary distributions.
+ *
+ * This software is provided "as is" without express or implied warranty
+ * to the extent permitted by applicable law.
+]*/
+
+ 
+/* 
+ * AUTHORS
+ * 
+ *   Original Software:
+ *       Gershon Elber
+ * 
+ *   Improvements to the numerical algorithms:
+ *        Hans-Martin Keller, 1995,1997 (hkeller@gwdg.de)
+ *
+ */
+
+#include "contour.h"
+
+#include "alloc.h"
+#include "myaxis.h"
+/*  #include "setshow.h" */
+
+/* exported variables (to be handled by the 'set' and friends): */
+
+char contour_format[32] = "%8.3g";	/* format for contour key entries */
+t_contour_kind contour_kind = CONTOUR_KIND_LINEAR;
+t_contour_levels_kind contour_levels_kind = LEVELS_AUTO;
+int contour_levels = DEFAULT_CONTOUR_LEVELS;
+int contour_order = DEFAULT_CONTOUR_ORDER;
+int contour_pts = DEFAULT_NUM_APPROX_PTS;
+
+dynarray dyn_contour_levels_list;/* storage for z levels to draw contours at */
+double * contour_levels_list=NULL;
+
+/* position of edge in mesh */
+typedef enum en_edge_position {
+    INNER_MESH=1,
+    BOUNDARY,
+    DIAGONAL
+} t_edge_position;
+
+
+/* FIXME HBB 2000052: yet another local copy of 'epsilon'. Why? */
+#define EPSILON  1e-5		/* Used to decide if two float are equal. */
+
+#ifndef TRUE
+#define TRUE     -1
+#define FALSE    0
+#endif
+
+#define MAX_POINTS_PER_CNTR 	100
+
+#define SQR(x)  ((x) * (x))
+
+/*
+ * struct vrtx_struct {
+ *      double X, Y, Z;
+ *      struct vrtx_struct *next;
+ * };
+ * 
+ * replaced by 'struct coordinate  ', see plot.h (HMK 1997) 
+ */
+
+struct edge_struct {
+    struct poly_struct *poly[2]; /* Each edge belongs to up to 2 polygons */
+    struct coordinate  *vertex[2]; /* The two extreme points of this edge. */
+    struct edge_struct *next;	/* To chain lists */
+    TBOOLEAN is_active;		/* is edge is 'active' at certain Z level? */
+    t_edge_position position;	/* position of edge in mesh */
+};
+
+struct poly_struct {
+    struct edge_struct *edge[3];	/* As we do triangolation here... */
+    struct poly_struct *next;	/* To chain lists. */
+};
+
+struct cntr_struct {		/* Contours are saved using this struct list. */
+    double X, Y;		/* The coordinates of this vertex. */
+    struct cntr_struct *next;	/* To chain lists. */
+};
+
+static struct gnuplot_contours *contour_list = NULL;
+static double crnt_cntr[MAX_POINTS_PER_CNTR * 2];
+static int crnt_cntr_pt_index = 0;
+static double contour_level = 0.0;
+
+/* Linear, Cubic interp., Bspline: */
+static t_contour_kind interp_kind = CONTOUR_KIND_LINEAR;
+
+static double x_min, y_min, z_min;	/* Minimum values of x, y, and z */
+static double x_max, y_max, z_max;	/* Maximum values of x, y, and z */
+
+static void add_cntr_point (double x, double y);
+static void end_crnt_cntr (void);
+static void gen_contours(struct edge_struct * p_edges, double z_level,
+				  double xx_min, double xx_max, double yy_min, double yy_max);
+static int update_all_edges(struct edge_struct * p_edges,
+				     double z_level);
+static struct cntr_struct *gen_one_contour (
+						       struct edge_struct * p_edges, double
+						       z_level, TBOOLEAN *contr_isclosed,
+						       int *num_active);
+static struct cntr_struct *trace_contour (
+						     struct edge_struct * pe_start, double
+						     z_level, int *num_active,
+						     TBOOLEAN contr_isclosed);
+static struct cntr_struct *update_cntr_pt (struct edge_struct * p_edge,
+						   double z_level);
+static int fuzzy_equal (struct cntr_struct * p_cntr1,
+				struct cntr_struct * p_cntr2);
+
+
+static void gen_triangle (int num_isolines,
+				  struct iso_curve * iso_lines, struct poly_struct ** p_polys,
+				  struct edge_struct ** p_edges);
+static void calc_min_max (int num_isolines,
+				  struct iso_curve * iso_lines, double *xx_min, double *yy_min,
+				  double *zz_min,
+				  double *xx_max, double *yy_max, double *zz_max);
+static struct edge_struct *add_edge (struct coordinate  * point0,
+					     struct coordinate  * point1, struct edge_struct
+					     ** p_edge,
+					     struct edge_struct ** pe_tail);
+static struct poly_struct *add_poly (struct edge_struct * edge0,
+					     struct edge_struct * edge1, struct edge_struct * edge2,
+					     struct poly_struct ** p_poly, struct poly_struct ** pp_tail);
+
+
+static void put_contour (struct cntr_struct * p_cntr, double z_level,
+				 double xx_min, double xx_max, double yy_min, double yy_max,
+		 		 TBOOLEAN contr_isclosed);
+static void put_contour_nothing (struct cntr_struct * p_cntr);
+static int chk_contour_kind (struct cntr_struct * p_cntr,
+				     TBOOLEAN contr_isclosed);
+static void put_contour_cubic (struct cntr_struct * p_cntr,
+				       double z_level, double xx_min, double xx_max, double
+				       yy_min, double yy_max,
+				       TBOOLEAN contr_isclosed);
+static void put_contour_bspline (struct cntr_struct * p_cntr,
+					 double z_level, double xx_min, double xx_max, double
+					 yy_min, double yy_max,
+					 TBOOLEAN contr_isclosed);
+static void free_contour (struct cntr_struct * p_cntr);
+static int count_contour (struct cntr_struct * p_cntr);
+static int gen_cubic_spline (int num_pts, struct cntr_struct * p_cntr,
+				     double d2x[], double d2y[], double delta_t[], TBOOLEAN contr_isclosed,
+				     double unit_x, double unit_y);
+static void intp_cubic_spline (int n, struct cntr_struct * p_cntr,
+				       double d2x[], double d2y[], double delta_t[], int n_intpol);
+static int solve_cubic_1 (tri_diag m[], int n);
+static void solve_cubic_2 (tri_diag m[], double x[], int n);
+static void gen_bspline_approx (struct cntr_struct * p_cntr,
+					int num_of_points, int order, TBOOLEAN contr_isclosed);
+static void eval_bspline (double t, struct cntr_struct * p_cntr,
+				  int num_of_points, int order, int j, TBOOLEAN contr_isclosed, double *x,
+				  double *y);
+static double fetch_knot (TBOOLEAN contr_isclosed, int num_of_points,
+				  int order, int i);
+
+
+static int num_of_z_levels;/*_____ OP ___________*/	/* # Z contour levels. */
+
+int Get_Num_Of_Z_Levels(){
+return num_of_z_levels;
+} /* OP       __________ */
+
+/*
+ * Entry routine to this whole set of contouring module.
+ */
+struct gnuplot_contours *
+contour(num_isolines, iso_lines)
+int num_isolines;
+struct iso_curve *iso_lines;
+{
+    int i;
+    /*OP    int num_of_z_levels;*/	/* # Z contour levels. */
+    struct poly_struct *p_polys, *p_poly;
+    struct edge_struct *p_edges, *p_edge;
+    double z = 0, dz = 0;
+    struct gnuplot_contours *save_contour_list;
+
+    num_of_z_levels = contour_levels;
+    interp_kind = contour_kind;
+
+    contour_list = NULL;
+
+    /* 
+     * Calculate min/max values :
+     */
+    calc_min_max(num_isolines, iso_lines,
+		 &x_min, &y_min, &z_min, &x_max, &y_max, &z_max);
+
+    dz = fabs(z_max - z_min);
+    /* 
+     *  printf(" contour z_max %g z_min %g dz=%g kind %d   \n",z_max, z_min , dz,contour_levels_kind);
+     * Generate list of edges (p_edges) and list of triangles (p_polys):
+     */
+
+    gen_triangle(num_isolines, iso_lines, &p_polys, &p_edges);
+    crnt_cntr_pt_index = 0;
+    /*AJOUT OP */
+    if (contour_levels_kind == LEVELS_NUM) {
+	dz = fabs(z_max - z_min)/(num_of_z_levels);
+	z = z_min - dz/2.;
+    }
+
+    if (contour_levels_kind == LEVELS_AUTO) {
+	dz = fabs(z_max - z_min);
+	/*printf(" contour z_max %g z_min %g dz=%g\n",z_max, z_min , dz);*/
+	if (dz == 0)
+	    return NULL;	/* empty z range ? */
+	/* what is the deeper sense of this ? (joze) */
+	dz = set_tic(log10(dz), ((int) contour_levels + 1) * 2);
+	z = floor(z_min / dz) * dz;
+	num_of_z_levels = (int) floor((z_max - z) / dz);
+        /*printf("contour() : num_of_z_levels %d\n",num_of_z_levels);*/
+    }
+    for (i = 0; i < num_of_z_levels; i++) {
+	switch (contour_levels_kind) {
+	case LEVELS_AUTO:
+	case LEVELS_NUM:
+	    z += dz;
+	    break;
+	case LEVELS_INCREMENTAL:
+	    z = contour_levels_list[0] + i * contour_levels_list[1];
+	    break;
+	case LEVELS_DISCRETE:
+	  /*printf("????? contour z=%f\n",contour_levels_list[i] );*/
+	    /*z = AXIS_LOG_VALUE(FIRST_Z_AXIS, contour_levels_list[i]); PAS BESOIN ? OP */
+	    z = contour_levels_list[i];
+	    break;
+	}
+	contour_level = z;
+	/*printf(" contour z=%f\n",z);*/
+	save_contour_list = contour_list;
+	gen_contours(p_edges, z, x_min, x_max, y_min, y_max);
+	if (contour_list != save_contour_list) {
+	    contour_list->isNewLevel = 1;
+	    sprintf(contour_list->label, contour_format, AXIS_DE_LOG_VALUE(FIRST_Z_AXIS,z));
+#ifdef PM3D
+	    contour_list->z = AXIS_DE_LOG_VALUE(FIRST_Z_AXIS, z);
+#endif
+	}
+    }
+
+    /* Free all contouring related temporary data. */ 
+    while (p_polys) {
+	p_poly = p_polys->next;
+	free(p_polys);
+	p_polys = p_poly; 
+    }
+    while (p_edges) {
+	p_edge = p_edges->next;
+	free(p_edges);
+	p_edges = p_edge;
+    }
+
+    return contour_list;
+}
+
+/*
+ * Adds another point to the currently build contour.
+ */
+static void
+add_cntr_point(x, y)
+double x, y;
+{
+    int index;
+
+    if (crnt_cntr_pt_index >= MAX_POINTS_PER_CNTR - 1) {
+	index = crnt_cntr_pt_index - 1;
+	end_crnt_cntr();
+	crnt_cntr[0] = crnt_cntr[index * 2];
+	crnt_cntr[1] = crnt_cntr[index * 2 + 1];
+	crnt_cntr_pt_index = 1;	/* Keep the last point as first of this one. */
+    }
+    crnt_cntr[crnt_cntr_pt_index * 2] = x;
+    crnt_cntr[crnt_cntr_pt_index * 2 + 1] = y;
+    crnt_cntr_pt_index++;
+}
+ 
+/*
+ * Done with current contour - create gnuplot data structure for it.
+ */
+static void
+end_crnt_cntr()
+{
+    int i;
+    struct gnuplot_contours *cntr = (struct gnuplot_contours *)
+    gp_alloc(sizeof(struct gnuplot_contours), "gnuplot_contour");
+    cntr->coords = (struct coordinate  *)
+	gp_alloc(sizeof(struct coordinate) * crnt_cntr_pt_index,
+		 "contour coords");
+
+    for (i = 0; i < crnt_cntr_pt_index; i++) {
+	cntr->coords[i].x = crnt_cntr[i * 2];
+	cntr->coords[i].y = crnt_cntr[i * 2 + 1];
+	cntr->coords[i].z = contour_level;
+    }
+    cntr->num_pts = crnt_cntr_pt_index;
+
+    cntr->next = contour_list;
+    contour_list = cntr;
+    contour_list->isNewLevel = 0;
+
+    crnt_cntr_pt_index = 0;
+}
+
+/*
+ * Generates all contours by tracing the intersecting triangles.
+ */
+static void
+gen_contours(p_edges, z_level, xx_min, xx_max, yy_min, yy_max)
+struct edge_struct *p_edges;
+double z_level, xx_min, xx_max, yy_min, yy_max;
+{
+    int num_active;		/* Number of edges marked ACTIVE. */
+    TBOOLEAN contr_isclosed;	/* Is this contour a closed line? */
+    struct cntr_struct *p_cntr;
+
+    num_active = update_all_edges(p_edges, z_level);	/* Do pass 1. */
+
+    contr_isclosed = FALSE;	/* Start to look for contour on boundaries. */
+    /*printf("<gen_contour> z=%g num_active %d \n",z_level,num_active);*/
+    while (num_active > 0) {	/* Do Pass 2. */
+	/* Generate One contour (and update MumActive as needed): */
+	p_cntr = gen_one_contour(p_edges, z_level, &contr_isclosed, &num_active);
+	if (p_cntr ==NULL) printf("<gen_contour> gen_one_contour retourne NULL \n");
+	/* Emit it in requested format: */
+	put_contour(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_isclosed);
+    }
+}
+
+/*
+ * Does pass 1, or marks the edges which are active (crosses this z_level)
+ * Returns number of active edges (marked ACTIVE).
+ */
+static int
+update_all_edges(p_edges, z_level)
+struct edge_struct *p_edges;
+double z_level;
+{
+    int count = 0;
+
+    while (p_edges) {
+	/* use the same test at both vertices to avoid roundoff errors */
+        
+	if ((p_edges->vertex[0]->z >= z_level) !=
+	    (p_edges->vertex[1]->z >= z_level)) {
+	    p_edges->is_active = TRUE;
+	    count++;
+	} else
+	    p_edges->is_active = FALSE;
+	p_edges = p_edges->next;
+    }
+
+    return count;
+}
+
+/*
+ * Does pass 2, or find one complete contour out of the triangulation
+ * data base:
+ *
+ * Returns a pointer to the contour (as linked list), contr_isclosed
+ * tells if the contour is a closed line or not, and num_active is
+ * updated.  
+ */
+static struct cntr_struct *
+gen_one_contour(p_edges, z_level, contr_isclosed, num_active)
+struct edge_struct *p_edges;	/* list of edges input */
+double z_level;			/* Z level of contour input */
+TBOOLEAN *contr_isclosed;	/* open or closed contour, in/out */
+int *num_active;		/* number of active edges     in/out */
+{
+    struct edge_struct *pe_temp;
+    
+    if (! *contr_isclosed) {
+      /*printf("<gen_one_contour> contr_isclosed FALSE \n");*/
+	/* Look for something to start with on boundary: */
+	pe_temp = p_edges;
+	while (pe_temp) {
+	    if (pe_temp->is_active && (pe_temp->position == BOUNDARY))
+		break;
+	    pe_temp = pe_temp->next;
+	}
+	if (!pe_temp)
+	    *contr_isclosed = TRUE;	/* No more contours on boundary. */
+	else {
+	  /*printf("<gen_one_contour> contr_isclosed FALSE return \n");*/
+	    return trace_contour(pe_temp, z_level, num_active, *contr_isclosed);
+	}
+    }
+    if (*contr_isclosed) {
+	/* Look for something to start with inside: */
+      /*printf("<gen_one_contour> contr_isclosed TRUE \n");*/
+	pe_temp = p_edges;
+	while (pe_temp) {
+	    if (pe_temp->is_active && (pe_temp->position != BOUNDARY))
+		break;
+	    pe_temp = pe_temp->next;
+	}
+	if (!pe_temp) {
+	    *num_active = 0;
+	    fprintf(stderr, "gen_one_contour: no contour found\n");
+	    return NULL;
+	} else {
+	    *contr_isclosed = TRUE;
+	    return trace_contour(pe_temp, z_level, num_active, *contr_isclosed);
+	}
+    }
+    printf("<gen_one_contour> We should never be here, but lint... \n"); 
+    return NULL;		/* We should never be here, but lint... */
+}
+
+/*
+ * Search the data base along a contour starts at the edge pe_start until
+ * a boundary edge is detected or until we close the loop back to pe_start.
+ * Returns a linked list of all the points on the contour
+ * Also decreases num_active by the number of points on contour.
+ */
+static struct cntr_struct *
+trace_contour(pe_start, z_level, num_active, contr_isclosed)
+    struct edge_struct *pe_start; /* edge to start contour input */
+    double z_level;		/* Z level of contour input */
+    int *num_active;		/* number of active edges in/out */
+    TBOOLEAN contr_isclosed;	/* open or closed contour line (input) */
+{
+    struct cntr_struct *p_cntr, *pc_tail;
+    struct edge_struct *p_edge, *p_next_edge;
+    struct poly_struct *p_poly, *PLastpoly = NULL;
+    int i;
+
+    p_edge = pe_start;		/* first edge to start contour */
+
+    /* Generate the header of the contour - the point on pe_start. */
+    if (! contr_isclosed) {
+	pe_start->is_active = FALSE;
+	(*num_active)--;
+    }
+    if (p_edge->poly[0] || p_edge->poly[1]) {	/* more than one point */
+
+	p_cntr = pc_tail = update_cntr_pt(pe_start, z_level);	/* first point */
+
+	do {
+	    /* Find polygon to continue (Not where we came from - PLastpoly): */
+	    if (p_edge->poly[0] == PLastpoly)
+		p_poly = p_edge->poly[1];
+	    else
+		p_poly = p_edge->poly[0];
+	    p_next_edge = NULL;	/* In case of error, remains NULL. */
+	    for (i = 0; i < 3; i++)	/* Test the 3 edges of the polygon: */
+		if (p_poly->edge[i] != p_edge)
+		    if (p_poly->edge[i]->is_active)
+			p_next_edge = p_poly->edge[i];
+	    if (!p_next_edge) {	/* Error exit */
+		pc_tail->next = NULL;
+		free_contour(p_cntr);
+		fprintf(stderr, "trace_contour: unexpected end of contour\n");
+		return NULL;
+	    }
+	    p_edge = p_next_edge;
+	    PLastpoly = p_poly;
+	    p_edge->is_active = FALSE;
+	    (*num_active)--;
+
+	    /* Do not allocate contour points on diagonal edges */
+	    if (p_edge->position != DIAGONAL) {
+
+		pc_tail->next = update_cntr_pt(p_edge, z_level);
+
+		/* Remove nearby points */
+		if (fuzzy_equal(pc_tail, pc_tail->next)) {
+
+		    free((char *) pc_tail->next);
+		} else
+		    pc_tail = pc_tail->next;
+	    }
+	} while ((p_edge != pe_start) && (p_edge->position != BOUNDARY));
+
+	pc_tail->next = NULL;
+
+	/* For closed contour the first and last point should be equal */
+	if (pe_start == p_edge) {
+	    (p_cntr->X) = (pc_tail->X);
+	    (p_cntr->Y) = (pc_tail->Y);
+	}
+    } else {			/* only one point, forget it */
+	p_cntr = NULL;
+    }
+
+    return p_cntr;
+}
+
+/*
+ * Allocates one contour location and update it to to correct position
+ * according to z_level and edge p_edge.
+ */
+static struct cntr_struct *
+update_cntr_pt(p_edge, z_level)
+struct edge_struct *p_edge;
+double z_level;
+{
+    double t;
+    struct cntr_struct *p_cntr;
+
+    t = (z_level - p_edge->vertex[0]->z) /
+	(p_edge->vertex[1]->z - p_edge->vertex[0]->z);
+
+    /* test if t is out of interval [0:1] (should not happen but who knows ...) */
+    /*if(t>1) printf(" <update_cntr_pt> t >1 !\n");*/
+    /*if(t<0) printf(" <update_cntr_pt> t negatif !\n");*/
+    t = (t < 0.0 ? 0.0 : t);
+    t = (t > 1.0 ? 1.0 : t);
+    /*printf(" <update_cntr_pt> Point 0 %g %g %g \n",p_edge->vertex[0]->x,p_edge->vertex[0]->y,p_edge->vertex[0]->z);*/
+    /*printf(" <update_cntr_pt> Point 1 %g %g %g \n",p_edge->vertex[1]->x,p_edge->vertex[1]->y,p_edge->vertex[1]->z);*/
+    p_cntr = (struct cntr_struct *)
+	gp_alloc(sizeof(struct cntr_struct), "contour cntr_struct");
+
+    p_cntr->X = p_edge->vertex[1]->x * t +
+	p_edge->vertex[0]->x * (1 - t);
+    p_cntr->Y = p_edge->vertex[1]->y * t +
+	p_edge->vertex[0]->y * (1 - t);
+    /*printf(" <update_cntr_pt> p_cntr X %g Y %g \n",p_cntr->X,p_cntr->Y);*/
+    return p_cntr;
+}
+
+/* Simple routine to decide if two contour points are equal by
+ * calculating the relative error (< EPSILON).  */
+/* HBB 20010121: don't use absolute value 'zero' to compare to data
+ * values. */
+static int
+fuzzy_equal(p_cntr1, p_cntr2)
+    struct cntr_struct *p_cntr1, *p_cntr2;
+{
+    double unit_x, unit_y;
+    unit_x = fabs(x_max - x_min);		/* reference */
+    unit_y = fabs(y_max - y_min);
+    return ((fabs(p_cntr1->X - p_cntr2->X) < unit_x * EPSILON)
+	    && (fabs(p_cntr1->Y - p_cntr2->Y) < unit_y * EPSILON));
+}
+
+/*
+ * Generate the triangles.
+ * Returns the lists (edges & polys) via pointers to their heads.
+ */
+static void
+gen_triangle(num_isolines, iso_lines, p_polys, p_edges)
+int num_isolines;		/* number of iso-lines input */
+struct iso_curve *iso_lines;	/* iso-lines input */
+struct poly_struct **p_polys;	/* list of polygons output */
+struct edge_struct **p_edges;	/* list of edges output */
+{
+    int i, j, grid_x_max = iso_lines->p_count;
+    struct edge_struct *p_edge1, *p_edge2, *edge0, *edge1, *edge2, *pe_tail,
+	   *pe_tail2, *pe_temp;
+    struct poly_struct *pp_tail, *lower_tri, *upper_tri;
+    /* HBB 980308: need to tag *each* of them as ! */
+    struct coordinate  *p_vrtx1,  * p_vrtx2;
+
+    (*p_polys) = pp_tail = NULL;	/* clear lists */
+    (*p_edges) = pe_tail = NULL;
+
+    p_vrtx1 = iso_lines->points;	/* first row of vertices */
+    p_edge1 = pe_tail = NULL;	/* clear list of edges */
+
+    /* Generate edges of first row */
+    for (j = 0; j < grid_x_max - 1; j++)
+	add_edge(p_vrtx1 + j, p_vrtx1 + j + 1, &p_edge1, &pe_tail);
+
+    (*p_edges) = p_edge1;	/* update main list */
+
+
+    /*
+     * Combines vertices to edges and edges to triangles:
+     * ==================================================
+     * The edges are stored in the edge list, referenced by p_edges
+     * (pe_tail points on last edge).
+     *
+     * Temporary pointers:
+     * 1. p_edge2: Top horizontal edge list:      +-----------------------+ 2
+     * 2. p_tail : end of middle edge list:       |\  |\  |\  |\  |\  |\  |
+     *                                            |  \|  \|  \|  \|  \|  \|
+     * 3. p_edge1: Bottom horizontal edge list:   +-----------------------+ 1
+     *
+     * pe_tail2  : end of list beginning at p_edge2
+     * pe_temp   : position inside list beginning at p_edge1
+     * p_edges   : head of the master edge list (part of our output)
+     * p_vrtx1   : start of lower row of input vertices
+     * p_vrtx2   : start of higher row of input vertices
+     *
+     * The routine generates two triangle            Lower      Upper 1  
+     * upper one and lower one:                     | \           ----   
+     * (Nums. are edges order in polys)            0|   \1       0\   |2 
+     * The polygons are stored in the polygon        ----           \ |  
+     * list (*p_polys) (pp_tail points on             2                  
+     * last polygon).
+     *                                                        1
+     *                                                   -----------
+     * In addition, the edge lists are updated -        | \   0     |
+     * each edge has two pointers on the two            |   \       |
+     * (one active if boundary) polygons which         0|1   0\1   0|1
+     * uses it. These two pointer to polygons           |       \   |
+     * are named: poly[0], poly[1]. The diagram         |    1    \ |
+     * on the right show how they are used for the       -----------
+     * upper and lower polygons (INNER_MESH polygons only).  0
+     */
+
+    for (i = 1; i < num_isolines; i++) {
+	/* Read next column and gen. polys. */
+	iso_lines = iso_lines->next;
+
+	p_vrtx2 = iso_lines->points;	/* next row of vertices */
+	p_edge2 = pe_tail2 = NULL;	/* clear top horizontal list */
+	pe_temp = p_edge1;	/* pointer in bottom list */
+
+	/*
+	 * Generate edges and triagles for next row:
+	 */
+
+	/* generate first vertical edge */
+	edge2 = add_edge(p_vrtx1, p_vrtx2, p_edges, &pe_tail);
+
+	for (j = 0; j < grid_x_max - 1; j++) {
+
+	    /* copy vertical edge for lower triangle */
+	    edge0 = edge2;
+
+	    if (pe_temp && pe_temp->vertex[0] == p_vrtx1 + j) {
+		/* test lower edge */
+		edge2 = pe_temp;
+		pe_temp = pe_temp->next;
+	    } else {
+		edge2 = NULL;	/* edge is undefined */
+	    }
+
+	    /* generate diagonal edge */
+	    edge1 = add_edge(p_vrtx1 + j + 1, p_vrtx2 + j, p_edges, &pe_tail);
+	    if (edge1)
+		edge1->position = DIAGONAL;
+
+	    /* generate lower triangle */
+	    lower_tri = add_poly(edge0, edge1, edge2, p_polys, &pp_tail);
+
+	    /* copy diagonal edge for upper triangle */
+	    edge0 = edge1;
+
+	    /* generate upper edge */
+	    edge1 = add_edge(p_vrtx2 + j, p_vrtx2 + j + 1, &p_edge2, &pe_tail2);
+
+	    /* generate vertical edge */
+	    edge2 = add_edge(p_vrtx1 + j + 1, p_vrtx2 + j + 1, p_edges, &pe_tail);
+
+	    /* generate upper triangle */
+	    upper_tri = add_poly(edge0, edge1, edge2, p_polys, &pp_tail);
+	}
+
+	if (p_edge2) {
+	    /* HBB 19991130 bugfix: if p_edge2 list is empty,
+	     * don't change p_edges list! Crashes by access
+	     * to NULL pointer pe_tail, the second time through,
+	     * otherwise */
+	    if ((*p_edges)) {	/* Chain new edges to main list. */
+		pe_tail->next = p_edge2;
+		pe_tail = pe_tail2;
+	    } else {
+		(*p_edges) = p_edge2;
+		pe_tail = pe_tail2;
+	    }
+	}
+
+	/* this row finished, move list heads up one row: */
+	p_edge1 = p_edge2;
+	p_vrtx1 = p_vrtx2;
+    }
+
+    /* Update the boundary flag, saved in each edge, and update indexes: */
+
+    pe_temp = (*p_edges);
+
+    while (pe_temp) {
+	if ((!(pe_temp->poly[0])) || (!(pe_temp->poly[1])))
+	    (pe_temp->position) = BOUNDARY;
+	pe_temp = pe_temp->next;
+    }
+}
+
+/*
+ * Calculate minimum and maximum values
+ */
+static void
+calc_min_max(num_isolines, iso_lines, xx_min, yy_min, zz_min, xx_max, yy_max, zz_max)
+    int num_isolines;		/* number of iso-lines input */
+    struct iso_curve *iso_lines; /* iso-lines input */
+    double *xx_min, *yy_min, *zz_min, *xx_max, *yy_max, *zz_max; /* min/max values in/out */
+{
+    int i, j, grid_x_max;
+    struct coordinate  *vertex;
+    /*printf("<calc_min_max> :  iso_lines->p_count %ld \n",iso_lines->p_count);*/
+    grid_x_max = iso_lines->p_count;	/* number of vertices per iso_line */
+
+    (*xx_min) = (*yy_min) = (*zz_min) = VERYLARGE;	/* clear min/max values */
+    (*xx_max) = (*yy_max) = (*zz_max) = -VERYLARGE;
+    /*printf(" <calc_min_max> %d \n",num_isolines);*/
+    for (j = 0; j < num_isolines; j++) {
+      /*printf(" <calc_min_max> iso_lines %lx %lx %d min %g max %g\n",      */
+      /*       iso_lines,iso_lines->points,grid_x_max,(*zz_min),(*zz_max)); */
+	vertex = iso_lines->points;
+
+	for (i = 0; i < grid_x_max; i++) {
+	    if (vertex[i].type != UNDEFINED) {
+		if (vertex[i].x > (*xx_max))
+		    (*xx_max) = vertex[i].x;
+		if (vertex[i].y > (*yy_max))
+		    (*yy_max) = vertex[i].y;
+		if (vertex[i].z > (*zz_max))
+		    (*zz_max) = vertex[i].z;
+		if (vertex[i].x < (*xx_min))
+		    (*xx_min) = vertex[i].x;
+		if (vertex[i].y < (*yy_min))
+		    (*yy_min) = vertex[i].y;
+		if (vertex[i].z < (*zz_min))
+		    (*zz_min) = vertex[i].z;
+		
+	    }
+	}
+	iso_lines = iso_lines->next;
+	/*printf(" End of loop calc_min_max %d \n", j);*/
+    }
+    /* HBB 20000426: this code didn't take into account that axes might
+     * be logscaled... */
+#if 0
+    /* HBB 20001220: DON'T. The values are actually already stored
+     * logarithmized, as should be! */
+    axis_unlog_interval(FIRST_X_AXIS, xx_min, xx_max, 0);
+    axis_unlog_interval(FIRST_Y_AXIS, yy_min, yy_max, 0);
+    axis_unlog_interval(FIRST_Z_AXIS, zz_min, zz_max, 0);
+#endif
+
+    /* 
+     * fprintf(stderr," x: %g, %g\n", (*xx_min), (*xx_max));
+     * fprintf(stderr," y: %g, %g\n", (*yy_min), (*yy_max));
+     * fprintf(stderr," z: %g, %g\n", (*zz_min), (*zz_max));
+     */
+}
+
+/*
+ * Generate new edge and append it to list, but only if both vertices are 
+ * defined. The list is referenced by p_edge and pe_tail (p_edge points on 
+ * first edge and pe_tail on last one).
+ * Note, the list may be empty (pe_edge==pe_tail==NULL) on entry and exit.
+ */
+static struct edge_struct *
+add_edge(point0, point1, p_edge, pe_tail)
+struct coordinate  *point0;	/* 2 vertices input */
+struct coordinate  *point1;
+struct edge_struct **p_edge, **pe_tail;		/* pointers to edge list in/out */
+{
+    struct edge_struct *pe_temp = NULL;
+
+#if 1
+    if (point0->type == INRANGE && point1->type == INRANGE) {
+#else
+    if (point0->type != UNDEFINED && point1->type != UNDEFINED) {
+#endif
+
+	pe_temp = (struct edge_struct *)
+	    gp_alloc(sizeof(struct edge_struct), "contour edge");
+
+	pe_temp->poly[0] = NULL;	/* clear links           */
+	pe_temp->poly[1] = NULL;
+	pe_temp->vertex[0] = point0;	/* First vertex of edge. */
+	pe_temp->vertex[1] = point1;	/* Second vertex of edge. */
+	pe_temp->next = NULL;
+	pe_temp->position = INNER_MESH;		/* default position in mesh */
+
+	if ((*pe_tail)) {
+	    (*pe_tail)->next = pe_temp;		/* Stick new record as last one. */
+	} else {
+	    (*p_edge) = pe_temp;	/* start new list if empty */
+	}
+	(*pe_tail) = pe_temp;	/* continue to last record. */
+
+    }
+    return pe_temp;		/* returns NULL, if no edge allocated */
+}
+
+/*
+ * Generate new triangle and append it to list, but only if all edges are defined.
+ * The list is referenced by p_poly and pp_tail (p_poly points on first ploygon
+ * and pp_tail on last one).
+ * Note, the list may be empty (pe_ploy==pp_tail==NULL) on entry and exit.
+ */
+static struct poly_struct *
+add_poly(edge0, edge1, edge2, p_poly, pp_tail)
+    struct edge_struct *edge0, *edge1, *edge2;	/* 3 edges input */
+    struct poly_struct **p_poly, **pp_tail;		/* pointers to polygon list in/out */
+{
+    struct poly_struct *pp_temp = NULL;
+
+    if (edge0 && edge1 && edge2) {
+
+	pp_temp = (struct poly_struct *)
+	    gp_alloc(sizeof(struct poly_struct), "contour polygon");
+
+	pp_temp->edge[0] = edge0;	/* First edge of triangle */
+	pp_temp->edge[1] = edge1;	/* Second one             */
+	pp_temp->edge[2] = edge2;	/* Third one              */
+	pp_temp->next = NULL;
+
+	if (edge0->poly[0])	/* update edge0 */
+	    edge0->poly[1] = pp_temp;
+	else
+	    edge0->poly[0] = pp_temp;
+
+	if (edge1->poly[0])	/* update edge1 */
+	    edge1->poly[1] = pp_temp;
+	else
+	    edge1->poly[0] = pp_temp;
+
+	if (edge2->poly[0])	/* update edge2 */
+	    edge2->poly[1] = pp_temp;
+	else
+	    edge2->poly[0] = pp_temp;
+
+	if ((*pp_tail))		/* Stick new record as last one. */
+	    (*pp_tail)->next = pp_temp;
+	else
+	    (*p_poly) = pp_temp;	/* start new list if empty */
+
+	(*pp_tail) = pp_temp;	/* continue to last record. */
+
+    }
+    return pp_temp;		/* returns NULL, if no edge allocated */
+}
+
+
+
+/*
+ * Calls the (hopefully) desired interpolation/approximation routine.
+ */
+static void
+put_contour(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_isclosed)
+    struct cntr_struct *p_cntr;	/* contour structure input */
+    double z_level;		/* Z level of contour input */
+    double xx_min, xx_max, yy_min, yy_max; /* minimum/maximum values input */
+    TBOOLEAN contr_isclosed;		/* contour line closed? (input) */
+{
+
+    if (!p_cntr)
+	return;			/* Nothing to do if it is empty contour. */
+
+    switch (interp_kind) {
+    case CONTOUR_KIND_LINEAR:	/* No interpolation/approximation. */
+	put_contour_nothing(p_cntr);
+	break;
+    case CONTOUR_KIND_CUBIC_SPL: /* Cubic spline interpolation. */
+	put_contour_cubic(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max,
+			  chk_contour_kind(p_cntr, contr_isclosed));
+
+	break;
+    case CONTOUR_KIND_BSPLINE:	/* Bspline approximation. */
+	put_contour_bspline(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max,
+			    chk_contour_kind(p_cntr, contr_isclosed));
+	break;
+    }
+    free_contour(p_cntr);
+}
+
+/*
+ * Simply puts contour coordinates in order with no interpolation or
+ * approximation.
+ */
+static void
+put_contour_nothing(p_cntr)
+struct cntr_struct *p_cntr;
+{
+    while (p_cntr) {
+	add_cntr_point(p_cntr->X, p_cntr->Y);
+	p_cntr = p_cntr->next;
+    }
+    end_crnt_cntr();
+}
+
+/*
+ * for some reason contours are never flagged as 'isclosed'
+ * if first point == last point, set flag accordingly
+ *
+ */
+
+static int
+chk_contour_kind(p_cntr, contr_isclosed)
+    struct cntr_struct *p_cntr;
+    TBOOLEAN contr_isclosed;
+{
+    struct cntr_struct *pc_tail = NULL;
+    TBOOLEAN current_contr_isclosed;
+
+    /*fprintf(stderr, "check_contour_kind: current contour_kind value is %d\n", contour_kind);*/
+
+    current_contr_isclosed = contr_isclosed;
+
+    if (! contr_isclosed) {
+	pc_tail = p_cntr;
+	while (pc_tail->next)
+	    pc_tail = pc_tail->next;	/* Find last point. */
+
+	/* test if first and last point are equal */
+	if (fuzzy_equal(pc_tail, p_cntr)) {
+	    current_contr_isclosed = TRUE;
+	    fprintf(stderr, "check_contour_kind: contr_isclosed changed to %d\n", current_contr_isclosed);
+	}
+    }
+    return (current_contr_isclosed);
+}
+
+/*
+ * Generate a cubic spline curve through the points (x_i,y_i) which are
+ * stored in the linked list p_cntr.
+ * The spline is defined as a 2d-function s(t) = (x(t),y(t)), where the
+ * parameter t is the length of the linear stroke.
+ */
+static void
+put_contour_cubic(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_isclosed)
+    struct cntr_struct *p_cntr;
+    double z_level, xx_min, xx_max, yy_min, yy_max;
+    TBOOLEAN contr_isclosed;
+{
+    int num_pts, num_intpol;
+    double unit_x, unit_y;	/* To define norm (x,y)-plane */
+    double *delta_t;		/* Interval length t_{i+1}-t_i */
+    double *d2x, *d2y;		/* Second derivatives x''(t_i), y''(t_i) */
+    struct cntr_struct *pc_tail;
+
+    num_pts = count_contour(p_cntr);	/* Number of points in contour. */
+
+    pc_tail = p_cntr;		/* Find last point. */
+    while (pc_tail->next)
+	pc_tail = pc_tail->next;
+
+    if (contr_isclosed) {
+	/* Test if first and last point are equal (should be) */
+	if (!fuzzy_equal(pc_tail, p_cntr)) {
+	    pc_tail->next = p_cntr;	/* Close contour list - make it circular. */
+	    num_pts++;
+	}
+    }
+    delta_t = (double *) gp_alloc(num_pts * sizeof(double), "contour delta_t");
+    d2x = (double *) gp_alloc(num_pts * sizeof(double), "contour d2x");
+    d2y = (double *) gp_alloc(num_pts * sizeof(double), "contour d2y");
+
+    /* Width and height of the grid is used as a unit length (2d-norm) */
+    unit_x = xx_max - x_min;
+    unit_y = yy_max - y_min;
+    /* FIXME HBB 20010121: 'zero' should not be used as an absolute
+     * figure to compare to data */
+    unit_x = (unit_x > zero ? unit_x : zero);	/* should not be zero */
+    unit_y = (unit_y > zero ? unit_y : zero);
+
+    if (num_pts > 2) {
+	/*
+	 * Calculate second derivatives d2x[], d2y[] and interval lengths delta_t[]:
+	 */
+	if (!gen_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t,
+			      contr_isclosed, unit_x, unit_y)) {
+	    free((char *) delta_t);
+	    free((char *) d2x);
+	    free((char *) d2y);
+	    if (contr_isclosed)
+		pc_tail->next = NULL;	/* Un-circular list */
+	    return;
+	}
+    }
+    /* If following (num_pts > 1) is TRUE then exactly 2 points in contour.  */
+    else if (num_pts > 1) {
+	/* set all second derivatives to zero, interval length to 1 */
+	d2x[0] = 0.;
+	d2y[0] = 0.;
+	d2x[1] = 0.;
+	d2y[1] = 0.;
+	delta_t[0] = 1.;
+    } else {			/* Only one point ( ?? ) - ignore it. */
+	free((char *) delta_t);
+	free((char *) d2x);
+	free((char *) d2y);
+	if (contr_isclosed)
+	    pc_tail->next = NULL;	/* Un-circular list */
+	return;
+    }
+
+    /* Calculate "num_intpol" interpolated values */
+    num_intpol = 1 + (num_pts - 1) * contour_pts;	/* global: contour_pts */
+    intp_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t, num_intpol);
+
+    free((char *) delta_t);
+    free((char *) d2x);
+    free((char *) d2y);
+
+    if (contr_isclosed)
+	pc_tail->next = NULL;	/* Un-circular list */
+
+    end_crnt_cntr();
+}
+
+
+/*
+ * Find Bspline approximation for this data set.
+ * Uses global variable contour_pts to determine number of samples per
+ * interval, where the knot vector intervals are assumed to be uniform, and
+ * global variable contour_order for the order of Bspline to use.
+ */
+static void
+put_contour_bspline(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_isclosed)
+struct cntr_struct *p_cntr;
+double z_level, xx_min, xx_max, yy_min, yy_max;
+TBOOLEAN contr_isclosed;
+{
+    int num_pts;
+    int order = contour_order - 1;
+
+    num_pts = count_contour(p_cntr);	/* Number of points in contour. */
+    if (num_pts < 2)
+	return;			/* Can't do nothing if empty or one points! */
+    /* Order must be less than number of points in curve - fix it if needed. */
+    if (order > num_pts - 1)
+	order = num_pts - 1;
+
+    gen_bspline_approx(p_cntr, num_pts, order, contr_isclosed);
+    end_crnt_cntr();
+}
+
+/*
+ * Free all elements in the contour list.
+ */
+static void
+free_contour(p_cntr)
+struct cntr_struct *p_cntr;
+{
+    struct cntr_struct *pc_temp;
+
+    while (p_cntr) {
+	pc_temp = p_cntr;
+	p_cntr = p_cntr->next;
+	free((char *) pc_temp);
+    }
+}
+
+/*
+ * Counts number of points in contour.
+ */
+static int
+count_contour(p_cntr)
+struct cntr_struct *p_cntr;
+{
+    int count = 0;
+
+    while (p_cntr) {
+	count++;
+	p_cntr = p_cntr->next;
+    }
+    return count;
+}
+
+/*
+ * Find second derivatives (x''(t_i),y''(t_i)) of cubic spline interpolation
+ * through list of points (x_i,y_i). The parameter t is calculated as the
+ * length of the linear stroke. The number of points must be at least 3.
+ * Note: For closed contours the first and last point must be equal.
+ */
+static int
+gen_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t, contr_isclosed, unit_x, unit_y)
+int num_pts;			/* Number of points (num_pts>=3), input */
+struct cntr_struct *p_cntr;	/* List of points (x(t_i),y(t_i)), input */
+double d2x[], d2y[],		/* Second derivatives (x''(t_i),y''(t_i)), output */
+ delta_t[];			/* List of interval lengths t_{i+1}-t_{i}, output */
+TBOOLEAN contr_isclosed;	/* Closed or open contour?, input  */
+double unit_x, unit_y;		/* Unit length in x and y (norm=1), input */
+{
+    int n, i;
+    double norm;
+    tri_diag *m;		/* The tri-diagonal matrix is saved here. */
+    struct cntr_struct *pc_temp;
+
+    m = (tri_diag *) gp_alloc(num_pts * sizeof(tri_diag), "contour tridiag m");
+
+    /*
+     * Calculate first differences in (d2x[i], d2y[i]) and interval lengths
+     * in delta_t[i]:
+     */
+    pc_temp = p_cntr;
+    for (i = 0; i < num_pts - 1; i++) {
+	d2x[i] = pc_temp->next->X - pc_temp->X;
+	d2y[i] = pc_temp->next->Y - pc_temp->Y;
+	/*
+	 * The norm of a linear stroke is calculated in "normal coordinates"
+	 * and used as interval length:
+	 */
+	delta_t[i] = sqrt(SQR(d2x[i] / unit_x) + SQR(d2y[i] / unit_y));
+
+	d2x[i] /= delta_t[i];	/* first difference, with unit norm: */
+	d2y[i] /= delta_t[i];	/*   || (d2x[i], d2y[i]) || = 1      */
+
+	pc_temp = pc_temp->next;
+    }
+
+    /*
+     * Setup linear system:  m * x = b
+     */
+    n = num_pts - 2;		/* Without first and last point */
+    if (contr_isclosed) {
+	/* First and last points must be equal for closed contours */
+	delta_t[num_pts - 1] = delta_t[0];
+	d2x[num_pts - 1] = d2x[0];
+	d2y[num_pts - 1] = d2y[0];
+	n++;			/* Add last point (= first point) */
+    }
+    for (i = 0; i < n; i++) {
+	/* Matrix M, mainly tridiagonal with cyclic second index ("j = j+n mod n") */
+	m[i][0] = delta_t[i];	/* Off-diagonal element M_{i,i-1} */
+	m[i][1] = 2. * (delta_t[i] + delta_t[i + 1]);	/* M_{i,i} */
+	m[i][2] = delta_t[i + 1];	/* Off-diagonal element M_{i,i+1} */
+
+	/* Right side b_x and b_y */
+	d2x[i] = (d2x[i + 1] - d2x[i]) * 6.;
+	d2y[i] = (d2y[i + 1] - d2y[i]) * 6.;
+
+	/*
+	 * If the linear stroke shows a cusps of more than 90 degree, the right
+	 * side is reduced to avoid oscillations in the spline:
+	 */
+	norm = sqrt(SQR(d2x[i] / unit_x) + SQR(d2y[i] / unit_y)) / 8.5;
+
+	if (norm > 1.) {
+	    d2x[i] /= norm;
+	    d2y[i] /= norm;
+	    /* The first derivative will not be continuous */
+	}
+    }
+
+    if (!contr_isclosed) {
+	/* Third derivative is set to zero at both ends */
+	m[0][1] += m[0][0];	/* M_{0,0}     */
+	m[0][0] = 0.;		/* M_{0,n-1}   */
+	m[n - 1][1] += m[n - 1][2];	/* M_{n-1,n-1} */
+	m[n - 1][2] = 0.;	/* M_{n-1,0}   */
+    }
+    /* Solve linear systems for d2x[] and d2y[] */
+
+
+    if (solve_cubic_1(m, n)) {	/* Calculate Cholesky decomposition */
+	solve_cubic_2(m, d2x, n);	/* solve M * d2x = b_x */
+	solve_cubic_2(m, d2y, n);	/* solve M * d2y = b_y */
+
+    } else {			/* Should not happen, but who knows ... */
+	free((char *) m);
+	return FALSE;
+    }
+
+    /* Shift all second derivatives one place right and abdate end points */
+    for (i = n; i > 0; i--) {
+	d2x[i] = d2x[i - 1];
+	d2y[i] = d2y[i - 1];
+    }
+    if (contr_isclosed) {
+	d2x[0] = d2x[n];
+	d2y[0] = d2y[n];
+    } else {
+	d2x[0] = d2x[1];	/* Third derivative is zero in */
+	d2y[0] = d2y[1];	/*     first and last interval */
+	d2x[n + 1] = d2x[n];
+	d2y[n + 1] = d2y[n];
+    }
+
+    free((char *) m);
+    return TRUE;
+}
+
+/*
+ * Calculate interpolated values of the spline function (defined via p_cntr
+ * and the second derivatives d2x[] and d2y[]). The number of tabulated
+ * values is n. On an equidistant grid n_intpol values are calculated.
+ */
+static void
+intp_cubic_spline(n, p_cntr, d2x, d2y, delta_t, n_intpol)
+int n;
+struct cntr_struct *p_cntr;
+double d2x[], d2y[], delta_t[];
+int n_intpol;
+{
+    double t, t_skip, t_max;
+    double x0, x1, x, y0, y1, y;
+    double d, hx, dx0, dx01, hy, dy0, dy01;
+    int i;
+
+    /* The length of the total interval */
+    t_max = 0.;
+    for (i = 0; i < n - 1; i++)
+	t_max += delta_t[i];
+
+    /* The distance between interpolated points */
+    t_skip = (1. - 1e-7) * t_max / (n_intpol - 1);
+
+    t = 0.;			/* Parameter value */
+    x1 = p_cntr->X;
+    y1 = p_cntr->Y;
+    add_cntr_point(x1, y1);	/* First point. */
+    t += t_skip;
+
+    for (i = 0; i < n - 1; i++) {
+	p_cntr = p_cntr->next;
+
+	d = delta_t[i];		/* Interval length */
+	x0 = x1;
+	y0 = y1;
+	x1 = p_cntr->X;
+	y1 = p_cntr->Y;
+	hx = (x1 - x0) / d;
+	hy = (y1 - y0) / d;
+	dx0 = (d2x[i + 1] + 2 * d2x[i]) / 6.;
+	dy0 = (d2y[i + 1] + 2 * d2y[i]) / 6.;
+	dx01 = (d2x[i + 1] - d2x[i]) / (6. * d);
+	dy01 = (d2y[i + 1] - d2y[i]) / (6. * d);
+	while (t <= delta_t[i]) {	/* t in current interval ? */
+	    x = x0 + t * (hx + (t - d) * (dx0 + t * dx01));
+	    y = y0 + t * (hy + (t - d) * (dy0 + t * dy01));
+	    add_cntr_point(x, y);	/* next point. */
+	    t += t_skip;
+	}
+	t -= delta_t[i];	/* Parameter t relative to start of next interval */
+    }
+}
+
+/*
+ * The following two procedures solve the special linear system which arise
+ * in cubic spline interpolation. If x is assumed cyclic ( x[i]=x[n+i] ) the
+ * equations can be written as (i=0,1,...,n-1):
+ *     m[i][0] * x[i-1] + m[i][1] * x[i] + m[i][2] * x[i+1] = b[i] .
+ * In matrix notation one gets M * x = b, where the matrix M is tridiagonal
+ * with additional elements in the upper right and lower left position:
+ *   m[i][0] = M_{i,i-1}  for i=1,2,...,n-1    and    m[0][0] = M_{0,n-1} ,
+ *   m[i][1] = M_{i, i }  for i=0,1,...,n-1
+ *   m[i][2] = M_{i,i+1}  for i=0,1,...,n-2    and    m[n-1][2] = M_{n-1,0}.
+ * M should be symmetric (m[i+1][0]=m[i][2]) and positiv definite.
+ * The size of the system is given in n (n>=1).
+ *
+ * In the first procedure the Cholesky decomposition M = C^T * D * C
+ * (C is upper triangle with unit diagonal, D is diagonal) is calculated.
+ * Return TRUE if decomposition exist.
+ */
+static int
+solve_cubic_1(m, n)
+tri_diag m[];
+int n;
+{
+    int i;
+    double m_ij, m_n, m_nn, d;
+
+    if (n < 1)
+	return FALSE;		/* Dimension should be at least 1 */
+
+    d = m[0][1];		/* D_{0,0} = M_{0,0} */
+    if (d <= 0.)
+	return FALSE;		/* M (or D) should be positiv definite */
+    m_n = m[0][0];		/*  M_{0,n-1}  */
+    m_nn = m[n - 1][1];		/* M_{n-1,n-1} */
+    for (i = 0; i < n - 2; i++) {
+	m_ij = m[i][2];		/*  M_{i,1}  */
+	m[i][2] = m_ij / d;	/* C_{i,i+1} */
+	m[i][0] = m_n / d;	/* C_{i,n-1} */
+	m_nn -= m[i][0] * m_n;	/* to get C_{n-1,n-1} */
+	m_n = -m[i][2] * m_n;	/* to get C_{i+1,n-1} */
+	d = m[i + 1][1] - m[i][2] * m_ij;	/* D_{i+1,i+1} */
+	if (d <= 0.)
+	    return FALSE;	/* Elements of D should be positiv */
+	m[i + 1][1] = d;
+    }
+    if (n >= 2) {		/* Complete last column */
+	m_n += m[n - 2][2];	/* add M_{n-2,n-1} */
+	m[n - 2][0] = m_n / d;	/* C_{n-2,n-1} */
+	m[n - 1][1] = d = m_nn - m[n - 2][0] * m_n;	/* D_{n-1,n-1} */
+	if (d <= 0.)
+	    return FALSE;
+    }
+    return TRUE;
+}
+
+/*
+ * The second procedure solves the linear system, with the Choleky
+ * decomposition calculated above (in m[][]) and the right side b given
+ * in x[]. The solution x overwrites the right side in x[].
+ */
+static void
+solve_cubic_2(m, x, n)
+tri_diag m[];
+double x[];
+int n;
+{
+    int i;
+    double x_n;
+
+    /* Division by transpose of C : b = C^{-T} * b */
+    x_n = x[n - 1];
+    for (i = 0; i < n - 2; i++) {
+	x[i + 1] -= m[i][2] * x[i];	/* C_{i,i+1} * x_{i} */
+	x_n -= m[i][0] * x[i];	/* C_{i,n-1} * x_{i} */
+    }
+    if (n >= 2)
+	x[n - 1] = x_n - m[n - 2][0] * x[n - 2];	/* C_{n-2,n-1} * x_{n-1} */
+
+    /* Division by D: b = D^{-1} * b */
+    for (i = 0; i < n; i++)
+	x[i] /= m[i][1];
+
+    /* Division by C: b = C^{-1} * b */
+    x_n = x[n - 1];
+    if (n >= 2)
+	x[n - 2] -= m[n - 2][0] * x_n;	/* C_{n-2,n-1} * x_{n-1} */
+    for (i = n - 3; i >= 0; i--) {
+	/*      C_{i,i+1} * x_{i+1} + C_{i,n-1} * x_{n-1} */
+	x[i] -= m[i][2] * x[i + 1] + m[i][0] * x_n;
+    }
+    return;
+}
+
+/*
+ * Solve tri diagonal linear system equation. The tri diagonal matrix is
+ * defined via matrix M, right side is r, and solution X i.e. M * X = R.
+ * Size of system given in n. Return TRUE if solution exist.
+ */
+/* not used any more in "contour.c", but in "spline.c" (21. Dec. 1995) ! */
+
+int
+solve_tri_diag(m, r, x, n)
+tri_diag m[];
+double r[], x[];
+int n;
+{
+    int i;
+    double t;
+
+    for (i = 1; i < n; i++) {	/* Eliminate element m[i][i-1] (lower diagonal). */
+	if (m[i - 1][1] == 0)
+	    return FALSE;
+	t = m[i][0] / m[i - 1][1];	/* Find ratio between the two lines. */
+/*      m[i][0] = m[i][0] - m[i-1][1] * t; */
+/* m[i][0] is not used any more (and set to 0 in the above line) */
+	m[i][1] = m[i][1] - m[i - 1][2] * t;
+	r[i] = r[i] - r[i - 1] * t;
+    }
+    /* Now do back subtitution - update the solution vector X: */
+    if (m[n - 1][1] == 0)
+	return FALSE;
+    x[n - 1] = r[n - 1] / m[n - 1][1];	/* Find last element. */
+    for (i = n - 2; i >= 0; i--) {
+	if (m[i][1] == 0)
+	    return FALSE;
+	x[i] = (r[i] - x[i + 1] * m[i][2]) / m[i][1];
+    }
+    return TRUE;
+}
+
+/*
+ * Generate a Bspline curve defined by all the points given in linked list p:
+ * Algorithm: using deBoor algorithm
+ * Note: if Curvekind is open contour than Open end knot vector is assumed,
+ *       else (closed contour) Float end knot vector is assumed.
+ * It is assumed that num_of_points is at least 2, and order of Bspline is less
+ * than num_of_points!
+ */
+static void
+gen_bspline_approx(p_cntr, num_of_points, order, contr_isclosed)
+struct cntr_struct *p_cntr;
+int num_of_points, order;
+TBOOLEAN contr_isclosed;
+{
+    int knot_index = 0, pts_count = 1;
+    double dt, t, next_t, t_min, t_max, x, y;
+    struct cntr_struct *pc_temp = p_cntr, *pc_tail = NULL;
+
+    /* If the contour is Closed one we must update few things:
+     * 1. Make the list temporary circular, so we can close the contour.
+     * 2. Update num_of_points - increase it by "order-1" so contour will be
+     *    closed. This will evaluate order more sections to close it!
+     */
+    if (contr_isclosed) {
+	pc_tail = p_cntr;
+	while (pc_tail->next)
+	    pc_tail = pc_tail->next;	/* Find last point. */
+
+	/* test if first and last point are equal */
+	if (fuzzy_equal(pc_tail, p_cntr)) {
+	    /* Close contour list - make it circular. */
+	    pc_tail->next = p_cntr->next;
+	    num_of_points += order - 1;
+	} else {
+	    pc_tail->next = p_cntr;
+	    num_of_points += order;
+	}
+    }
+    /* Find first (t_min) and last (t_max) t value to eval: */
+    t = t_min = fetch_knot(contr_isclosed, num_of_points, order, order);
+    t_max = fetch_knot(contr_isclosed, num_of_points, order, num_of_points);
+    next_t = t_min + 1.0;
+    knot_index = order;
+    dt = 1.0 / contour_pts;	/* Number of points per one section. */
+
+
+    while (t < t_max) {
+	if (t > next_t) {
+	    pc_temp = pc_temp->next;	/* Next order ctrl. pt. to blend. */
+	    knot_index++;
+	    next_t += 1.0;
+	}
+	eval_bspline(t, pc_temp, num_of_points, order, knot_index,
+		     contr_isclosed, &x, &y);	/* Next pt. */
+	add_cntr_point(x, y);
+	pts_count++;
+	/* As we might have some real number round off problems we do      */
+	/* the last point outside the loop                                 */
+	if (pts_count == contour_pts * (num_of_points - order) + 1)
+	    break;
+	t += dt;
+    }
+
+    /* Now do the last point */
+    eval_bspline(t_max - EPSILON, pc_temp, num_of_points, order, knot_index,
+		 contr_isclosed, &x, &y);
+    add_cntr_point(x, y);	/* Complete the contour. */
+
+    if (contr_isclosed)	/* Update list - un-circular it. */
+	pc_tail->next = NULL;
+}
+
+/*
+ * The routine to evaluate the B-spline value at point t using knot vector
+ * from function fetch_knot(), and the control points p_cntr.
+ * Returns (x, y) of approximated B-spline. Note that p_cntr points on the
+ * first control point to blend with. The B-spline is of order order.
+ */
+static void
+eval_bspline(t, p_cntr, num_of_points, order, j, contr_isclosed, x, y)
+double t;
+struct cntr_struct *p_cntr;
+int num_of_points, order, j;
+TBOOLEAN contr_isclosed;
+double *x, *y;
+{
+    int i, p;
+    double ti, tikp, *dx, *dy;	/* Copy p_cntr into it to make it faster. */
+
+    dx = (double *) gp_alloc((order + j) * sizeof(double), "contour b_spline");
+    dy = (double *) gp_alloc((order + j) * sizeof(double), "contour b_spline");
+
+    /* Set the dx/dy - [0] iteration step, control points (p==0 iterat.): */
+    for (i = j - order; i <= j; i++) {
+	dx[i] = p_cntr->X;
+	dy[i] = p_cntr->Y;
+	p_cntr = p_cntr->next;
+    }
+
+    for (p = 1; p <= order; p++) {	/* Iteration (b-spline level) counter. */
+	for (i = j; i >= j - order + p; i--) {	/* Control points indexing. */
+	    ti = fetch_knot(contr_isclosed, num_of_points, order, i);
+	    tikp = fetch_knot(contr_isclosed, num_of_points, order, i + order + 1 - p);
+	    if (ti == tikp) {	/* Should not be a problems but how knows... */
+	    } else {
+		dx[i] = dx[i] * (t - ti) / (tikp - ti) +	/* Calculate x. */
+		    dx[i - 1] * (tikp - t) / (tikp - ti);
+		dy[i] = dy[i] * (t - ti) / (tikp - ti) +	/* Calculate y. */
+		    dy[i - 1] * (tikp - t) / (tikp - ti);
+	    }
+	}
+    }
+    *x = dx[j];
+    *y = dy[j];
+    free((char *) dx);
+    free((char *) dy);
+}
+
+/*
+ * Routine to get the i knot from uniform knot vector. The knot vector
+ * might be float (Knot(i) = i) or open (where the first and last "order"
+ * knots are equal). contr_isclosed determines knot kind - open contour means
+ * open knot vector, and closed contour selects float knot vector.
+ * Note the knot vector is not exist and this routine simulates it existance
+ * Also note the indexes for the knot vector starts from 0.
+ */
+static double
+fetch_knot(contr_isclosed, num_of_points, order, i)
+    TBOOLEAN contr_isclosed;
+    int num_of_points, order, i;
+{
+    if(! contr_isclosed) {
+	if (i <= order)
+	    return 0.0;
+	else if (i <= num_of_points)
+	    return (double) (i - order);
+	else
+	    return (double) (num_of_points - order);
+    } else {
+	return (double) i;
+    }
+}
Index: /tags/PIGC_V1/SophyaPI/PIGcont/contour.h
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/contour.h	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/contour.h	(revision 1830)
@@ -0,0 +1,100 @@
+/*
+ * $Id: contour.h,v 1.1.1.1 2001-12-19 18:01:18 perderos Exp $
+ */
+
+/* GNUPLOT - contour.h */
+
+/*[
+ * Copyright 1999   Thomas Williams, Colin Kelley
+ *
+ * Permission to use, copy, and distribute this software and its
+ * documentation for any purpose with or 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.
+ *
+ * Permission to modify the software is granted, but not the right to
+ * distribute the complete modified source code.  Modifications are to
+ * be distributed as patches to the released version.  Permission to
+ * distribute binaries produced by compiling modified sources is granted,
+ * provided you
+ *   1. distribute the corresponding source modifications from the
+ *    released version in the form of a patch file along with the binaries,
+ *   2. add special version identification to distinguish your version
+ *    in addition to the base release version number,
+ *   3. provide your name and address as the primary contact for the
+ *    support of your modified version, and
+ *   4. retain our contact information in regard to use of the base
+ *    software.
+ * Permission to distribute the released version of the source code along
+ * with corresponding source modifications in the form of a patch file is
+ * granted with same provisions 2 through 4 for binary distributions.
+ *
+ * This software is provided "as is" without express or implied warranty
+ * to the extent permitted by applicable law.
+]*/
+
+#ifndef GNUPLOT_CONTOUR_H
+# define GNUPLOT_CONTOUR_H
+
+/* #if... / #include / #define collection: */
+
+/***
+#include "syscfg.h"
+***/
+#include "gp_types.h"
+
+#include "dynarray.h"
+#include "graph3d.h"
+#include "mymisc.h"
+#include <math.h>
+
+#define DEFAULT_CONTOUR_LEVELS 5
+#define DEFAULT_NUM_APPROX_PTS 5
+#define DEFAULT_CONTOUR_ORDER  4
+#define MAX_BSPLINE_ORDER      10
+
+/* Type definitions */
+
+typedef enum en_contour_kind {
+    /* Method of drawing the contour lines found */
+    CONTOUR_KIND_LINEAR,
+    CONTOUR_KIND_CUBIC_SPL,
+    CONTOUR_KIND_BSPLINE
+} t_contour_kind;
+
+typedef enum en_contour_levels_kind {
+    /* How contour levels are set */
+    LEVELS_AUTO,		/* automatically selected */
+
+    LEVELS_INCREMENTAL,		/* user specified start & incremnet */
+    LEVELS_DISCRETE,		/* user specified discrete levels */
+    LEVELS_NUM                  /* User sets contour number - AJOUT OP */
+     
+} t_contour_levels_kind;
+
+
+/* Used to allocate the tri-diag matrix. */
+typedef double tri_diag[3];
+
+/* Variables of contour.c needed by other modules: */
+
+extern char contour_format[];
+extern t_contour_kind contour_kind;
+extern t_contour_levels_kind contour_levels_kind;
+extern int contour_levels;
+extern int contour_order;
+extern int contour_pts;
+
+/* storage for z levels to draw contours at */
+extern dynarray dyn_contour_levels_list; 
+/*#define contour_levels_list ((double *)dyn_contour_levels_list.v) */
+extern double * contour_levels_list;
+
+/* Prototypes of functions exported by contour.c */
+
+struct gnuplot_contours *contour(int num_isolines, struct iso_curve *iso_lines);
+int solve_tri_diag(tri_diag m[], double r[], double x[], int n);
+int Get_Num_Of_Z_Levels();
+
+#endif /* GNUPLOT_CONTOUR_H */
Index: /tags/PIGC_V1/SophyaPI/PIGcont/dynarray.h
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/dynarray.h	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/dynarray.h	(revision 1830)
@@ -0,0 +1,63 @@
+/*
+ * $Id: dynarray.h,v 1.1.1.1 2001-12-19 18:01:18 perderos Exp $
+ */
+
+/* GNUPLOT - dynarray.h */
+
+/*[
+ * Copyright 1999   Thomas Williams, Colin Kelley
+ *
+ * Permission to use, copy, and distribute this software and its
+ * documentation for any purpose with or 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.
+ *
+ * Permission to modify the software is granted, but not the right to
+ * distribute the complete modified source code.  Modifications are to
+ * be distributed as patches to the released version.  Permission to
+ * distribute binaries produced by compiling modified sources is granted,
+ * provided you
+ *   1. distribute the corresponding source modifications from the
+ *    released version in the form of a patch file along with the binaries,
+ *   2. add special version identification to distinguish your version
+ *    in addition to the base release version number,
+ *   3. provide your name and address as the primary contact for the
+ *    support of your modified version, and
+ *   4. retain our contact information in regard to use of the base
+ *    software.
+ * Permission to distribute the released version of the source code along
+ * with corresponding source modifications in the form of a patch file is
+ * granted with same provisions 2 through 4 for binary distributions.
+ *
+ * This software is provided "as is" without express or implied warranty
+ * to the extent permitted by applicable law.
+]*/
+
+#ifndef DYNARRAY__H
+# define DYNARRAY__H
+
+/***
+# include "syscfg.h"
+# include "stdfn.h"
+***/
+#include <stdio.h>
+#include "mymisc.h"
+
+typedef struct dynarray {
+    long size;			/* alloced size of the array */
+    long end;			/* index of first unused entry */
+    long increment;		/* amount to increment size by, on realloc */
+    size_t entry_size;		/* size of the entries in this array */
+    void *v;		/* the vector itself */
+} dynarray;
+
+/* Prototypes */
+void init_dynarray (dynarray * array, size_t element, long size, long increment);
+void free_dynarray (dynarray * array);
+void extend_dynarray(dynarray * array, long increment);
+void resize_dynarray (dynarray * array, long newsize);
+void *nextfrom_dynarray (dynarray * array);
+void droplast_dynarray (dynarray * array);
+
+#endif /* DYNARRAY_H */
Index: /tags/PIGC_V1/SophyaPI/PIGcont/gp_types.h
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/gp_types.h	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/gp_types.h	(revision 1830)
@@ -0,0 +1,191 @@
+/*
+ * $Id: gp_types.h,v 1.1.1.1 2001-12-19 18:01:18 perderos Exp $
+ */
+
+/* GNUPLOT - gp_types.h */
+
+/*[
+ * Copyright 2000   Thomas Williams, Colin Kelley
+ *
+ * Permission to use, copy, and distribute this software and its
+ * documentation for any purpose with or 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.
+ *
+ * Permission to modify the software is granted, but not the right to
+ * distribute the complete modified source code.  Modifications are to
+ * be distributed as patches to the released version.  Permission to
+ * distribute binaries produced by compiling modified sources is granted,
+ * provided you
+ *   1. distribute the corresponding source modifications from the
+ *    released version in the form of a patch file along with the binaries,
+ *   2. add special version identification to distinguish your version
+ *    in addition to the base release version number,
+ *   3. provide your name and address as the primary contact for the
+ *    support of your modified version, and
+ *   4. retain our contact information in regard to use of the base
+ *    software.
+ * Permission to distribute the released version of the source code along
+ * with corresponding source modifications in the form of a patch file is
+ * granted with same provisions 2 through 4 for binary distributions.
+ *
+ * This software is provided "as is" without express or implied warranty
+ * to the extent permitted by applicable law.
+]*/
+
+#ifndef GNUPLOT_GPTYPES_H
+#define GNUPLOT_GPTYPES_H
+
+/***
+#include "syscfg.h"
+***/
+#include "mymisc.h"
+
+#define MAX_ID_LEN 50		/* max length of an identifier */
+#define MAX_LINE_LEN 1024	/* maximum number of chars allowed on line */
+
+#define DEG2RAD (M_PI / 180.0)
+
+/* These used to be #defines in plot.h. For easier debugging, I've
+ * converted most of them into enum's. */
+enum DATA_TYPES {
+	INTGR, CMPLX
+};
+
+enum PLOT_TYPE {
+	FUNC, DATA, FUNC3D, DATA3D
+};
+
+/* we explicitly assign values to the types, such that we can
+ * perform bit tests to see if the style involves points and/or lines
+ * bit 0 (val 1) = line, bit 1 (val 2) = point, bit 2 (val 4)= error
+ * This allows rapid decisions about the sample drawn into the key,
+ * for example.
+ */
+/* HBB 20010610: new enum, to make mnemonic names for these flags
+ * accessible everywhere */
+typedef enum e_PLOT_STYLE_FLAGS {
+    PLOT_STYLE_HAS_LINE      = (1<<0),
+    PLOT_STYLE_HAS_POINT     = (1<<1),
+    PLOT_STYLE_HAS_ERRORBAR  = (1<<2)
+} PLOT_STYLE_FLAGS;
+
+typedef enum PLOT_STYLE {
+    LINES        =  0*(1<<3) + PLOT_STYLE_HAS_LINE,
+    POINTSTYLE   =  1*(1<<3) + PLOT_STYLE_HAS_POINT,
+    IMPULSES     =  2*(1<<3) + PLOT_STYLE_HAS_LINE,
+    LINESPOINTS  =  3*(1<<3) + (PLOT_STYLE_HAS_POINT | PLOT_STYLE_HAS_LINE),
+    DOTS         =  4*(1<<3) + 0,
+    XERRORBARS   =  5*(1<<3) + (PLOT_STYLE_HAS_POINT | PLOT_STYLE_HAS_ERRORBAR),
+    YERRORBARS   =  6*(1<<3) + (PLOT_STYLE_HAS_POINT | PLOT_STYLE_HAS_ERRORBAR),
+    XYERRORBARS  =  7*(1<<3) + (PLOT_STYLE_HAS_POINT | PLOT_STYLE_HAS_ERRORBAR),
+    BOXXYERROR   =  8*(1<<3) + PLOT_STYLE_HAS_LINE,
+    BOXES        =  9*(1<<3) + PLOT_STYLE_HAS_LINE,
+    BOXERROR     = 10*(1<<3) + PLOT_STYLE_HAS_LINE,
+    STEPS        = 11*(1<<3) + PLOT_STYLE_HAS_LINE,
+    FSTEPS       = 12*(1<<3) + PLOT_STYLE_HAS_LINE,
+    HISTEPS      = 13*(1<<3) + PLOT_STYLE_HAS_LINE,
+    VECTOR       = 14*(1<<3) + PLOT_STYLE_HAS_LINE,
+    CANDLESTICKS = 15*(1<<3) + PLOT_STYLE_HAS_ERRORBAR,
+    /* FIXME HBB 20010214: shouldn't fincancebars have 4 or 5?
+     * They behave very much like errorbars, with the sole
+     * exception of the key entry ... */
+    FINANCEBARS  = 16*(1<<3) + PLOT_STYLE_HAS_LINE,
+    XERRORLINES  = 17*(1<<3) + (PLOT_STYLE_HAS_LINE | PLOT_STYLE_HAS_POINT | PLOT_STYLE_HAS_ERRORBAR),
+    YERRORLINES  = 18*(1<<3) + (PLOT_STYLE_HAS_LINE | PLOT_STYLE_HAS_POINT | PLOT_STYLE_HAS_ERRORBAR),
+    XYERRORLINES = 19*(1<<3) + (PLOT_STYLE_HAS_LINE | PLOT_STYLE_HAS_POINT | PLOT_STYLE_HAS_ERRORBAR)
+} PLOT_STYLE;
+
+typedef enum PLOT_SMOOTH { 
+    SMOOTH_NONE,
+    SMOOTH_ACSPLINES,
+    SMOOTH_BEZIER,
+    SMOOTH_CSPLINES,
+    SMOOTH_SBEZIER,
+    SMOOTH_UNIQUE
+} PLOT_SMOOTH;
+
+/* FIXME HBB 20000521: 'struct value' and its part, 'cmplx', should go
+ * into one of scanner/internal/standard/util .h, but I've yet to
+ * decide which of them */
+
+#if !(defined(ATARI)&&defined(__GNUC__)&&defined(_MATH_H)) &&  !(defined(MTOS)&&defined(__GNUC__)&&defined(_MATH_H)) /* FF's math.h has the type already */
+struct cmplx {
+	double real, imag;
+};
+#endif
+
+typedef struct value {
+    enum DATA_TYPES type;
+    union {
+	int int_val;
+	struct cmplx cmplx_val;
+    } v;
+} t_value;
+
+#if 0 /* HBB 20000521: move to command.h */
+typedef struct lexical_unit {	/* produced by scanner */
+	TBOOLEAN is_token;	/* true if token, false if a value */ 
+	struct value l_val;
+	int start_index;	/* index of first char in token */
+	int length;			/* length of token in chars */
+} lexical_unit;
+#endif /* 0 */
+
+#if 0 /* HBB 20000521: move to eval.h */
+typedef struct udft_entry {	/* user-defined function table entry */
+  struct udft_entry *next_udf;	/* pointer to next udf in linked list */
+  char *udf_name;		/* name of this function entry */
+  struct at_type *at;		/* pointer to action table to execute */
+  char *definition;		/* definition of function as typed */
+  struct value dummy_values[MAX_NUM_VAR]; /* current value of dummy variables */
+} udft_entry;
+
+
+typedef struct udvt_entry {	/* user-defined value table entry */
+	struct udvt_entry *next_udv; /* pointer to next value in linked list */
+	char *udv_name;			/* name of this value entry */
+	TBOOLEAN udv_undef;		/* true if not defined yet */
+	struct value udv_value;	/* value it has */
+} udvt_entry;
+
+typedef union argument {	/* p-code argument */
+	int j_arg;		/* offset for jump */
+	struct value v_arg;	/* constant value */
+	struct udvt_entry *udv_arg; /* pointer to dummy variable */
+	struct udft_entry *udf_arg; /* pointer to udf to execute */
+} argument;
+
+
+/* This type definition has to come after union argument has been declared. */
+#ifdef __ZTC__
+typedef int (*FUNC_PTR)(...);
+#else
+typedef int (*FUNC_PTR) __PROTO((union argument *arg));
+#endif
+
+typedef struct ft_entry {	/* standard/internal function table entry */
+	const char *f_name;	/* pointer to name of this function */
+	FUNC_PTR func;		/* address of function to call */
+} ft_entry;
+#endif /* 0 */
+
+
+/* Defines the type of a coordinate */
+/* INRANGE and OUTRANGE points have an x,y point associated with them */
+typedef enum coord_type {
+    INRANGE,			/* inside plot boundary */
+    OUTRANGE,			/* outside plot boundary, but defined */
+    UNDEFINED			/* not defined at all */
+} coord_type;
+
+
+typedef struct coordinate {
+    enum coord_type type;	/* see above */
+    double x, y, z;
+    double ylow, yhigh;	/* ignored in 3d */
+    double xlow, xhigh;	/* also ignored in 3d */
+} coordinate;
+
+#endif /* GNUPLOT_GPTYPES_H */
Index: /tags/PIGC_V1/SophyaPI/PIGcont/graph3d.h
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/graph3d.h	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/graph3d.h	(revision 1830)
@@ -0,0 +1,151 @@
+/*
+ * $Id: graph3d.h,v 1.1.1.1 2001-12-19 18:01:18 perderos Exp $
+ */
+
+/* GNUPLOT - graph3d.h */
+
+/*[
+ * Copyright 1999   Thomas Williams, Colin Kelley
+ *
+ * Permission to use, copy, and distribute this software and its
+ * documentation for any purpose with or 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.
+ *
+ * Permission to modify the software is granted, but not the right to
+ * distribute the complete modified source code.  Modifications are to
+ * be distributed as patches to the released version.  Permission to
+ * distribute binaries produced by compiling modified sources is granted,
+ * provided you
+ *   1. distribute the corresponding source modifications from the
+ *    released version in the form of a patch file along with the binaries,
+ *   2. add special version identification to distinguish your version
+ *    in addition to the base release version number,
+ *   3. provide your name and address as the primary contact for the
+ *    support of your modified version, and
+ *   4. retain our contact information in regard to use of the base
+ *    software.
+ * Permission to distribute the released version of the source code along
+ * with corresponding source modifications in the form of a patch file is
+ * granted with same provisions 2 through 4 for binary distributions.
+ *
+ * This software is provided "as is" without express or implied warranty
+ * to the extent permitted by applicable law.
+]*/
+
+#ifndef GNUPLOT_GRAPH3D_H
+# define GNUPLOT_GRAPH3D_H
+
+/* #if... / #include / #define collection: */
+
+/***
+#include "syscfg.h"
+
+#include "gadgets.h"
+#include "term_api.h"
+****/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include "gp_types.h"
+#include "mymisc.h"
+/* Function macros to map from user 3D space into normalized -1..1 */
+#define map_x3d(x) ((x-X_AXIS.min)*xscale3d-1.0)
+#define map_y3d(y) ((y-Y_AXIS.min)*yscale3d-1.0)
+#define map_z3d(z) ((z-floor_z)*zscale3d-1.0)
+
+/* Type definitions */
+
+typedef enum en_contour_placement {
+    /* Where to place contour maps if at all. */
+    CONTOUR_NONE,
+    CONTOUR_BASE,
+    CONTOUR_SRF,
+    CONTOUR_BOTH
+} t_contour_placement;
+
+typedef double transform_matrix[4][4]; /* HBB 990826: added */
+
+typedef struct gnuplot_contours {
+    struct gnuplot_contours *next;
+    struct coordinate *coords;
+    char isNewLevel;
+    char label[32];
+    int num_pts;
+#ifdef PM3D
+    double z;
+#endif
+} gnuplot_contours;
+
+typedef struct iso_curve {
+    struct iso_curve *next;
+    int p_max;			/* how many points are allocated */
+    int p_count;			/* count of points in points */
+    struct coordinate  *points;
+} iso_curve;
+
+typedef struct surface_points {
+    struct surface_points *next_sp; /* linked list */
+    int token;			/* last token nr, for second pass */
+    enum PLOT_TYPE plot_type;
+    enum PLOT_STYLE plot_style;
+    char *title;
+    int title_no_enhanced;	/* don't typeset title in enhanced mode */
+  /*struct lp_style_type lp_properties;*/
+    int has_grid_topology;
+    
+    /* Data files only - num of isolines read from file. For
+     * functions, num_iso_read is the number of 'primary' isolines (in
+     * x direction) */
+    int num_iso_read;		
+    struct gnuplot_contours *contours; /* NULL if not doing contours. */
+    struct iso_curve *iso_crvs;	/* the actual data */
+
+#ifdef PM3D_COLUMN
+    int pm3d_color_from_column;
+    coordval color_min;
+    coordval color_max;
+#endif
+    
+} surface_points;
+
+/* Variables of graph3d.c needed by other modules: */
+
+extern int xmiddle, ymiddle, xscaler, yscaler;
+extern double floor_z;
+extern double ceiling_z, base_z; /* made exportable for PM3D */
+extern transform_matrix trans_mat;
+extern double xscale3d, yscale3d, zscale3d;
+
+extern t_contour_placement draw_contour;
+extern TBOOLEAN	label_contours;
+
+extern TBOOLEAN	draw_surface;
+
+/* is hidden3d display wanted? */
+extern TBOOLEAN	hidden3d;
+
+extern float surface_rot_z;
+extern float surface_rot_x;
+extern float surface_scale;
+extern float surface_zscale;
+
+extern float ticslevel;
+
+#define ISO_SAMPLES 10		/* default number of isolines per splot */
+extern int iso_samples_1;
+extern int iso_samples_2;
+
+#ifdef USE_MOUSE
+extern int axis3d_o_x, axis3d_o_y, axis3d_x_dx, axis3d_x_dy, axis3d_y_dx, axis3d_y_dy;
+#endif
+
+/* Prototypes from file "graph3d.c" */
+
+void do_3dplot (struct surface_points *plots, int pcount, int quick);
+void map3d_position (struct position * pos, unsigned int *x,
+				  unsigned int *y, const char *what);
+
+
+#endif /* GNUPLOT_GRAPH3D_H */
Index: /tags/PIGC_V1/SophyaPI/PIGcont/myalloc.c
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/myalloc.c	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/myalloc.c	(revision 1830)
@@ -0,0 +1,104 @@
+#include "alloc.h"
+#include "graph3d.h"
+#include <math.h>
+#define TRACE_ALLOC(x)	
+struct frame_struct {
+    char *use;
+    int requested_size;
+    int pad;			/* preserve 8-byte alignment */
+    int checksum;
+};
+ 
+#define CHECKSUM_INT 0xcaac5e1f
+#define RESERVED_SIZE sizeof(struct frame_struct)
+#define CHECKSUM_CHAR 0xc5
+
+static long bytes_allocated = 0;
+
+static void
+mark(p, size, usage)
+struct frame_struct *p;
+unsigned long size;
+char *usage;
+{
+    p->use = usage;
+    p->requested_size = size;
+    p->checksum = (CHECKSUM_INT ^ (int) (p->use) ^ size);
+    ((unsigned char *) (p + 1))[size] = CHECKSUM_CHAR;
+}
+
+void *
+gp_alloc(size, usage)
+size_t size;
+const char *usage;
+{
+    struct frame_struct *p;
+    size_t total_size = size + RESERVED_SIZE + 1;
+
+    /*TRACE_ALLOC(("gp_alloc %d for %s\n", (int) size, usage ? usage : "<unknown>"));*/
+    /*printf("gp_alloc %d for %s\n", (int) total_size, usage ? usage : "<unknown>");*/
+    p = malloc(total_size);
+    if (!p){ printf("gp_alloc :  Out of memory !!!! \n"); exit(-99);}
+    /*int_error(NO_CARET, "Out of memory");*/
+
+    bytes_allocated += size;
+
+    /*mark(p, size, usage);*/
+
+    /*return (char *) (p+1);*/
+    return (char *) (p);
+}
+ 
+void gp_free(void *p){
+  free(p);
+}
+
+void contour_free(struct gnuplot_contours *ptr){
+  char *c;
+  struct coordinate *coords;
+  struct gnuplot_contours *old;
+  struct gnuplot_contours *cur;
+ 
+  cur = ptr; 
+  while(cur){
+    old = cur;
+    cur = cur->next;
+    c=old->label;
+    coords = old->coords;
+    if(c){free(c); c=NULL;}
+    if(coords){free(coords) ; coords = NULL;}
+    free(old);    
+  }
+
+}
+/*
+ * iso_alloc() allocates a iso_curve structure that can hold 'num'
+ * points.
+ */
+struct iso_curve *
+iso_alloc(num)
+int num;
+{
+    struct iso_curve *ip;
+    ip = (struct iso_curve *) gp_alloc(sizeof(struct iso_curve), "iso curve");
+    ip->p_max = (num >= 0 ? num : 0);
+    if (num > 0) {
+	ip->points = (struct coordinate  *)
+	    gp_alloc(num * sizeof(struct coordinate), "iso curve points");
+    } else
+	ip->points = (struct coordinate  *) NULL;
+    ip->next = NULL;
+    return (ip);
+}
+/*
+ * iso_free() releases any memory which was previously malloc()'d to hold
+ *   iso curve points.
+ */
+void iso_free(struct iso_curve *ip)
+{
+    if (ip) {
+        if (ip->points)
+            free((char *) ip->points);
+        free((char *) ip);
+    }
+}
Index: /tags/PIGC_V1/SophyaPI/PIGcont/myaxis.c
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/myaxis.c	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/myaxis.c	(revision 1830)
@@ -0,0 +1,97 @@
+#include "myaxis.h"
+#include <math.h>
+AXIS axis_array[AXIS_ARRAY_SIZE] = AXIS_ARRAY_INITIALIZER(DEFAULT_AXIS_STRUCT);
+
+void dummy(){
+
+  int k=4;
+  
+}
+/*{{{  dbl_raise() used by set_tics */
+/* FIXME HBB 20000426: is this really useful? */
+static double
+dbl_raise(x, y)
+double x;
+int y;
+{
+    register int i = abs(y);
+    double val = 1.0;
+
+    while (--i >= 0)
+	val *= x;
+
+    if (y < 0)
+	return (1.0 / val);
+    return (val);
+}
+
+/*}}} */
+
+/*{{{  set_tic() */
+/* the guide parameter was intended to allow the number of tics
+ * to depend on the relative sizes of the plot and the font.
+ * It is the approximate upper limit on number of tics allowed.
+ * But it did not go down well with the users.
+ * A value of 20 gives the same behaviour as 3.5, so that is
+ * hardwired into the calls to here. Maybe we will restore it
+ * to the automatic calculation one day
+ */
+
+double
+set_tic(l10, guide)
+     double l10;
+     int guide;
+{
+    double xnorm, tics, posns;
+
+    int fl = (int) floor(l10);
+    xnorm = pow(10.0, l10 - fl);	/* approx number of decades */
+
+    posns = guide / xnorm;	/* approx number of tic posns per decade */
+
+    if (posns > 40)
+	tics = 0.05;		/* eg 0, .05, .10, ... */
+    else if (posns > 20)
+	tics = 0.1;		/* eg 0, .1, .2, ... */
+    else if (posns > 10)
+	tics = 0.2;		/* eg 0,0.2,0.4,... */
+    else if (posns > 4)
+	tics = 0.5;		/* 0,0.5,1, */
+    else if (posns > 1)
+	tics = 1;		/* 0,1,2,.... */
+    else if (posns > 0.5)
+	tics = 2;		/* 0, 2, 4, 6 */
+    else
+	/* getting desperate... the ceil is to make sure we
+	 * go over rather than under - eg plot [-10:10] x*x
+	 * gives a range of about 99.999 - tics=xnorm gives
+	 * tics at 0, 99.99 and 109.98  - BAD !
+	 * This way, inaccuracy the other way will round
+	 * up (eg 0->100.0001 => tics at 0 and 101
+	 * I think latter is better than former
+	 */
+	tics = ceil(xnorm);
+
+    return (tics * dbl_raise(10.0, fl));
+}
+/* this is used in a few places all over the code: undo logscaling of
+ * a given range if necessary. If checkrange is TRUE, will int_error() if
+ * range is invalid */
+void
+axis_unlog_interval(axis, min, max, checkrange)
+    AXIS_INDEX axis;
+    double *min, *max;
+    TBOOLEAN checkrange;
+{
+    if (axis_array[axis].log) {
+      if (checkrange && (*min<= 0.0 || *max <= 0.0))printf(" axis_unlog_interval : range must be > 0 !! \n");
+	  /*int_error(NO_CARET,
+		      "%s range must be greater than 0 for log scale",
+		      axis_defaults[axis].name);*/
+	*min = (*min<=0) ? -VERYLARGE : AXIS_DO_LOG(axis,*min);
+	*max = (*max<=0) ? -VERYLARGE : AXIS_DO_LOG(axis,*max);
+    }
+} 
+
+/*}}} */
+
Index: /tags/PIGC_V1/SophyaPI/PIGcont/myaxis.h
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/myaxis.h	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/myaxis.h	(revision 1830)
@@ -0,0 +1,615 @@
+/* 
+ * $Id: myaxis.h,v 1.1.1.1 2001-12-19 18:01:18 perderos Exp $
+ *
+ */
+#include <stdio.h>
+/*[
+ * Copyright 2000   Thomas Williams, Colin Kelley
+ *
+ * Permission to use, copy, and distribute this software and its
+ * documentation for any purpose with or 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.
+ *
+ * Permission to modify the software is granted, but not the right to
+ * distribute the complete modified source code.  Modifications are to
+ * be distributed as patches to the released version.  Permission to
+ * distribute binaries produced by compiling modified sources is granted,
+ * provided you
+ *   1. distribute the corresponding source modifications from the
+ *    released version in the form of a patch file along with the binaries,
+ *   2. add special version identification to distinguish your version
+ *    in addition to the base release version number,
+ *   3. provide your name and address as the primary contact for the
+ *    support of your modified version, and
+ *   4. retain our contact information in regard to use of the base
+ *    software.
+ * Permission to distribute the released version of the source code along
+ * with corresponding source modifications in the form of a patch file is
+ * granted with same provisions 2 through 4 for binary distributions.
+ *
+ * This software is provided "as is" without express or implied warranty
+ * to the extent permitted by applicable law.
+]*/
+
+#ifndef GNUPLOT_AXIS_H
+#define GNUPLOT_AXIS_H
+
+#include "mymisc.h"
+/*#include "gp_types.h"	*/	/* for TBOOLEAN */
+
+/*#include "gadgets.h"*/
+/*#include "parse.h"	*/	/* for const_express() */
+/*#include "tables.h"	*/	/* for the axis name parse table */
+/*#include "term_api.h"	*/	/* for lp_style_type */
+/*#include "util.h"	*/	/* for int_error() */
+/****/
+
+/* typedefs / #defines */
+
+/* give some names to some array elements used in command.c and grap*.c
+ * maybe one day the relevant items in setshow will also be stored
+ * in arrays. 
+ *
+ * Always keep the following conditions alive:
+ * SECOND_X_AXIS = FIRST_X_AXIS + SECOND_AXES
+ * FIRST_X_AXIS & SECOND_AXES == 0
+ */
+
+typedef enum AXIS_INDEX {
+#define FIRST_AXES 0
+    FIRST_Z_AXIS, 
+    FIRST_Y_AXIS,
+    FIRST_X_AXIS,
+    T_AXIS,			/* fill gap */
+#define SECOND_AXES 4
+    SECOND_Z_AXIS,		/* not used, yet */
+    SECOND_Y_AXIS,
+    SECOND_X_AXIS,
+    R_AXIS,			/* never used ? */
+    U_AXIS,			/* dito */
+    V_AXIS			/* dito */
+#ifdef PM3D
+    ,COLOR_AXIS
+#endif
+} AXIS_INDEX ;
+
+#ifdef PM3D
+# define AXIS_ARRAY_SIZE 11
+#else
+# define AXIS_ARRAY_SIZE 10
+#endif
+
+
+/* What kind of ticmarking is wanted? */
+typedef enum en_ticseries_type {
+    TIC_COMPUTED=1, 		/* default; gnuplot figures them */
+    TIC_SERIES,			/* user-defined series */
+    TIC_USER,			/* user-defined points */
+    TIC_MONTH,   		/* print out month names ((mo-1)%12)+1 */
+    TIC_DAY      		/* print out day of week */
+} t_ticseries_type;
+
+/* Defines one ticmark for TIC_USER style.
+ * If label==NULL, the value is printed with the usual format string.
+ * else, it is used as the format string (note that it may be a constant
+ * string, like "high" or "low").
+ */
+typedef struct ticmark {
+    double position;		/* where on axis is this */
+    char *label;		/* optional (format) string label */
+    struct ticmark *next;	/* linked list */
+} ticmark;
+
+/* Tic-mark labelling definition; see set xtics */
+typedef struct ticdef {
+    t_ticseries_type type;	
+    union {
+	   struct ticmark *user;	/* for TIC_USER */
+	   struct {			/* for TIC_SERIES */
+		  double start, incr;
+		  double end;		/* ymax, if VERYLARGE */
+	   } series;
+    } def;
+} t_ticdef;
+
+/* we want two auto modes for minitics - default where minitics are
+ * auto for log/time and off for linear, and auto where auto for all
+ * graphs I've done them in this order so that logscale-mode can
+ * simply test bit 0 to see if it must do the minitics automatically.
+ * similarly, conventional plot can test bit 1 to see if minitics are
+ * required */
+enum en_minitics_status {
+    MINI_OFF,
+    MINI_DEFAULT,
+    MINI_USER,
+    MINI_AUTO
+};
+
+/* Function pointer type for callback functions doing operations for a
+ * single ticmark */
+/*typedef void (*tic_callback) __PROTO((AXIS_INDEX, double, char *, struct lp_style_type ));*/
+
+/* Values to put in the axis_tics[] variables that decides where the
+ * ticmarks should be drawn: not at all, on one or both plot borders,
+ * or the zeroaxes. These look like a series of values, but TICS_MASK
+ * shows that they're actually bit masks --> don't turn into an enum
+ * */
+#define NO_TICS        0
+#define TICS_ON_BORDER 1
+#define TICS_ON_AXIS   2
+#define TICS_MASK      3
+#define TICS_MIRROR    4
+
+
+/* Need to allow user to choose grid at first and/or second axes tics.
+ * Also want to let user choose circles at x or y tics for polar grid.
+ * Also want to allow user rectangular grid for polar plot or polar
+ * grid for parametric plot. So just go for full configurability.
+ * These are bitmasks
+ */
+#define GRID_OFF    0
+#define GRID_X      (1<<0)
+#define GRID_Y      (1<<1)
+#define GRID_Z      (1<<2)
+#define GRID_X2     (1<<3)
+#define GRID_Y2     (1<<4)
+#define GRID_MX     (1<<5)
+#define GRID_MY     (1<<6)
+#define GRID_MZ     (1<<7)
+#define GRID_MX2    (1<<8)
+#define GRID_MY2    (1<<9)
+/* GRID_{M}CB only for PM3D, but defined always for source code readability */
+#define GRID_CB     (1<<10)
+#define GRID_MCB    (1<<11)
+
+/* HBB 20010610: new type for storing autoscale activity. Effectively
+ * two booleans (bits) in a single variable, so I'm using an enum with
+ * all 4 possible bit masks given readable names. */
+typedef enum e_autoscale {
+    AUTOSCALE_NONE = 0,
+    AUTOSCALE_MIN = 1<<0,
+    AUTOSCALE_MAX = 1<<1,
+    AUTOSCALE_BOTH = (1<<0 | 1 << 1),
+    AUTOSCALE_FIXMIN = 1<<2,
+    AUTOSCALE_FIXMAX = 1<<3
+} t_autoscale;
+
+
+/* HBB 20000725: gather all per-axis variables into a struct, and set up
+ * a single large array of such structs */
+/* FIXME 20000725: collect some of those various TBOOLEAN fields into
+ * a larger int (or -- shudder -- a bitfield?) */
+typedef struct axis {
+/* range of this axis */
+    t_autoscale autoscale;	/* Which end(s) are autoscaled? */
+    t_autoscale set_autoscale;	/* what does 'set' think autoscale to be? */
+    int range_flags;		/* flag bits about autoscale/writeback: */
+    /* write auto-ed ranges back to variables for autoscale */
+#define RANGE_WRITEBACK 1	
+    /* allow auto and reversed ranges */
+#define RANGE_REVERSE   2	
+    TBOOLEAN reverse_range;	/* range [high:low] silently reverted? */
+    double min;			/* 'transient' axis extremal values */
+    double max;
+    double set_min;		/* set/show 'permanent' values */
+    double set_max;
+    double writeback_min;	/* ULIG's writeback implementation */
+    double writeback_max;
+
+/* output-related quantities */
+    int term_lower;		/* low and high end of the axis on output, */
+    int term_upper;		/* ... (in terminal coordinates)*/
+    double term_scale;		/* scale factor: plot --> term coords */
+    unsigned int term_zero;	/* position of zero axis */
+
+/* log axis control */
+    TBOOLEAN log;		/* log axis stuff: flag "islog?" */
+    double base;		/* logarithm base value */
+    double log_base;		/* ln(base), for easier computations */
+
+/* time/date axis control */
+    TBOOLEAN is_timedata;	/* is this a time/date axis? */
+    TBOOLEAN format_is_numeric;	/* format string looks like numeric??? */
+    char timefmt[MAX_ID_LEN+1];	/* format string for input */
+    char formatstring[MAX_ID_LEN+1];
+				/* the format string for output */
+
+/* ticmark control variables */
+    int ticmode;		/* tics on border/axis? mirrored? */
+    struct ticdef ticdef;	/* tic series definition */
+    TBOOLEAN tic_rotate;	/* ticmarks rotated (vertical text)? */
+    int minitics;		/* minor tic mode (none/auto/user)? */
+    double mtic_freq;		/* minitic stepsize */
+
+/* other miscellaneous fields */
+    label_struct label;		/* label string and position offsets */
+  /*lp_style_type zeroaxis;*/	/* drawing style for zeroaxis, if any */
+} AXIS;
+
+#define DEFAULT_AXIS_TICDEF {TIC_COMPUTED, {NULL} }
+#ifdef PM3D
+# define DEFAULT_AXIS_ZEROAXIS {0, -3, 0, 1.0, 1.0, 0}
+#else
+# define DEFAULT_AXIS_ZEROAXIS {0, -3, 0, 1.0, 1.0}
+#endif
+#define DEFAULT_AXIS_STRUCT {						    \
+	AUTOSCALE_BOTH, AUTOSCALE_BOTH, /* auto, set_auto */		    \
+	0, FALSE,		/* range_flags, rev_range */		    \
+	-10.0, 10.0,		/* 3 pairs of min/max */		    \
+	-10.0, 10.0,							    \
+	-10.0, 10.0,							    \
+	0, 0, 0, 0,		/* terminal dependents */		    \
+	FALSE, 0.0, 0.0,	/* log, base, log(base) */		    \
+	0, 1,			/* is_timedata, format_numeric */	    \
+	DEF_FORMAT, TIMEFMT,	/* output format, timefmt */		    \
+	NO_TICS,		/* tic output positions (border, mirror) */ \
+	DEFAULT_AXIS_TICDEF,	/* tic series definition */		    \
+	FALSE, MINI_DEFAULT, 10, /* tic_rotate, minitics, mtic_freq */	    \
+	EMPTY_LABELSTRUCT/**,**/	/* axis label */			    \
+	/**DEFAULT_AXIS_ZEROAXIS**/}	/* zeroaxis line style */
+
+/* Table of default behaviours --- a subset of the struct above. Only
+ * those fields are present that differ from axis to axis. */
+typedef struct axis_defaults {
+    double min;			/* default axis endpoints */
+    double max;
+    char name[4];		/* axis name, like in "x2" or "t" */
+    int ticmode;		/* tics on border/axis? mirrored? */
+} AXIS_DEFAULTS;
+
+
+
+/* global variables in axis.c */
+
+extern AXIS axis_array[AXIS_ARRAY_SIZE];
+extern AXIS_DEFAULTS axis_defaults[AXIS_ARRAY_SIZE];
+
+/* A parsing table for mapping axis names into axis indices. For use
+ * by the set/show machinery, mainly */
+extern struct gen_table axisname_tbl[AXIS_ARRAY_SIZE+1];
+
+
+extern const struct ticdef default_axis_ticdef;
+
+extern double ticscale;		/* scale factor for tic marks (was (0..1])*/
+extern double miniticscale;	/* and for minitics */
+extern TBOOLEAN	tic_in;		/* tics to be drawn inward?  */
+
+
+/* default format for tic mark labels */
+#define DEF_FORMAT "% g"
+
+/* default parse timedata string */
+#define TIMEFMT "%d/%m/%y,%H:%M"
+
+/* axis labels */
+extern const label_struct default_axis_label;
+
+/* zeroaxis linetype (flag type==-3 if none wanted) */
+
+/***extern const lp_style_type default_axis_zeroaxis;***/
+
+/* grid drawing control variables */
+extern int grid_selection;
+/***extern const struct lp_style_type default_grid_lp;***/ /* default for 'unset' */
+/***extern struct lp_style_type grid_lp;***/ /* linestyle for major grid lines */
+/***extern struct lp_style_type mgrid_lp;***/ /* linestyle for minor grid lines */
+/***extern double polar_grid_angle;***/ /* angle step in polar grid in radians */
+
+/* global variables for communication with the tic callback functions */
+
+extern int tic_start, tic_direction, tic_text,
+    rotate_tics, tic_hjust, tic_vjust, tic_mirror;
+
+/* axes being used by the current plot */
+extern AXIS_INDEX x_axis, y_axis, z_axis;
+/* macros to reduce code clutter caused by the array notation, mainly
+ * in graphics.c and fit.c */
+#define X_AXIS axis_array[x_axis]
+#define Y_AXIS axis_array[y_axis]
+#define Z_AXIS axis_array[z_axis]
+#ifdef PM3D
+#define CB_AXIS axis_array[COLOR_AXIS]
+#endif
+
+/* -------- macros using these variables: */
+
+/* Macros to map from user to terminal coordinates and back */
+#define AXIS_MAP(axis, variable)		\
+  (int) ((axis_array[axis].term_lower)		\
+	 + ((variable) - axis_array[axis].min)	\
+	 * axis_array[axis].term_scale + 0.5)
+#define AXIS_MAPBACK(axis, pos)						   \
+  (((double)(pos)-axis_array[axis].term_lower)/axis_array[axis].term_scale \
+   + axis_array[axis].min)
+
+/* these are the old names for these: */
+#define map_x(x) AXIS_MAP(x_axis, x)
+#define map_y(y) AXIS_MAP(y_axis, y)
+
+#define AXIS_SETSCALE(axis, out_low, out_high)			\
+    axis_array[axis].term_scale = ((out_high) - (out_low))	\
+        / (axis_array[axis].max - axis_array[axis].min)
+
+/* write current min/max_array contents into the set/show status
+ * variables */
+#define AXIS_WRITEBACK(axis)			\
+do {						\
+    AXIS *this = axis_array + axis;		\
+						\
+    if (this->range_flags & RANGE_WRITEBACK) {	\
+	if (this->autoscale & AUTOSCALE_MIN)	\
+	    this->set_min = this->min;		\
+	if (this->autoscale & AUTOSCALE_MAX)	\
+	    this->set_max = this->max;		\
+    }						\
+} while(0)
+
+/* HBB 20000430: New macros, logarithmize a value into a stored
+ * coordinate*/ 
+#define AXIS_DO_LOG(axis,value) (log(value) / axis_array[axis].log_base)
+#define AXIS_UNDO_LOG(axis,value) exp((value) * axis_array[axis].log_base)
+
+/* HBB 20000430: same, but these test if the axis is log, first: */
+#define AXIS_LOG_VALUE(axis,value)				\
+    (axis_array[axis].log ? AXIS_DO_LOG(axis,value) : (value))
+#define AXIS_DE_LOG_VALUE(axis,coordinate)				  \
+    (axis_array[axis].log ? AXIS_UNDO_LOG(axis,coordinate): (coordinate))
+ 
+
+/* copy scalar data to arrays. The difference between 3D and 2D
+ * versions is: dont know we have to support ranges [10:-10] - lets
+ * reverse it for now, then fix it at the end.  */
+/* FIXME HBB 20000426: unknown if this distinction makes any sense... */
+#define AXIS_INIT3D(axis, islog_override, infinite)			\
+do {									\
+    AXIS *this = axis_array + axis;					\
+									\
+    if ((this->autoscale = this->set_autoscale) == AUTOSCALE_NONE	\
+	&& this->set_max < this->set_min) {				\
+	this->min = this->set_max;					\
+	this->max = this->set_min;					\
+        /* we will fix later */						\
+    } else {								\
+	this->min = (infinite && (this->set_autoscale & AUTOSCALE_MIN))	\
+	    ? VERYLARGE : this->set_min;				\
+	this->max = (infinite && (this->set_autoscale & AUTOSCALE_MAX))	\
+	    ? -VERYLARGE : this->set_max;				\
+    }									\
+    if (islog_override) {						\
+	this->log = 0;							\
+	this->base = 1;							\
+	this->log_base = 0;						\
+    } else {								\
+	this->log_base = this->log ? log(this->base) : 0;		\
+    }									\
+} while(0)
+
+#define AXIS_INIT2D(axis, infinite)					\
+do {									\
+    AXIS *this = axis_array + axis;					\
+									\
+    this->autoscale = this->set_autoscale;				\
+    this->min = (infinite && (this->set_autoscale & AUTOSCALE_MIN))	\
+	? VERYLARGE : this->set_min;					\
+    this->max = (infinite && (this->set_autoscale & AUTOSCALE_MAX))	\
+	? -VERYLARGE : this->set_max;					\
+    this->log_base = this->log ? log(this->base) : 0;			\
+} while(0)
+
+/* handle reversed ranges */
+#define CHECK_REVERSE(axis)						\
+do {									\
+    AXIS *this = axis_array + axis;					\
+									\
+    if ((this->autoscale == AUTOSCALE_NONE)				\
+	&& (this->max < this->min)					\
+	) {								\
+	double temp = this->min;					\
+	this->min = this->max;						\
+	this->max = temp;						\
+	this->reverse_range = 1;					\
+    } else								\
+	this->reverse_range = (this->range_flags & RANGE_REVERSE);	\
+} while(0)
+
+/* HBB 20000725: new macro, built upon ULIG's SAVE_WRITEBACK(axis),
+ * but easier to use. Code like this occured twice, in plot2d and
+ * plot3d: */
+#define SAVE_WRITEBACK_ALL_AXES					\
+do {								\
+    AXIS_INDEX axis;						\
+								\
+    for (axis = 0; axis < AXIS_ARRAY_SIZE; axis++)		\
+	if(axis_array[axis].range_flags & RANGE_WRITEBACK) {	\
+	    set_writeback_min(axis);				\
+	    set_writeback_max(axis);				\
+	}							\
+} while(0)
+	    
+/* get optional [min:max] */
+#define PARSE_RANGE(axis)						   \
+do {									   \
+    if (equals(c_token, "[")) {						   \
+	c_token++;							   \
+	axis_array[axis].autoscale =					   \
+	    load_range(axis, &axis_array[axis].min, &axis_array[axis].max, \
+		       axis_array[axis].autoscale);			   \
+	if (!equals(c_token, "]"))					   \
+	    int_error(c_token, "']' expected");				   \
+	c_token++;							   \
+    }									   \
+} while (0)
+
+/* HBB 20000430: new macro, like PARSE_RANGE, but for named ranges as
+ * in 'plot [phi=3.5:7] sin(phi)' */
+#define PARSE_NAMED_RANGE(axis, dummy_token)				     \
+do {									     \
+    if (equals(c_token, "[")) {						     \
+	c_token++;							     \
+	if (isletter(c_token)) {					     \
+	    if (equals(c_token + 1, "=")) {				     \
+		dummy_token = c_token;					     \
+		c_token += 2;						     \
+	    }								     \
+		/* oops; probably an expression with a variable: act	     \
+		 * as if no variable name had been seen, by		     \
+		 * fallthrough */					     \
+	}								     \
+	axis_array[axis].autoscale = load_range(axis, &axis_array[axis].min, \
+				      &axis_array[axis].max,		     \
+				      axis_array[axis].autoscale);	     \
+	if (!equals(c_token, "]"))					     \
+	    int_error(c_token, "']' expected");				     \
+	c_token++;							     \
+    }				/* first '[' */				     \
+} while (0)
+
+/* parse a position of the form
+ *    [coords] x, [coords] y {,[coords] z}
+ * where coords is one of first,second.graph,screen
+ * if first or second, we need to take axis_is_timedata into account
+ */
+#define GET_NUMBER_OR_TIME(store,axes,axis)				\
+do {									\
+    if (((axes) >= 0) && (axis_array[(axes)+(axis)].is_timedata)	\
+	&& isstring(c_token)) {						\
+	char ss[80];							\
+	struct tm tm;							\
+	quote_str(ss,c_token, 80);					\
+	++c_token;							\
+	if (gstrptime(ss,axis_array[axis].timefmt,&tm))			\
+	    (store) = (double) gtimegm(&tm);				\
+    } else {								\
+	struct value value;						\
+	(store) = real(const_express(&value));				\
+    }									\
+} while(0)
+
+/* This is one is very similar to GET_NUMBER_OR_TIME, but has slightly
+ * different usage: it writes out '0' in case of inparsable time data,
+ * and it's used when the target axis is fixed without a 'first' or
+ * 'second' keyword in front of it. */
+#define GET_NUM_OR_TIME(store,axis)			\
+do {							\
+    (store) = 0;					\
+    GET_NUMBER_OR_TIME(store, FIRST_AXES, axis);	\
+} while (0);
+
+/* store VALUE or log(VALUE) in STORE, set TYPE as appropriate
+ * Do OUT_ACTION or UNDEF_ACTION as appropriate
+ * adjust range provided type is INRANGE (ie dont adjust y if x is outrange
+ * VALUE must not be same as STORE
+ */
+
+#define STORE_WITH_LOG_AND_UPDATE_RANGE(STORE, VALUE, TYPE, AXIS,	  \
+				       OUT_ACTION, UNDEF_ACTION)	  \
+do {									  \
+    /* HBB 20000726: new check, to avoid crashes with axis index -1 */	  \
+    if (AXIS==-1)							  \
+	break;								  \
+    if (axis_array[AXIS].log) {						  \
+	if (VALUE<0.0) {						  \
+	    TYPE = UNDEFINED;						  \
+	    UNDEF_ACTION;						  \
+	    break;							  \
+	} else if (VALUE == 0.0) {					  \
+	    STORE = -VERYLARGE;						  \
+	    TYPE = OUTRANGE;						  \
+	    OUT_ACTION;							  \
+	    break;							  \
+	} else {							  \
+	    STORE = AXIS_DO_LOG(AXIS,VALUE);				  \
+	}								  \
+    } else								  \
+	STORE = VALUE;							  \
+    if (TYPE != INRANGE)						  \
+	break;  /* don't set y range if x is outrange, for example */	  \
+    if ((int)AXIS < 0)							  \
+	break;	/* HBB 20000507: don't check range if not a coordinate */ \
+    if ( VALUE<axis_array[AXIS].min ) {					  \
+	if (axis_array[AXIS].autoscale & AUTOSCALE_MIN)			  \
+	    axis_array[AXIS].min = VALUE;				  \
+	else {								  \
+	    TYPE = OUTRANGE;						  \
+	    OUT_ACTION;							  \
+	    break;							  \
+	}								  \
+    }									  \
+    if ( VALUE>axis_array[AXIS].max ) {					  \
+	if (axis_array[AXIS].autoscale & AUTOSCALE_MAX)			  \
+	    axis_array[AXIS].max = VALUE;				  \
+	else {								  \
+	    TYPE = OUTRANGE;						  \
+	    OUT_ACTION;							  \
+	}								  \
+    }									  \
+} while(0)
+
+/* use this instead empty macro arguments to work around NeXT cpp bug */
+/* if this fails on any system, we might use ((void)0) */
+#define NOOP			/* */
+
+/* HBB 20000506: new macro, initializes one variable to the same
+ * value, for all axes. */
+#define INIT_AXIS_ARRAY(field, value)		\
+do {						\
+    int tmp;					\
+    for (tmp=0; tmp<AXIS_ARRAY_SIZE; tmp++)	\
+	axis_array[tmp].field=(value);		\
+} while(0)
+
+/* HBB 20000506: new macro to automatically build intializer lists
+ * for arrays of AXIS_ARRAY_SIZE equal elements */
+#ifdef PM3D
+#define AXIS_ARRAY_INITIALIZER(value) {			\
+    value, value, value, value, value,			\
+	value, value, value, value, value /*, value*/ }
+#else
+#define AXIS_ARRAY_INITIALIZER(value) {		\
+    value, value, value, value, value,		\
+	value, value, value, value, value }
+#endif
+
+
+/* used by set.c */
+#define SET_DEFFORMAT(axis, flag_array)				\
+	if (flag_array[axis]) {					\
+	    (void) strcpy(axis_array[axis].formatstring,DEF_FORMAT);	\
+	    axis_array[axis].format_is_numeric = 1;		\
+	}
+
+
+/* 'roundoff' check tolerance: less than one hundredth of a tic mark */
+#define SIGNIF (0.01)
+/* (DFK) Watch for cancellation error near zero on axes labels */
+/* FIXME HBB 20000521: these seem not to be used much, anywhere... */
+#define CheckZero(x,tic) (fabs(x) < ((tic) * SIGNIF) ? 0.0 : (x))
+#define NearlyEqual(x,y,tic) (fabs((x)-(y)) < ((tic) * SIGNIF))
+
+
+	
+/* ------------ functions exported by axis.c */
+/*t_autoscale load_range __PROTO((AXIS_INDEX, double *, double *, t_autoscale));*/
+void axis_unlog_interval(AXIS_INDEX, double *, double *, TBOOLEAN);
+/*void axis_revert_and_unlog_range __PROTO((AXIS_INDEX));*/
+/*double axis_log_value_checked __PROTO((AXIS_INDEX, double, const char *));*/
+/*void axis_checked_extend_empty_range __PROTO((AXIS_INDEX, const char *mesg));*/
+/*void   timetic_format __PROTO((AXIS_INDEX, double, double));*/
+double set_tic(double, int); 
+/*void setup_tics __PROTO((AXIS_INDEX, int));*/
+/*void gen_tics __PROTO((AXIS_INDEX, int, tic_callback));*/
+/*void gprintf __PROTO((char *, size_t, char *, double, double));*/
+/*void axis_output_tics __PROTO((AXIS_INDEX, int *, AXIS_INDEX, int, tic_callback));*/
+/*void axis_set_graphical_range __PROTO((AXIS_INDEX, unsigned int lower, unsigned int upper));*/
+/*TBOOLEAN axis_position_zeroaxis __PROTO((AXIS_INDEX));*/
+/*void axis_draw_2d_zeroaxis __PROTO((AXIS_INDEX, AXIS_INDEX));
+
+/*double get_writeback_min __PROTO((AXIS_INDEX));*/
+/*double get_writeback_max __PROTO((AXIS_INDEX));*/
+/*void set_writeback_min __PROTO((AXIS_INDEX));*/
+/*void set_writeback_max __PROTO((AXIS_INDEX));*/
+
+#endif /* GNUPLOT_AXIS_H */
Index: /tags/PIGC_V1/SophyaPI/PIGcont/mymisc.h
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/mymisc.h	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/mymisc.h	(revision 1830)
@@ -0,0 +1,43 @@
+#ifndef MYMISC_H
+#define MYMISC_H
+
+#ifndef __cplusplus
+#undef bool
+typedef unsigned int bool;
+#endif
+
+/* TRUE or FALSE */
+#define TBOOLEAN bool
+#define EMPTY_LABELSTRUCT {"", "", 0.0, 0.0}
+#define NO_CARET (-1)
+#define VERYLARGE 1e37
+
+#define MAX_LINE_LEN 1024       /* maximum number of chars allowed on line */
+#define MAX_ID_LEN 50  
+
+#define ZERO 1e-8               /* default for 'zero' set option */
+extern double zero;             /* zero threshold, not 0! */
+
+/** labels **/
+/* The data describing an axis label, or the plot title */
+typedef struct {
+    char text[MAX_LINE_LEN+1];
+    char font[MAX_LINE_LEN+1];
+    double xoffset, yoffset;
+} label_struct;
+
+typedef enum position_type {
+    first_axes,
+    second_axes,
+    graph,
+    screen
+} position_type; 
+
+/* A full 3D position, with all 3 coordinates of different axes,
+ * possibly. Used for 'set label' and 'set arrow' positions: */
+typedef struct position {
+    enum position_type scalex,scaley,scalez;
+    double x,y,z;
+} t_position;
+
+#endif /*ifndef MYMISC_H*/
Index: /tags/PIGC_V1/SophyaPI/PIGcont/picntools.cc
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/picntools.cc	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/picntools.cc	(revision 1830)
@@ -0,0 +1,775 @@
+// This may look like C code, but it is really -*- C++ -*-
+// Outils de manipulation des Contours  O. Perdereau 11/2001
+// LAL (Orsay) / IN2P3-CNRS  DAPNIA/SPP (Saclay) / CEA
+
+#include <stdlib.h>
+
+
+#include "picntools.h"
+#include <iostream.h>
+#include "piscdrawwdg.h"
+#include "pi3ddrw.h"
+#include "pigncont.h"
+#include "perrors.h"
+#include "exceptions.h"
+
+//++
+// Class	PICnTools
+// Lib		PI
+// include	pidrwtools.h
+//
+//	Classe de fenêtre de dialogue permettant de modifier interactivement
+//	Les attributs graphique des "Drawer" attachés aux "PIBaseWdg".
+//--
+// Links	Parents
+// PIWindow
+//--
+//++
+// Links	Voir aussi
+// PIDrawer
+// PIBaseWdg
+//--
+//++
+// Links	Ainsi que
+// PIScDrawWdg
+// PIImage
+// PIDraw3DWdg
+//--
+
+//++
+// Titre	Méthodes statiques
+//--
+//++
+// void  ShowPICnTools() 	
+//	Méthode statique permettant l'affichage de la fenêtre de contrôle des
+//	attributs graphiques des "Drawer". L'instance de la classe
+//	 "PICnTools" est créée si nécessaire.
+// void  HidePICnTools() 
+//	Désacative (cache) la fenêtre "PICnTools" (méthode statique).
+// void  SetCurrentBaseWdg(PIBaseWdgGen* cbw)
+//	Méthode statique permettant de sélectionner le "PIBaseWdg" courant.
+// PIBaseWdgGen*  GetCurrentBaseWdg()
+//	Retourne le "PIBaseWdg" courant (méthode statique).
+//--
+
+
+PICnTools* PICnTools::cwdrwt = NULL;
+PIBaseWdgGen* PICnTools::mCurBW = NULL;
+PIContourDrawer* PICnTools::mCurCDrw = NULL;
+
+/* --Methode-- */
+void PICnTools::ShowPICnTools() 
+{
+if (cwdrwt == NULL)  cwdrwt = new PICnTools(PIApplicationGetApp());
+// cout << " from PICnTools::ShowPICnTools() - hello "<<endl;
+cwdrwt->Show();
+}
+
+/* --Methode-- */
+void PICnTools::HidePICnTools() 
+{
+if (cwdrwt == NULL)  return; 
+cwdrwt->Hide();
+}
+
+void PICnTools::SetCurrentBaseWdg(PIBaseWdgGen* cbw)
+{
+if (mCurBW == cbw) return;
+mCurBW = cbw;
+}
+
+/* --Methode-- */
+PIBaseWdgGen* PICnTools::GetCurrentBaseWdg()
+{
+return(mCurBW);
+}
+/* --Methode-- */
+void PICnTools::SetCurrentCnDrw( PIContourDrawer* pcd) 
+{ 
+ 
+    mCurCDrw = pcd ; 
+}
+
+/* --Methode-- */
+PIContourDrawer*  PICnTools::GetCurrentCnDrw()
+{
+return(mCurCDrw);
+}
+
+/* --Methode-- */
+PICnTools::PICnTools(PIApplication* par)
+: PIWindow((PIMsgHandler*)par, "PIContourTools", PIWK_normal, 240, 240, 150, 150)
+{
+int bsx, bsy, spx, spy;
+
+// On definit la taille a partir de la taille par defaut des composantes 
+// PIApplicationPrefCompSize(bsx, bsy);
+par->PrefCompSize(bsx, bsy);
+
+// Espacements entre elts 
+spx = bsx/10;
+spy = bsy/4;
+
+int wszx = 5*spx+5.5*bsx; 
+int wszy = (1+8*0.85)*bsy+12*(1.85)*spy+15;
+SetSize(wszx, wszy);
+
+//int cpx = spx*2;
+//int offy = cpy+bsy+spy*2;
+
+//  Option-Menu  Marker, Line, couleur, etc...
+int cpx = 2*spx; // current position 
+int cpy = 2*spy;
+
+int bsyo = bsy;  
+int spyo = spy; 
+bsy = 0.85*bsy;   spy = 0.85*spy;
+
+mLab[0] = new PILabel(this, " Couleurs des  contours ", 2*bsx, 1.25*bsy, cpx, cpy);
+mLab[0]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
+
+
+cpx += spx+2*bsx; 
+mOpt[1] = new PIOptMenu(this, "cwt-opt-1", bsx*1.75, bsy, cpx, cpy);
+mOpt[1]->AppendItem("No.CMap", 200);
+int kcc;
+//  Groupement des tables de couleurs
+mCasc[0] = new PIMenu(mOpt[1]->Menu(), "PIStd-128Col");
+mCasc[1] = new PIMenu(mOpt[1]->Menu(), "MIDAS-CMap");
+
+int nsct1,nsct2,nsct3;
+//  D'abord les tables standard de PI  32 couleurs
+nsct1 = 5;   // Les 5 premieres tables
+for(kcc=0; kcc<nsct1; kcc++) 
+  mOpt[1]->AppendItem(PIColorMap::GetStandardColorMapName(kcc).c_str(), 201+kcc);
+//  Ensuite les tables standard de PI  128 couleurs
+nsct2 = 9;   // Les 4 tables suivantes
+for(kcc=nsct1; kcc<nsct2; kcc++) 
+  mCasc[0]->AppendItem(PIColorMap::GetStandardColorMapName(kcc).c_str(), 201+kcc);
+mOpt[1]->AppendPDMenu(mCasc[0]);
+//  Apres les tables importees de MIDAS
+nsct3 = PIColorMap::NumberStandardColorMaps()-1; // Les reste jusqu'a l'avant derniere
+for(kcc=nsct2; kcc<nsct3; kcc++) 
+  mCasc[1]->AppendItem(PIColorMap::GetStandardColorMapName(kcc).c_str(), 201+kcc);
+
+mOpt[1]->AppendPDMenu(mCasc[1]);
+// Les tables qui restent ( Col16 , ... )
+for(kcc=nsct3; kcc<PIColorMap::NumberStandardColorMaps(); kcc++) 
+mOpt[1]->AppendItem(PIColorMap::GetStandardColorMapName(kcc).c_str(), 201+kcc);
+
+mOpt[1]->SetValue(200);
+mOpt[1]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+cpx += spx+1.75*bsx;
+
+mOpt[0] = new PIOptMenu(this, "cwt-opt-0", bsx*1.75, bsy, cpx, cpy);
+mOpt[0]->AppendItem("Def.Col", 100);
+mOpt[0]->AppendItem("Black", 101);
+mOpt[0]->AppendItem("White", 102);
+mOpt[0]->AppendItem("Grey", 103);
+mOpt[0]->AppendItem("Red", 104);
+mOpt[0]->AppendItem("Blue", 105);
+mOpt[0]->AppendItem("Green", 106);
+mOpt[0]->AppendItem("Yellow", 107);
+mOpt[0]->AppendItem("Magenta", 108);
+mOpt[0]->AppendItem("Cyan", 109);
+mOpt[0]->AppendItem("Turquoise", 110);
+mOpt[0]->AppendItem("NavyBlue", 111);
+mOpt[0]->AppendItem("Orange", 112);
+mOpt[0]->AppendItem("SiennaRed", 113);
+mOpt[0]->AppendItem("Purple", 114);
+mOpt[0]->AppendItem("LimeGreen", 115);
+mOpt[0]->AppendItem("Gold", 116);
+mOpt[0]->AppendItem("Violet", 117);
+mOpt[0]->AppendItem("VioletRed", 118);
+mOpt[0]->AppendItem("BlueViolet", 119);
+mOpt[0]->AppendItem("DarkViolet", 120);
+mOpt[0]->SetValue(100);
+mOpt[0]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+
+
+cpx = 2*spx;
+cpy += spy+bsyo;
+
+cmapv = new PICMapView(this, "CMapView",wszx-4*spx,20 ,cpx,cpy); 
+cmapv->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+
+
+
+cpx = 2*spx; // current position 
+cpy += spy+bsy;
+
+// check boxes : lines on/off ; marker on/off 
+mCkb[0] = new PICheckBox(this,"LineON", 2001, bsx, bsy, cpx, cpy);
+mCkb[0]->SetState(false); 
+mCkb[0]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
+
+// type de ligne pour les niveaux
+cpx += spx+1.75*bsx; 
+mOpt[5] = new PIOptMenu(this, "bwt-opt-5", bsx*1.75, bsy, cpx, cpy);
+mOpt[5]->AppendItem("Def.Line", 600);
+mOpt[5]->AppendItem("Thin-Line", 601);
+mOpt[5]->AppendItem("Normal-Line", 602);
+mOpt[5]->AppendItem("Thick-Line", 603);
+mOpt[5]->AppendItem("ThinDashed", 604);
+mOpt[5]->AppendItem("Dashed-Line", 605);
+mOpt[5]->AppendItem("ThickDashed", 606);
+mOpt[5]->AppendItem("ThinDotted", 607);
+mOpt[5]->AppendItem("Dotted-Line", 608);
+mOpt[5]->AppendItem("ThickDotted", 609);
+mOpt[5]->SetValue(600);
+mOpt[5]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+
+if(mCkb[0]->GetState()==false) 
+  mOpt[5]->SetUnSensitive();
+
+
+cpx = 2*spx;
+cpy += spy+bsy;
+//cpx += spx+bsx; 
+
+mCkb[1] = new PICheckBox(this,"MarkerON", 2002, bsx, bsy, cpx, cpy);
+mCkb[1]->SetState(true); 
+mCkb[1]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
+cpx += spx+1.75*bsx; 
+
+mOpt[2] = new PIOptMenu(this, "bwt-opt-2", bsx*1.75, bsy, cpx, cpy);
+mOpt[2]->AppendItem("Def.Mrk", 300);
+mOpt[2]->AppendItem("Mrk Dot(.)", 301);
+mOpt[2]->AppendItem("Mrk Plus(+)", 302);
+mOpt[2]->AppendItem("Mrk Cross(x)", 303);
+mOpt[2]->AppendItem("Mrk Circle", 304);
+mOpt[2]->AppendItem("Mrk FCircle", 305);
+mOpt[2]->AppendItem("Mrk Box", 306);
+mOpt[2]->AppendItem("Mrk FBox", 307);
+mOpt[2]->AppendItem("Mrk Triangle", 308);
+mOpt[2]->AppendItem("Mrk FTriangle", 309);
+mOpt[2]->AppendItem("Mrk Star", 310);
+mOpt[2]->AppendItem("Mrk FStar", 311);
+mOpt[2]->SetValue(300);
+mOpt[2]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+if(mCkb[1]->GetState()==false) 
+  mOpt[2]->SetUnSensitive();
+
+cpx += spx+1.75*bsx; 
+mOpt[3] = new PIOptMenu(this, "bwt-opt-3", bsx*1.75, bsy, cpx, cpy);
+mOpt[3]->AppendItem("Mrk 1 pt", 401);
+mOpt[3]->AppendItem("Mrk 3 pt", 403);
+mOpt[3]->AppendItem("Mrk 5 pt", 405);
+mOpt[3]->AppendItem("Mrk 7 pt", 407);
+mOpt[3]->AppendItem("Mrk 9 pt", 409);
+mOpt[3]->AppendItem("Mrk 11 pt", 411);
+mOpt[3]->AppendItem("Mrk 13 pt", 413);
+mOpt[3]->AppendItem("Mrk 15 pt", 415);
+mOpt[3]->SetValue(405);
+mOpt[3]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+if(mCkb[1]->GetState()==false) 
+  mOpt[3]->SetUnSensitive();
+
+
+
+cpx = 2*spx;
+cpy += spy+bsy;
+mLab[1] = new PILabel(this, " Label des  contours ", 2*bsx, 1.25*bsy, cpx, cpy);
+mLab[1]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
+cpx += spx+2*bsx; 
+//  check boxes : label on/off 
+mCkb[2] = new PICheckBox(this,"LabelON", 2003, bsx, bsy, cpx, cpy);
+mCkb[2]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
+mCkb[2]->SetState(false);
+
+// fonte pour les niveaux 
+cpx = 2*spx;
+cpy += spy+bsy;
+mOpt[4] = new PIOptMenu(this, "bwt-opt-4", bsx*1.75, bsy, cpx, cpy);
+mOpt[4]->AppendItem("Def.Font", 500);
+mOpt[4]->AppendItem("Courier", 501);
+mOpt[4]->AppendItem("Helvetica", 502);
+mOpt[4]->AppendItem("Times", 503);
+mOpt[4]->AppendItem("Symbol-Font", 504);
+mOpt[4]->SetValue(500);
+mOpt[4]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+if(mCkb[2]->GetState()==false) mOpt[4]->SetUnSensitive();
+
+
+cpx += spx+1.75*bsx;
+mOpt[6] = new PIOptMenu(this, "bwt-opt-6", bsx*1.75, bsy, cpx, cpy);
+mOpt[6]->AppendItem("Def.FontAtt", 700);
+mOpt[6]->AppendItem("Roman-Font", 701);
+mOpt[6]->AppendItem("Bold-Font", 702);
+mOpt[6]->AppendItem("Italic-Font", 703);
+mOpt[6]->AppendItem("BoldItalic", 704);
+mOpt[6]->SetValue(700);
+mOpt[6]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+if(mCkb[2]->GetState()==false) mOpt[6]->SetUnSensitive();
+
+cpx += spx+1.75*bsx; 
+mOpt[7] = new PIOptMenu(this, "bwt-opt-7", bsx*1.75, bsy, cpx, cpy);
+mOpt[7]->AppendItem("Def.FontSz", 800);
+mOpt[7]->AppendItem("Tiny FontSz", 801);
+mOpt[7]->AppendItem("Small FontSz", 802);
+mOpt[7]->AppendItem("Normal FontSz", 803);
+mOpt[7]->AppendItem("Big FontSz", 804);
+mOpt[7]->AppendItem("Huge FontSz", 805);
+mOpt[7]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+if(mCkb[2]->GetState()==false) mOpt[7]->SetUnSensitive();
+
+cpy += spyo+bsyo;
+cpx = 2*spx;
+mButcf[0] = new PIButton(this, "GetAtt", 2800, bsx, bsy, cpx, cpy);
+mButcf[0]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+
+cpx += bsx+spx;
+mButcf[1] = new PIButton(this, "SetAtt", 2900, bsx, bsy, cpx, cpy);
+mButcf[1]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+
+cpx += bsx+spx;
+mButcf[4] = new PIButton(this, "Refresh", 2950, bsx, bsy, cpx, cpy);
+mButcf[4]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+
+
+cpx = 2*spx;
+cpy += spy+bsyo;
+mLab[2] = new PILabel(this, " Options de calcul des contours ", 3*bsx, bsy, cpx, cpy);
+mLab[2]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+
+cpy += spy+bsyo;
+
+mOpt[8] = new PIOptMenu(this, "bwt-opt-8", bsx*1.75, bsy, cpx, cpy);
+mOpt[8]->AppendItem("Int. Lin.", 900);
+mOpt[8]->AppendItem("CubicCpl", 901);
+mOpt[8]->AppendItem("B-Spline", 902);
+mOpt[8]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+
+
+cpx += spx+1.75*bsx;
+mOpt[9] = new PIOptMenu(this, "bwt-opt-9", bsx*1.75, bsy, cpx, cpy);
+mOpt[9]->AppendItem("LevelAuto", 910);
+mOpt[9]->AppendItem("LevelNum ", 911);
+mOpt[9]->AppendItem("LevelInc ", 912);
+mOpt[9]->AppendItem("LevelDisc", 913);
+mOpt[9]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+
+cpy += spy+bsyo;
+cpx = 2*spx;
+
+mLab[3] = new PILabel(this, "# Contours", bsx, bsy, cpx, cpy);
+mLab[3]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+cpx += spx+bsx;
+
+mTxt[0] = new PIText(this, "Nlev",bsx, (int)(bsy*1.3), cpx, cpy);
+mTxt[0]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+mTxt[0]->SetText("5");
+cpx += spx+bsx;
+mLab[4] = new PILabel(this, "Niveaux", bsx, bsy, cpx, cpy);
+mLab[4]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+cpx += spx+bsx;
+mTxt[1] = new PIText(this, "levels",(int)(bsx*2), (int)(bsy*1.3), cpx, cpy);
+mTxt[1]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+
+
+bsy = bsyo; spy = spyo;
+
+
+cpx = (wszx-3.2*bsx-2*spx)/2;
+cpy += 2*spy+bsy;
+mButcf[2] = new PIButton(this, "SetPArm", 3100, bsx*1.2, bsyo, cpx, cpy);
+mButcf[2]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+cpx += bsx*1.2+spx;
+
+/******
+mButcf[3] = new PIButton(this, "Calc", 3200, bsx*0.8, bsyo, cpx, cpy);
+mButcf[3]->ActivatePress(true);
+mButcf[3]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+*****/ 
+cpx += bsx*0.8+spx; 
+
+
+mButcf[3] = new PIButton(this, "Dismiss", 3777, bsx*1.2, bsyo, cpx, cpy);
+mButcf[3]->SetBinding(PIBK_elastic, PIBK_elastic, PIBK_elastic, PIBK_elastic);
+
+cwdrwt = this;
+
+}
+
+/* --Methode-- */
+PICnTools::~PICnTools()
+{
+
+if (cwdrwt == this)  cwdrwt = NULL; 
+
+mCurCDrw = NULL;
+int i;
+for(i=0; i<5; i++)  
+  delete mLab[i];
+for(i=0; i<2; i++)  
+  delete mCkb[i];
+
+for(i=0; i<5; i++) delete mButcf[i]; 
+for(i=0; i<2; i++) delete mButdr[i]; 
+for(i=0; i<10; i++)  delete mOpt[i];
+for(i=0; i<2; i++)  delete mCasc[i];
+for(i=0; i<2; i++)  delete mTxt[i];
+
+delete mNlb;
+
+delete cmapv;
+}
+
+//  Declaration des differents attributs graphiques 
+#define MYNMXCOLORS 21
+static PIColors cols[MYNMXCOLORS] = { PI_NotDefColor, 
+                     PI_Black, PI_White, PI_Grey, 
+                     PI_Red, PI_Blue, PI_Green, 
+                     PI_Yellow, PI_Magenta, 
+                     PI_Cyan , PI_Turquoise, PI_NavyBlue, 
+                     PI_Orange, PI_SiennaRed, PI_Purple, 
+                     PI_LimeGreen, PI_Gold, PI_Violet, PI_VioletRed,
+                     PI_BlueViolet, PI_DarkViolet };
+
+static PIMarker mrk[12] = { PI_NotDefMarker, PI_DotMarker, 
+                            PI_PlusMarker, PI_CrossMarker,
+                            PI_CircleMarker, PI_FCircleMarker,
+                            PI_BoxMarker, PI_FBoxMarker,
+                            PI_TriangleMarker, PI_FTriangleMarker,
+                            PI_StarMarker, PI_FStarMarker };
+static PILineAtt line[10] = { PI_NotDefLineAtt, PI_ThinLine, PI_NormalLine, PI_ThickLine, 
+                              PI_ThinDashedLine, PI_DashedLine, PI_ThickDashedLine, 
+                              PI_ThinDottedLine, PI_DottedLine, PI_ThickDottedLine };
+
+static PIFontName fntname[5] = { PI_DefaultFont, PI_CourierFont, PI_HelveticaFont, 
+				 PI_TimesFont, PI_SymbolFont};
+static PIFontAtt fntatt[5] = { PI_NotDefFontAtt, PI_RomanFont, PI_BoldFont, 
+			       PI_ItalicFont, PI_BoldItalicFont };
+static PIFontSize fntsz[6] = { PI_NotDefFontSize, PI_TinySizeFont, PI_SmallSizeFont, 
+			       PI_NormalSizeFont, PI_BigSizeFont, PI_HugeSizeFont };
+
+/* --Methode-- */
+void PICnTools::Show()
+{
+PIWindow::Show();
+
+mNDr = 0;
+
+
+UpdateAttFromDrawer();
+return;
+}
+
+static t_contour_kind crbe_algo[3]={CONTOUR_KIND_LINEAR,CONTOUR_KIND_CUBIC_SPL,CONTOUR_KIND_BSPLINE};
+static t_contour_levels_kind lev_kind[4]={LEVELS_AUTO,LEVELS_NUM,LEVELS_INCREMENTAL,LEVELS_DISCRETE};
+
+/* --Methode-- */
+void PICnTools::UpdateAttFromDrawer()
+{
+
+PIContourDrawer* dr = mCurCDrw;
+
+mOpt[0]->SetValue(100);
+mOpt[1]->SetValue(200);
+mOpt[2]->SetValue(300);
+mOpt[3]->SetValue(405);
+mOpt[4]->SetValue(500);
+mOpt[5]->SetValue(600);
+mOpt[6]->SetValue(700);
+mOpt[7]->SetValue(800);
+mOpt[8]->SetValue(800);
+mOpt[9]->SetValue(800);
+
+// etats des check
+
+
+
+if (!dr)  return;
+
+mCkb[0]->SetState(mCurCDrw->IsLineOn());
+mCkb[1]->SetState(mCurCDrw->IsMarkOn());
+mCkb[2]->SetState(mCurCDrw->IsLabelOn());
+
+
+int ii,jj,kk;
+PIColors fc = dr->GetFgColAtt();
+for(kk=0; kk<MYNMXCOLORS; kk++) 
+  if (fc == cols[kk])  { mOpt[0]->SetValue(100+kk);  break; }    
+CMapId cmi = dr->GetColMapId();
+if (cmi == CMAP_OTHER)  mOpt[1]->SetValue(200);
+else for(kk=0; kk<PIColorMap::NumberStandardColorMaps(); kk++)
+  if (cmi == PIColorMap::GetStandardColorMapId(kk)) { mOpt[1]->SetValue(201+kk);  break; }  
+PIMarker mk = dr->GetMarker();   
+for(kk=0; kk<12; kk++) 
+  if (mk == mrk[kk])   { mOpt[2]->SetValue(300+kk);  break; }
+kk = dr->GetMarkerSz();
+if (kk < 1)  kk = 1;
+if (kk > 15) kk = 15;
+if (kk%2 == 0) kk++;
+mOpt[3]->SetValue(400+kk);
+PIFontName fnm = dr->GetFontName();
+PIFontAtt fat = dr->GetFontAtt();
+PIFontSize fsz = dr->GetFontSz();
+for(ii=0; ii<5; ii++)
+  if (fnm ==  fntname[ii]) break;
+for(kk=0; kk<5; kk++)
+  if (fat ==  fntatt[kk]) break;
+for(jj=0; jj<5; jj++) 
+  if (fsz == fntsz[jj])  break;
+
+if (ii >= 5) ii = 1;
+if (kk >= 5) kk = 1;
+if (jj >= 6) jj = 3;
+mOpt[4]->SetValue(500+ii);
+mOpt[6]->SetValue(700+kk);
+mOpt[7]->SetValue(800+jj);
+
+PILineAtt lat = dr->GetLineAtt();
+for(kk=0; kk<10; kk++)
+  if (lat == line[kk])   { mOpt[5]->SetValue(600+kk);  break; }
+
+// parametres pour calcul des niveaux 
+t_contour_kind kind = dr->GetCntKind()  ;
+for (ii=0 ; ii<3 ; ii++)
+  if(crbe_algo[ii]==kind)
+    {mOpt[8]->SetValue(900+ii); break;}
+
+t_contour_levels_kind lvkind = dr->GetCntLevelKind()  ;
+for (ii=0 ; ii<4 ; ii++)
+  if(crbe_algo[ii]==kind)
+    {mOpt[9]->SetValue(910+ii); break;}
+
+int nlv =  dr->NLevels();
+char txt[24];
+sprintf(txt,"%d",nlv);
+mTxt[0]->SetText(txt);  
+// cout << " nlv "<<nlv<< endl;
+if(dr->Levels()!=NULL){
+  string strg=" ";
+  for(int kk=0 ; kk<nlv ; kk++) {
+    cout << " dr->Level("<<kk<<") "<<endl;
+    cout <<dr->Level(kk)<<endl;
+    sprintf(txt,"%g",dr->Level(kk)); 
+    strg += txt;
+  }
+  mTxt[1]->SetText(strg);  
+  
+}
+
+
+return;
+}
+
+
+
+
+
+
+
+/* --Methode-- */
+//setting des parms du calcul des crbes de niveau
+void PICnTools::SetCalcParm(){
+
+PIContourDrawer* dr = mCurCDrw;
+if(dr==NULL)return;
+int k,ii,jj,kk;
+
+k = mOpt[8]->GetValue()-900; // type de ligne 
+dr->SetCntKind(crbe_algo[k]);
+
+k = mOpt[9]->GetValue()-910; // choix des niveaux des courbes 
+dr->SetCntLevelKind(lev_kind[k]);
+
+string str=mTxt[0]->GetText();
+ii=atoi(str.c_str());
+dr->SetNLevel(ii);
+char *buff;
+char *tmp;
+str=mTxt[1]->GetText();
+buff=strdup(str.c_str());
+double * zlev=NULL;
+//out << " k= "<<k<<endl;
+if(k==2) {
+  zlev = new double[2];
+  tmp = strtok(buff,",");
+  zlev[0] = atof(tmp);
+  tmp = strtok(NULL,",");
+  zlev[1]= atof(tmp);
+  dr->SetMyLevels(zlev,2);
+  // cas LEVELS_INCREMENTAL
+}else if (k==3){
+  // cas LEVELS_DISCRETE
+  zlev = new double[ii];
+  tmp = strtok(buff,",");
+  for(jj=0 ; jj<ii ; jj++){
+      zlev[jj]= atof(tmp);
+      cout << " jj "<<zlev[jj]<<endl;
+      tmp = strtok(NULL,",");
+  }
+  dr->SetMyLevels(zlev,ii);
+   
+
+
+} 
+
+//
+//
+
+dr->CalcContour();
+dr->Refresh();
+
+if(zlev!=NULL) {
+  delete[] zlev;
+  zlev=NULL;
+}
+
+}
+
+/* --Methode-- */
+void PICnTools::SetGraphicAtt(bool refr)
+{
+
+  if ( mCurCDrw== NULL )  return;
+  // etats des boutons 
+  
+  mCurCDrw->SetLineOn(mCkb[0]->GetState());
+  mCurCDrw->SetMarkOn(mCkb[1]->GetState());
+  mCurCDrw->SetLabelOn(mCkb[2]->GetState());
+
+  
+  int k,ii,jj,kk;
+  k = mOpt[0]->GetValue()-100; // -> COULEUR 
+  if ( (k < 0) || (k >= MYNMXCOLORS) )  k = 0;
+  mCurCDrw->SetColAtt(cols[k]);
+  double zmin=0;
+  double zmax=1.;
+  if(mCurCDrw!=NULL){
+    zmin = mCurCDrw->Zmin();
+    zmax = mCurCDrw->Zmax();
+  }
+  
+  k = mOpt[1]->GetValue()-200; // --> COLOR MAP 
+  if (k == 0){
+    mCurCDrw->SetColMapId(CMAP_OTHER);
+    
+  }else{
+    mCurCDrw->SetColMapId(PIColorMap::GetStandardColorMapId(k-1));
+    cmapv->SetColMapId(PIColorMap::GetStandardColorMapId(k-1),false, zmin , zmax,true);
+  }
+  k =  mOpt[2]->GetValue()-300; // Marker 
+  if ( (k < 0) || (k > 11) )  k = 0;
+  kk =  mOpt[3]->GetValue()-400; // marker SIZE
+  if ( (kk < 1) || (kk > 27) )  kk = 1;
+  mCurCDrw->SetMarkerAtt(kk, mrk[k]);
+
+  ii =  (mOpt[4]->GetValue()-500);
+  jj =  (mOpt[6]->GetValue()-700);
+  kk =  (mOpt[7]->GetValue()-800);
+  mCurCDrw->SetFont(fntname[ii], fntsz[kk], fntatt[jj]);
+
+  k =  mOpt[5]->GetValue()-600;
+  if ( (k < 0) || (k > 9) )  k = 0;
+  mCurCDrw->SetLineAtt(line[k]);
+
+  if (refr)
+    if (mCurBW !=NULL) 
+      mCurBW->Refresh();
+    else  
+      mCurCDrw->Refresh();
+  else  
+    mCurCDrw->Refresh();
+
+}
+
+/* --Methode-- */
+void PICnTools::Process(PIMessage msg, PIMsgHandler* /*sender*/, void* /*data*/)
+{
+  char strg[64];
+  PIMessage mod = ModMsg(msg);
+  msg = UserMsg(msg);
+  double zmin=0;
+  double zmax=1.;
+  if(mCurCDrw!=NULL){
+    zmin = mCurCDrw->Zmin();
+    zmax = mCurCDrw->Zmax();
+  }
+  
+  if(msg>=200&&msg<300){
+    int k = mOpt[1]->GetValue()-200; // _> COLOR MAP 
+    if (k == 0){
+      cmapv->SetColMapId(CMAP_OTHER,false, zmin , zmax,true);
+    }else{
+      cmapv->SetColMapId(PIColorMap::GetStandardColorMapId(k-1),false , zmin , zmax,true);
+    }
+    
+  }
+  switch (msg) {
+    case 2001 :
+      if(mCkb[0]->GetState()==false) 
+	mOpt[5]->SetUnSensitive();
+      else
+	mOpt[5]->SetSensitive();
+
+      break;
+
+    case 2002 :
+      if(mCkb[1]->GetState()==false){ 	
+	mOpt[2]->SetUnSensitive();
+	mOpt[3]->SetUnSensitive();
+      }else{	
+	mOpt[2]->SetSensitive();
+	mOpt[3]->SetSensitive();
+      }
+
+    case 2003 :
+      if(mCkb[2]->GetState()==false){ 
+	mOpt[4]->SetUnSensitive();
+	mOpt[6]->SetUnSensitive();
+	mOpt[7]->SetUnSensitive();
+      }else{
+	mOpt[4]->SetSensitive();
+	mOpt[6]->SetSensitive();
+	mOpt[7]->SetSensitive();
+      }
+      
+      
+      break;
+    case 2800 :
+      UpdateAttFromDrawer();
+      break;
+
+    case 2900 :
+      SetGraphicAtt(false);
+      break;
+
+    case 2950 :
+      SetGraphicAtt(true);
+      break;
+
+    case 3100 :
+      SetCalcParm();
+      break;
+
+
+    case 3200 :
+ 
+      break;
+
+    case 3777 :
+      Hide();
+      break;
+
+      /*   On ne change pas les options choisies   Reza+cmv 2/10/98  
+	   mOpt[0]->SetValue(100);
+	   mOpt[1]->SetValue(200);
+	   mOpt[2]->SetValue(300);
+	   mOpt[3]->SetValue(405);
+	   mOpt[4]->SetValue(500);
+	   mOpt[5]->SetValue(600);
+    */
+      break;
+
+    default :
+      break;
+  }
+
+return;
+}
+
+
+
Index: /tags/PIGC_V1/SophyaPI/PIGcont/picntools.h
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/picntools.h	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/picntools.h	(revision 1830)
@@ -0,0 +1,62 @@
+// This may look like C code, but it is really -*- C++ -*-
+// Outils de manipulation des PIContDrawer  O. Perdereau 11/2001
+// Relatif a un SEUL contourdrawer 
+// LAL (Orsay) / IN2P3-CNRS  DAPNIA/SPP (Saclay) / CEA
+
+#ifndef PICNTOOLS_H_SEEN
+#define PICNTOOLS_H_SEEN
+
+#include "pisysdep.h"
+#include PISTDWDG_H
+#include PIOPTMENU_H
+#include PIWIN_H
+#include PIAPP_H
+#include "pibwdggen.h"
+#include "pigncont.h"
+
+// --------- fenetre de dialogue pour modifs  options de trace des contours  ----------
+// --------------------------------------------------------------------------------
+
+
+class PICnTools : public PIWindow 
+{
+public :
+  PICnTools(PIApplication* par);
+  ~PICnTools();
+  virtual void  Show();
+
+  static PICnTools* curlutw;
+  static void ShowPICnTools();
+  static void HidePICnTools();
+  static void SetCurrentBaseWdg(PIBaseWdgGen* cbw);
+  static void SetCurrentCnDrw(PIContourDrawer* drw);
+  static PIBaseWdgGen* GetCurrentBaseWdg();
+  static PIContourDrawer* GetCurrentCnDrw();
+  inline PICMapView* CMapVW() { return cmapv; }
+protected:
+  virtual void  UpdateAttFromDrawer();
+  virtual void  SetGraphicAtt(bool );
+  virtual void  SetCalcParm();
+  virtual void  Process(PIMessage msg, PIMsgHandler* sender, void* data=NULL);
+
+private:
+// Pour couleur, fontes, etc 
+
+  PICMapView* cmapv;
+  PICheckBox*  mCkb[3];
+  PILabel* mLab[5];
+  PIOptMenu* mOpt[10];
+  PIMenu* mCasc[5];
+  PILabel* mNlb;
+  PIText* mTxt[2];
+  PIButton* mButdr[2];
+  PIButton* mButcf[5];
+  int mNDr;
+  int mtyp, msz;
+//
+  static PIBaseWdgGen* mCurBW;
+  static PIContourDrawer* mCurCDrw;
+  static PICnTools* cwdrwt;
+};
+
+#endif
Index: /tags/PIGC_V1/SophyaPI/PIGcont/pigncont.cc
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/pigncont.cc	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/pigncont.cc	(revision 1830)
@@ -0,0 +1,664 @@
+#include "machdefs.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream.h>
+#include <math.h>
+
+#include <typeinfo>
+
+#include <vector>
+#include <string>
+
+#include "ctimer.h"
+#include "piacmd.h"
+#include "nobjmgr.h"
+#include "pistdimgapp.h"
+
+#include "pidrawer.h"
+#include "nomgadapter.h"
+
+#include "histos.h"
+#include "histos2.h"
+#include "hisprof.h"
+#include "ntuple.h"
+
+#include "servnobjm.h"
+#include <nbtri.h> 
+
+#include "picntools.h"
+#include "pigncont.h"
+//A-Virer #include "pidrwtools.h"
+//A-Virer #include "piaxestools.h"
+
+extern "C" { 
+#include "contour.h"
+#include "alloc.h"
+#include "set.h"
+}
+
+TBOOLEAN multiplot;
+double zero =0;
+TBOOLEAN polar;
+int xleft, xright, ybot, ytop;
+float xsize; 
+float ysize; 
+float xoffset, yoffset;
+
+
+GNUPlotContour::GNUPlotContour(struct iso_curve * iso){
+  int k=0;
+  iso1 = iso;
+  struct iso_curve *iso_cur = iso ;
+  struct iso_curve *iso_old;
+  _xmin = 1.e37;
+  _ymin = 1.e37;
+  _xmax = -1.e37;
+  _ymax = -1.e37;
+  _myiso = NULL;
+  _zmin = 1.e37;
+  _zmax = -1.e37;
+  while(iso_cur){
+    k++;
+    iso_cur = iso_cur->next;
+    _ny = iso_cur->p_count;
+    for(int k=0 ; k<_ny ; k++){
+      double xx=iso->points[k].x;
+      double yy=iso->points[k].y;
+      double zz=iso->points[k].z;
+      if(xx<_xmin) _xmin = xx;
+      if(xx>_xmax) _xmax = xx;
+      if(yy<_ymin) _ymin = yy;
+      if(yy>_ymax) _ymax = yy;
+      if(zz<_zmin) _zmin = zz;
+      if(zz>_zmax) _zmax = zz;
+
+    }
+  }
+
+  _nx = k;
+  _contours = NULL;
+  My_Levels = NULL;
+  
+  _npolys=-1;  
+  _Dxp =_Dyp =1.;
+  _Exy =_Invy =_Invy = FALSE ;
+
+}
+
+
+GNUPlotContour::GNUPlotContour(NTupleInterface *nt,int *nz=NULL,double *z=NULL,double *dz=NULL){
+  _xmin = 1.e37;
+  _ymin = 1.e37;
+  _xmax = -1.e37;
+  _ymax = -1.e37;
+  _zmin = 1.e37;
+  _zmax = -1.e37;
+  char *nom[4]={"x","y","z","niso"};
+  _myiso = new NTuple(4,nom);  
+  //_myiso->Show();
+  float xnt[4];
+  double mn,mx;
+
+  nt->GetMinMax(0,_xmin,_xmax);
+  nt->GetMinMax(1,_ymin,_ymax);
+  nt->GetMinMax(2,_zmin,_zmax);
+
+  
+  struct iso_curve *istmp = NULL ;
+  struct iso_curve *iso = NULL ;
+  struct iso_curve *iso_old = NULL ;
+  double x,y;
+  int nen = nt->NbLines();
+  
+  // histo des y
+  for(int k= 0 ; k<1000 ; k++){
+    istmp = iso_alloc(100);
+    iso_free(istmp);
+  }
+
+#define NUMPTMIN 25
+#define MAXISO 100
+#define MINISO 5
+
+  Histo prep(_ymin,_ymax,nen);  
+  for(int ix=0 ; ix<nen ; ix++) prep.Add(nt->GetCell(ix,1),1.);
+  Histo Intg = prep;
+  Intg.HInteg(0);
+  //prep.Print();
+  //Intg.Print();
+  float *bin;
+  int *cnt;
+
+  int niso = sqrt( (double)nen);
+  bin = new float[niso+1];
+  cnt = new int[niso+1];
+
+  for(int i=0 ; i<niso ; i++){    
+    bin[i] = prep.BinHighEdge(prep.BinPercent( ((double)i)/(double)niso));
+    cnt[i] = Intg(prep.BinPercent( ((double)i)/(double)niso)); 
+  }
+
+  bin[niso]=_ymax;
+  
+  cnt[niso] = prep.NEntries();
+  int numx = 0;
+  int numi = 9999;
+  for(int i=0;i<niso ; i++){
+    int numy = cnt[i+1]-cnt[i];
+    
+    if(numy>0){
+      numx = numx>numy ? numx : numy ;
+      numi = numy>numi ? numi : numy ;
+    }
+  }
+
+  
+  if (numi<=0) {
+    
+    return;
+  }
+  _nx = 0;
+  
+  for(int i=0;i<niso ; i++){
+    double ymi = bin[i];
+    double ymx = bin[i+1];
+    int numy = cnt[i+1]-cnt[i];
+    
+    if(numy<=0) continue;
+    double *xp;
+    xp = new double[numy];
+    
+    int *indx;
+    indx = new int[numy];
+
+    istmp = iso_alloc(numy) ; // structure temporaire 
+    
+    int ip=0;
+    
+    for(int ix=0 ; ix<nen ; ix++){
+      
+      if(nt->GetCell(ix,1)>=ymi&&nt->GetCell(ix,1)<ymx){
+        
+
+	istmp->points[ip].type = INRANGE;
+	istmp->points[ip].x = nt->GetCell(ix,0);
+	istmp->points[ip].y = nt->GetCell(ix,1);
+	istmp->points[ip].z = nt->GetCell(ix,2);
+	indx[ip] = ip;
+	xp[ip++] = nt->GetCell(ix,0);
+      }      
+    }
+    
+    // tri des points de l'iso_curve 
+    //
+    istmp->p_count = ip;
+    
+    tri_double(xp,indx,ip);
+    
+    // remplissage structure finale 
+    
+    iso = iso_alloc(numx) ;
+    if( _nx == 0 ) iso1 = iso; //premiere structure remplie 
+    if (iso_old) iso_old->next = iso;
+    for(int ix=0 ; ix<numy   ; ix++){
+      int jx = indx[ix];
+    
+      
+      iso->points[ix].type = INRANGE ;
+      iso->points[ix].x = istmp->points[jx].x ;
+      iso->points[ix].y = istmp->points[jx].y ;
+      iso->points[ix].z = istmp->points[jx].z ;
+      xnt[0] = iso->points[ix].x;
+      xnt[1] = iso->points[ix].y;
+      xnt[2] = iso->points[ix].z;
+      xnt[3] = _nx;
+      //cout << " fill avec "<<xnt[0]<<"  "<<xnt[1] <<" "<<xnt[2] << " "<<xnt[3]<<endl;
+      // _myiso->Show();
+ 
+      _myiso->Fill(xnt);
+    }
+    //cout << " completion "<<ip<<" =>> "<<numx<<endl;
+    // on complete 
+    for( int ix=ip; ix<numx ; ix++ ){
+      int jmx = indx[ip-1];
+      iso->points[ix].type = INRANGE ;
+      iso->points[ix].x = istmp->points[jmx].x ;
+      iso->points[ix].y = istmp->points[jmx].y ;
+      iso->points[ix].z = istmp->points[jmx].z ;
+      xnt[0] = iso->points[ix].x;
+      xnt[1] = iso->points[ix].y;
+      xnt[2] = iso->points[ix].z;
+      xnt[3] = _nx;
+      //cout << " fill avec "<<xnt[0]<<"  "<<xnt[1] <<" "<<xnt[2] << " "<<xnt[3]<<endl;
+      _myiso->Fill(xnt);
+      //_myiso->Show();
+    }
+    
+    iso->p_count = numx;
+
+    iso_old = iso;
+    _nx++;
+    //cout << " destruction structure tempo "<<_nx <<endl;
+    if(istmp){
+      iso_free( istmp );// destruction de la structure temporaire 
+      istmp = NULL;
+    }
+    if(xp){delete[] xp ; xp=NULL;}
+    if(indx) {delete[] indx ; indx = NULL;}
+
+  }
+  //if(_myiso!=NULL)_myiso->Write("myiso.ppf"); 
+  //cout << " fin du remplissage des iso_curves xmin,max "<< _xmin<<","<<_xmax<<" ymin,max "<<_ymin<<","<<_ymax<<endl;
+  _contours = NULL;
+  My_Levels = NULL;
+  _npolys=-1;  
+  Contour_Levels_kind = LEVELS_AUTO ;
+  Contour_kind = CONTOUR_KIND_LINEAR ;
+  set_contour_kind(Contour_kind);
+  Contour_Levels = 5;
+  set_contour_levels(Contour_Levels);
+  set_contour_levels_kind(Contour_Levels_kind);
+  
+ 
+  if(dz!=NULL){
+    My_Levels = new double[2];
+    Contour_Levels = *nz;
+    My_Levels[0] = *z ;
+    My_Levels[1] = *dz;
+    Contour_Levels_kind = LEVELS_INCREMENTAL ;
+    set_contour_levels(Contour_Levels);
+    set_contour_levels_kind(Contour_Levels_kind);
+    set_contour_levels_list(My_Levels);//,2);
+
+  }else  if(nz!=NULL){
+    
+    Contour_Levels = *nz;
+    Contour_Levels_kind = LEVELS_NUM ;
+    set_contour_levels(Contour_Levels);
+    set_contour_levels_kind(Contour_Levels_kind);
+  }
+  
+}
+
+GNUPlotContour::GNUPlotContour(P2DArrayAdapter*arr,int *nz=NULL,double *z=NULL,double *dz=NULL){
+  _xmin = 1.e37;
+  _ymin = 1.e37;
+  _xmax = -1.e37;
+  _ymax = -1.e37;
+  _zmin = 1.e37;
+  _zmax = -1.e37;
+  // _nx = arr->XSize()-1;
+  // _ny = arr->YSize()-1;
+   _nx = arr->XSize()  ;
+   _ny = arr->YSize()  ;
+   //cout << " size x,y "<<_nx<<" , "<<_ny <<endl;
+   char *nom[4]={"x","y","z","niso"};
+   
+   _myiso = new NTuple(4,nom);  
+   float xnt[4];
+   struct iso_curve *iso =NULL ;
+   struct iso_curve *iso_old =NULL ;
+   double x,y;
+   // Calcul de la taille en X et en Y d'un pixel
+   double dxp, dyp;
+   double x1,y1;
+   arr->XYfromxy(0,0,x,y);
+   arr->XYfromxy(1,1,x1,y1);
+   _Dxp = (x1-x);
+   _Dyp = (y1-y);
+   double zz;
+   for(int ix=0 ; ix<_nx ; ix++){
+     iso = iso_alloc(_ny) ;
+     if(ix==0) iso1 = iso;
+     if (iso_old) iso_old->next = iso;
+     for(int iy = 0 ; iy<_ny ; iy++){
+       /*i = ix*100+iy;*/
+       arr->XYfromxy(ix,iy,x,y);
+
+       iso->points[iy].type = INRANGE;
+       iso->points[iy].x = x + _Dxp/2.;
+       iso->points[iy].y = y + _Dyp/2.;
+       zz=arr->Value(ix,iy);
+       iso->points[iy].z = zz;
+       //iso->points[iy].z = ((x-25)*(x-25) + (y-25)*(y-25))*.05 ;
+       if(x <_xmin) _xmin = x;
+       if(x >_xmax) _xmax = x;
+       if(y <_ymin) _ymin = y;
+       if(y >_ymax) _ymax = y;
+       if(zz <_zmin) _zmin = zz;
+       if(zz >_zmax) _zmax = zz;
+       //printf("ix=%d iy=%d x=%f y=%f z=%f \n",ix,iy,iso->points[iy].x,iso->points[iy].y,iso->points[iy].z);
+       xnt[0] = iso->points[iy].x;
+       xnt[1] = iso->points[iy].y;
+       xnt[2] = iso->points[iy].z;
+       xnt[3] = _nx;
+
+       _myiso->Fill(xnt);
+     }
+     iso->p_count = _ny;
+     iso_old = iso;
+     /*printf(" remplissage %lx old %lx old next %lx \n",iso,iso_old,iso->next,iso_old->next);*/
+     /*iso->next = iso ;*/
+   }
+   //cout << " fin du remplissage des iso_curves xmin,max "<< _xmin<<","<<_xmax<<" ymin,max "<<_ymin<<","<<_ymax<<endl;
+   _contours = NULL;
+   My_Levels = NULL;
+  _npolys=-1;  
+  Contour_kind = CONTOUR_KIND_LINEAR ;
+  set_contour_kind(Contour_kind);
+  Contour_Levels_kind = LEVELS_AUTO ;
+  Contour_Levels = 5;
+  set_contour_levels(Contour_Levels);
+  set_contour_levels_kind(Contour_Levels_kind);
+  
+
+  if(dz!=NULL){
+    My_Levels = new double[2];
+    Contour_Levels = *nz;
+    My_Levels[0] = *z ;
+    My_Levels[1] = *dz;
+    Contour_Levels_kind = LEVELS_INCREMENTAL ;
+    set_contour_levels(Contour_Levels);
+    set_contour_levels_kind(Contour_Levels_kind);
+    set_contour_levels_list(My_Levels);//,2);
+
+  }else  if(nz!=NULL){
+    Contour_Levels = *nz;
+    Contour_Levels_kind = LEVELS_NUM ;
+    set_contour_levels(Contour_Levels);
+    set_contour_levels_kind(Contour_Levels_kind);
+  }
+  
+}
+
+GNUPlotContour::~GNUPlotContour(){
+ cout << " destructeur de GNUPlotContour "<<endl;
+  
+ struct iso_curve *iso_cur;
+ struct iso_curve *iso_nx;
+ iso_cur = iso1;
+ 
+  while(iso_cur){
+    // cout << " ~GNUPlotContour() : destruction de  iso1 "<< iso_cur << endl;
+   iso_nx = iso_cur->next;  
+   iso_free(iso_cur);iso_cur = NULL;
+   
+   iso_cur = iso_nx;
+   
+   }
+ 
+ iso1 = NULL;
+ 
+ struct gnuplot_contours *cntcur = _contours;
+ struct gnuplot_contours *cntold;
+
+ while(cntcur) {
+   //cout << " ~GNUPlotContour() : destruction de _contours" << _contours <<endl; 
+   cntold = cntcur;
+   cntcur = cntold->next;
+   gp_free(cntold);
+   cntold=NULL;
+ } 
+ 
+ _contours = NULL;
+
+ if(My_Levels) {
+   //cout << " ~GNUPlotContour() : destruction de MyLevels "<<My_Levels <<endl;
+   delete[] My_Levels; My_Levels=NULL;}
+
+ if(_myiso!=NULL){ delete _myiso; _myiso=NULL;}
+
+} 
+
+void GNUPlotContour::CalcContour(){
+ cout << " GNUPlotContour::CalcContour(): determination des contours "<<endl;
+
+ set_contour_kind(Contour_kind);
+ if(Contour_Levels_kind == LEVELS_INCREMENTAL){
+    
+    set_contour_levels(Contour_Levels);
+    set_contour_levels_kind(Contour_Levels_kind);
+    //free_contour_levels_list();
+    set_contour_levels_list(My_Levels);//,2);
+
+  }else  if(Contour_Levels_kind == LEVELS_DISCRETE ){
+    set_contour_levels(Contour_Levels);
+    set_contour_levels_kind(Contour_Levels_kind);
+    //free_contour_levels_list();
+    set_contour_levels_list(My_Levels);//,Contour_Levels);
+    
+  }else  if(Contour_Levels_kind == LEVELS_NUM ){
+    
+    set_contour_levels(Contour_Levels);
+    set_contour_levels_kind(Contour_Levels_kind);
+  }else if(Contour_Levels_kind == LEVELS_AUTO){
+    Contour_Levels = 5;
+    set_contour_levels(Contour_Levels);
+    set_contour_levels_kind(Contour_Levels_kind);
+  }
+
+ if(_contours) {
+   //cout << " GNUPlotContour::CalcContour(): destruction des contours "<<endl;
+   struct gnuplot_contours *cntcur = _contours;
+   struct gnuplot_contours *cntold;
+   while(cntcur) {
+     cntold = cntcur;
+     cntcur = cntold->next;
+     gp_free(cntold);
+   }
+   _contours = NULL;
+   
+ }
+
+   
+ 
+ //struct gnuplot_contours *cntcur;
+ //struct gnuplot_contours *cntold;
+ 
+ _contours = contour (_nx,iso1);
+ //free_contour_levels_list();
+
+ 
+}
+//_+_ Methode 
+void 
+GNUPlotContour::SetMyLevels(double *ptr,int k){
+  
+  if(My_Levels!=NULL) delete My_Levels;
+  My_Levels = new double[k];
+  for(int i=0 ; i<k ; i++)My_Levels[i] = ptr[i];
+}
+
+// Class PIContourDrawer
+/* --Methode-- */
+PIContourDrawer::PIContourDrawer(P2DArrayAdapter* arr,bool autodel,int *nz,double *z0,double *dz )
+  :GNUPlotContour(arr,nz,z0,dz)
+
+{
+  _arr = arr;
+  _nti = NULL;
+  _nz=5;
+  if(nz!=NULL)_nz = *nz;
+  _autodel = autodel;
+  InitAtts();
+}
+
+PIContourDrawer::PIContourDrawer(NTupleInterface* nti,bool autodel,int *nz,double *z0,double *dz )
+  :GNUPlotContour(nti,nz,z0,dz)
+
+{
+  _arr = NULL;
+  _nti = nti;
+  _nz=5;
+  if(nz!=NULL)_nz = *nz;
+  _autodel = autodel;
+  InitAtts();
+}
+
+
+
+/* --Methode-- */
+PIContourDrawer::~PIContourDrawer()
+{
+  if(_autodel){ 
+    if(_arr!=NULL){delete _arr; _arr=NULL;}
+    if(_nti!=NULL){delete _nti ; _nti=NULL;}
+  }
+  //if (mColorMap!=NULL){delete  mColorMap;  mColorMap=NULL;}
+}
+
+/* --Methode-- */
+void PIContourDrawer::Draw(PIGraphicUC* g, double xmin, double ymin, double xmax, double ymax)
+{
+  // double xmin, double ymin, double xmax, double ymax LIMITES DE LA ZONE A REFRESHER
+  //PIGrCoord * xx = new PIGrCoord[10];
+  //PIGrCoord * yy = new PIGrCoord[10];
+  PIGrCoord * xx =NULL;
+  PIGrCoord * yy =NULL;
+  PIGrCoord xa;
+  PIGrCoord ya;
+  char buff[64];
+  // On met les bons attributs graphiques
+  //SelGraAtt(g);
+  g->SaveGraphicAtt();
+  struct gnuplot_contours *cntcur;
+  struct gnuplot_contours *cntold;
+  cntcur = _contours;
+
+  double a,b;
+  int tot=0;
+  bool Invx, Invy, Exy;
+  PIColorMap * cmap = NULL;
+  if(mCmapid !=CMAP_OTHER ){
+    cmap = new PIColorMap(mCmapid);
+  }
+  int numdr = 0;
+  while(cntcur){
+    int npts = cntcur->num_pts ;
+    if(npts<0) continue;
+
+    xx = new PIGrCoord[npts];
+    yy = new PIGrCoord[npts];
+    double lev = (cntcur->coords)[0].z ;
+    for(int l=0 ; l< npts ; l++){
+      a = (cntcur->coords)[l].x;
+      b = (cntcur->coords)[l].y;
+      
+      xx[l] = a ;
+      yy[l] = b ;
+
+    }
+
+    // choix de la couleur 
+    
+    g->SelForeground(mFCol);
+    if(mCmapid !=CMAP_OTHER ){
+      if(_zmax-_zmin!=0){
+	int kc = ((lev - _zmin)/(_zmax-_zmin))*cmap->NCol();	
+	
+	g->SelForeground(*cmap,kc);
+      }
+    }    
+    // traitement des  des markers 
+    if(IsMarkOn()==true){
+      
+      g->SelMarker(mMSz,mMrk);
+      g->DrawMarkers(xx,yy,npts);
+    }
+    // traitement des lignes 
+    if(mLineOn==true){
+      
+      g->SelLine(mLAtt);
+      g->DrawPolygon(xx,yy,npts,false);
+    }
+    // traitement des labels     
+    // SIMPLISTE POUR L'INSTANT : UN AFFICHAGE 
+    if(mLabelOn==true){
+      // 
+      
+      char strg[10];
+      sprintf(strg,"%g",lev);
+      PIFont myfont(mFName);
+      myfont.SetFontAtt(mFAtt);
+      myfont.SetFontSz(mFSz);
+      
+      g->SelFont(myfont);
+      double px,py;
+      px = (cntcur->coords)[0].x+2*(Xmax()-Xmin())/100.;
+      py = (cntcur->coords)[0].y;
+      g->DrawString(px,py,strg, PI_HorizontalCenter&&PI_VerticalCenter );
+      
+    }
+    numdr++;
+    delete[] xx;
+    xx=NULL;
+    delete[] yy;
+    yy=NULL;
+    cntcur = cntcur->next;
+    tot++;
+  }    
+  
+  //cout << "PIContourDrawer::Draw  fin de trace des "<<tot<< " polygones "<<endl; 
+ 
+  // 
+  if(cmap !=NULL)delete cmap;
+  
+  g->RestoreGraphicAtt();
+}
+
+/* --Methode-- */
+void PIContourDrawer::UpdateLimits()
+{
+  // Doit calculer les limites
+
+  double xmn = _xmin;
+  double xmx = _xmax;  
+  double ymn = _ymin;
+  double ymx = _ymax;
+  
+  SetLimits(xmn, xmx, ymn, ymx);
+}
+
+
+/* --Methode-- */
+void PIContourDrawer::ShowControlWindow(PIBaseWdgGen* wdg)
+{
+  PICnTools::SetCurrentBaseWdg(wdg);
+  PICnTools::SetCurrentCnDrw(this);
+  PICnTools::ShowPICnTools();
+}
+
+
+/* --Methode-- */
+bool PIContourDrawer::IsLabelOn(){
+  return (mLabelOn);
+}
+
+bool PIContourDrawer::IsLineOn(){
+  return (mLineOn);
+}
+
+bool PIContourDrawer::IsMarkOn(){
+  return (mMarkOn);
+}
+
+void PIContourDrawer::SetLabelOn(bool state){
+  mLabelOn = state;
+}
+
+void PIContourDrawer::SetMarkOn(bool state){
+  mMarkOn = state;
+}
+
+void PIContourDrawer::SetLineOn(bool state){
+  mLineOn = state;
+}
+
+void PIContourDrawer::InitAtts(){
+  mLineOn = false;
+  mMarkOn = true;
+  mLabelOn = false;
+}
+
+
+
Index: /tags/PIGC_V1/SophyaPI/PIGcont/pigncont.h
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/pigncont.h	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/pigncont.h	(revision 1830)
@@ -0,0 +1,106 @@
+// This may look like C code, but it is really -*- C++ -*-
+// Classe  GNUPlotContour   : interface avec le code de gnuplot pour tracer les contours
+// Classe  PIContourDrawer  : trace des contours ds PI 
+// (c) DAPNIA (CEA)                 LAL (IN2P3/CNRS)
+//                           R. Ansari O.Perdereau  11/2001
+//=========================================================//
+
+#ifndef PIGNCONT_H_SEEN 
+#define PIGNCONT_H_SEEN 
+#include "pistdimgapp.h"
+#include "piscdrawwdg.h"
+
+#include "ntuple.h"
+extern "C" { 
+#include "contour.h"
+#include "alloc.h"
+#include "set.h"
+}
+
+class GNUPlotContour{
+public:
+  GNUPlotContour(P2DArrayAdapter*,int *nz,double *z0,double *dz);
+  GNUPlotContour(NTupleInterface *nt,int *nz,double *z0,double *dz);
+  GNUPlotContour(struct iso_curve *);
+  ~GNUPlotContour();
+  void CalcContour();
+  double Level(int i){return My_Levels[i];};
+  const double *Levels(){return My_Levels;};
+  void SetMyLevels(double *ptr,int k);
+  int NLevels(){return Contour_Levels;};
+  void SetNLevel(int n){Contour_Levels=n;};
+  int NContours(){return _npolys;};
+  void SetCntLevelKind(t_contour_levels_kind kind){Contour_Levels_kind = kind;};
+  t_contour_levels_kind GetCntLevelKind(){return Contour_Levels_kind ;};
+  void SetCntKind(t_contour_kind kind){Contour_kind = kind;};
+  t_contour_kind GetCntKind(){return Contour_kind ;};
+
+  NTuple * & MyIso(){return  _myiso;};
+
+  double Xmin(){return _xmin;};
+  double Xmax(){return _xmax;};
+  double Ymin(){return _ymin;};
+  double Ymax(){return _ymax;};
+  double Zmin(){return _zmin;};
+  double Zmax(){return _zmax;};
+
+private: 
+
+protected: 
+  NTuple * _myiso;
+  int _nx,_ny;
+  struct iso_curve *iso1 ;
+  bool _Invx,_Invy,_Exy;
+  double _Dxp, _Dyp;
+  struct gnuplot_contours *_contours;
+  int Contour_Levels;
+  
+  double *My_Levels; 
+  
+  double _xmin,_xmax,_ymin,_ymax,_zmin,_zmax;
+  int _npolys;
+  t_contour_levels_kind Contour_Levels_kind;
+  t_contour_kind Contour_kind;
+  
+  
+};
+
+class PIContourDrawer : public PIDrawer, public GNUPlotContour{
+public:
+                    PIContourDrawer(P2DArrayAdapter* arr,bool autodel ,int *nz=NULL,double *z0=NULL,double *dz=NULL);
+                    PIContourDrawer(NTupleInterface* nti,bool autodel ,int *nz=NULL,double *z0=NULL,double *dz=NULL);
+  virtual           ~PIContourDrawer();
+  virtual void      Draw(PIGraphicUC* g, double xmin, double ymin, double xmax, double ymax);
+  virtual void      UpdateLimits();
+  virtual bool      IsLineOn();
+  virtual bool      IsMarkOn();  
+  virtual bool      IsLabelOn();
+  virtual void      SetLineOn(bool );
+  virtual void      SetMarkOn(bool);
+  virtual void      SetLabelOn(bool);
+
+  //Init  attributs graphiques 
+
+  void InitAtts();  
+  
+  
+//  Methode permettant l'affichage d'une fenetre de controle specialisee 
+  virtual void      ShowControlWindow(PIBaseWdgGen* wdg);
+
+
+protected:
+  bool _autodel;
+  int _nz;
+  P2DArrayAdapter* _arr;
+  NTupleInterface* _nti;
+  bool mLineOn,mMarkOn,mLabelOn;
+  
+
+  
+
+};
+
+
+//
+
+#endif
Index: /tags/PIGC_V1/SophyaPI/PIGcont/set.c
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/set.c	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/set.c	(revision 1830)
@@ -0,0 +1,55 @@
+#include "set.h"
+
+void set_contour_kind(t_contour_kind in){
+  contour_kind = in;
+}
+
+t_contour_kind get_contour_kind(){
+  return (contour_kind);
+}
+
+void set_contour_levels_kind(t_contour_levels_kind in){
+ contour_levels_kind = in;
+
+}
+
+t_contour_levels_kind get_contour_levels_kind(){
+
+  return(contour_levels_kind);
+}
+
+void set_contour_levels(int num){
+  contour_levels = num;
+}
+
+int get_contour_levels(){
+  return(contour_levels);
+}
+void set_contour_levels_list(double *vec){
+
+  contour_levels_list=vec;
+}
+
+/*******
+void set_contour_levels_list(double *vec,int sz){
+  int i;
+  printf(" <set_contour_levels_list> sz %d \n",sz);
+  contour_levels_list = (double *)  malloc(sz*sizeof(double));
+  for (i=0 ; i<sz ; i++){
+    printf(" <set_contour_levels_list> i %d vec %g\n",i,vec[i]);
+    contour_levels_list[i] = vec[i];
+    
+  }
+}
+void free_contour_levels_list(){
+  if(contour_levels_list != NULL){
+    free(contour_levels_list);
+    contour_levels_list = NULL;
+  }
+}
+*****/
+double * get_contour_levels_list(){
+
+return(contour_levels_list);
+}
+ 
Index: /tags/PIGC_V1/SophyaPI/PIGcont/set.h
===================================================================
--- /tags/PIGC_V1/SophyaPI/PIGcont/set.h	(revision 1830)
+++ /tags/PIGC_V1/SophyaPI/PIGcont/set.h	(revision 1830)
@@ -0,0 +1,23 @@
+#ifndef SET_H
+#define SET_H
+
+#include <stdio.h>
+#include <stdlib.h>  
+#include "contour.h"
+
+void set_contour_kind(t_contour_kind);
+t_contour_kind get_contour_kind();
+
+void set_contour_levels_kind(t_contour_levels_kind);
+t_contour_levels_kind get_contour_levels_kind();
+
+void set_contour_levels(int);
+int get_contour_levels();
+void set_contour_levels_list(double *);
+/***
+void set_contour_levels_list(double *,int);
+void free_contour_levels_list();
+***/
+double * get_contour_levels_list();
+
+#endif
