Logo Search packages:      
Sourcecode: v-sim version File versions  Download package

pot2surf.c

/*   EXTRAITS DE LA LICENCE
      Copyright CEA, contributeurs : Luc BILLARD, Damien
      CALISTE, Olivier D'Astier, laboratoire L_Sim, (2001-2005)
  
      Adresses mèl :
      BILLARD, non joignable par mèl ;
      CALISTE, damien P caliste AT cea P fr.
      D'ASTIER, dastier AT iie P cnam P fr.

      Ce logiciel est un programme informatique servant à visualiser des
      structures atomiques dans un rendu pseudo-3D. 

      Ce logiciel est régi par la licence CeCILL soumise au droit français et
      respectant les principes de diffusion des logiciels libres. Vous pouvez
      utiliser, modifier et/ou redistribuer ce programme sous les conditions
      de la licence CeCILL telle que diffusée par le CEA, le CNRS et l'INRIA 
      sur le site "http://www.cecill.info".

      Le fait que vous puissiez accéder à cet en-tête signifie que vous avez 
      pris connaissance de la licence CeCILL, et que vous en avez accepté les
      termes (cf. le fichier Documentation/licence.fr.txt fourni avec ce logiciel).
*/

/*   LICENCE SUM UP
      Copyright CEA, contributors : Luc BILLARD and Damien
      CALISTE and Olivier D'Astier, laboratoire L_Sim, (2001-2005)

      E-mail addresses :
      BILLARD, not reachable any more ;
      CALISTE, damien P caliste AT cea P fr.
      D'ASTIER, dastier AT iie P cnam P fr.

      This software is a computer program whose purpose is to visualize atomic
      configurations in 3D.

      This software is governed by the CeCILL  license under French law and
      abiding by the rules of distribution of free software.  You can  use, 
      modify and/ or redistribute the software under the terms of the CeCILL
      license as circulated by CEA, CNRS and INRIA at the following URL
      "http://www.cecill.info". 

      The fact that you are presently reading this means that you have had
      knowledge of the CeCILL license and that you accept its terms. You can
      find a copy of this licence shipped with this software at Documentation/licence.en.txt.
*/
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>

#include <math.h>

#include <visu_object.h>
#include <visu_tools.h>
#include "surfaces.h"
#include "pot2surf.h"

#define ISOSURFACES_FLAG_POTENTIAL "# potentialValue"


static int nx, ny, nz;
static double dxx, dyx, dyy, dzx, dzy, dzz;
static double exx, eyx, eyy, ezx, ezy, ezz;
static int nxm1, nym1, nzm1;
static int nz3, nyz3, nxyz3;
static double dxx1, dyx1, dyy1, dzx1, dzy1, dzz1;
static double ***f = NULL;
static double isovalue;
static FILE *in = NULL, *out = NULL;
static int n1, n2;
static int n1_tot, n2_tot;

static int *itab = NULL;
static double *xs = NULL, *ys = NULL, *zs = NULL, *xns = NULL, *yns = NULL, *zns = NULL;

/* Tables from :
   http://www.swin.edu.au/astronomy/pbourke/modelling/polygonise/index.html
   (valid link on February 25, 2002)
   
   YOU WILL FIND MANY VERY INTERESTING PAGES at Paul Bourke's site:
   http://astronomy.swin.edu.au/pbourke/
   
   Other routines by Luc Billard, on March 1, 2002
*/

static int edgeTable[256]={
0x0  , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0   };

static int triTable[256][16] =
{{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
{3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
{3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
{3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
{9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
{2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
{8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
{4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
{3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
{1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
{4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
{4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
{5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
{2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
{9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
{0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
{2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
{10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
{5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
{5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
{9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
{10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
{8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
{2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
{7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
{2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
{11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
{5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
{11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
{11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
{9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
{2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
{6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
{3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
{6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
{10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
{6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
{8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
{7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
{3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
{0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
{9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
{8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
{5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
{0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
{6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
{10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
{10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
{8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
{1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
{0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
{10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
{3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
{6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
{9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
{8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
{3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
{6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
{0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
{10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
{10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
{2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
{7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
{7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
{2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
{1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
{11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
{8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
{0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
{7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
{10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
{2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
{6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
{7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
{2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
{10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
{10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
{0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
{7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
{6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
{8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
{9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
{6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
{4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
{10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
{8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
{0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
{1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
{8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
{10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
{4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
{10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
{11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
{9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
{6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
{7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
{3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
{7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
{3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
{6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
{9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
{1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
{4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
{7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
{6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
{3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
{0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
{6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
{0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
{11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
{6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
{5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
{9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
{1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
{10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
{0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
{5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
{10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
{11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
{9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
{7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
{2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
{8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
{9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
{9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
{1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
{9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
{5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
{0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
{10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
{2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
{0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
{0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
{9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
{5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
{3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
{5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
{8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
{0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
{9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
{1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
{3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
{4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
{9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
{11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
{11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
{2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
{9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
{3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
{1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
{4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
{3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
{0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
{1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}};

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

static double gx(int i, int j, int k) {
   if(i > 0 && i < nxm1) {
      return (f[i+1][j][k]-f[i-1][j][k]) / 2.0;
   }
   else if(i == 0) {
      return (f[i+1][j][k]-f[i][j][k]);
   }
   else {
      return (f[i][j][k]-f[i-1][j][k]);
   }
}

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

static double gy(int i, int j, int k) {
   if(j > 0 && j < nym1) {
      return (f[i][j+1][k]-f[i][j-1][k]) / 2.0;
   }
   else if(j == 0) {
      return (f[i][j+1][k]-f[i][j][k]);
   }
   else {
      return (f[i][j][k]-f[i][j-1][k]);
   }
}

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

static double gz(int i, int j, int k) {
   if(k > 0 && k < nzm1) {
      return (f[i][j][k+1]-f[i][j][k-1]) / 2.0;
   }
   else if(k == 0) {
      return (f[i][j][k+1]-f[i][j][k]);
   }
   else {
      return (f[i][j][k]-f[i][j][k-1]);
   }   
}

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

static int calc() {

/*  I label local cell vertices as in ref. above 
    and assigns x y z names like this :
    
                 y
                 |

                 4----------5
                /.         /|
               7----------6 |
               | .        | |
               | .        | |
               | 0........|.1  --x
               |/         |/
               3----------2
             
             /
            z
            
*/         
   register int i, j, k;
   int n;
   double fac;
   double xu, yu, zu, xui, xuij, yuj;
   int cubeindex;
   int e;
   int vertlist[12];
   int m1, m2, m3;
   double ux, uy, uz, vx, vy, vz, wx, wy, wz, s;
   double xnu, ynu, znu;

   /* edges are labelled from each corner:
      x+ then z+ then y+
      for k varying fastest
      then j varying fastest
      then i varying fastest
   */
   n = 0;
   for(i=0; i<nx; i++) {
      xui = i*dxx1;
      for(j=0; j<ny; j++) {
         xuij = xui+j*dyx1;
         yuj = j*dyy1;
         for(k=0; k<nz; k++) {
            xu = xuij+k*dzx1;
            yu = yuj+k*dzy1;
            zu = k*dzz1;
            xnu = gx(i, j, k);
            ynu = gy(i, j, k);
            znu = gz(i, j, k);
            if(i<nxm1) {
               if((isovalue-f[i][j][k] < 0.0 && isovalue-f[i+1][j][k] >= 0) ||
                  (isovalue-f[i][j][k] >= 0.0 && isovalue-f[i+1][j][k] < 0)) {
                  fac = (isovalue-f[i][j][k])/(f[i+1][j][k]-f[i][j][k]);
                  itab[n] = n2;
                  xs[n2] = xu + fac*dxx1;
                  ys[n2] = yu;
                  zs[n2] = zu;
                  ux = xnu + fac*(gx(i+1, j, k) - xnu);
                  uy = ynu + fac*(gy(i+1, j, k) - ynu);
                  uz = znu + fac*(gz(i+1, j, k) - znu);
                  ux = ux*exx+uy*eyx+uz*ezx;
                  uy =        uy*eyy+uz*ezy;
                  uz =               uz*ezz;
                  s = sqrt(ux*ux+uy*uy+uz*uz);
                  if(s <= 0.0) {
                     (void)fprintf(stderr, 
                        "\aWARNING: normal direction assumed\n");
                     xns[n2] = 1.0;
                     yns[n2] = 0.0;
                     zns[n2] = 0.0;
                  }
                  else {
                     xns[n2] = ux/s;
                     yns[n2] = uy/s;
                     zns[n2] = uz/s;
                  }
                  n2++;
               }
               else {
                  itab[n] = -1;
               }
            }
            else {
               itab[n] = -1;
            }
            n++;
            if(k<nzm1) {
               if((isovalue-f[i][j][k] < 0.0 && isovalue-f[i][j][k+1] >= 0) ||
                  (isovalue-f[i][j][k] >= 0.0 && isovalue-f[i][j][k+1] < 0)) {
                  fac = (isovalue-f[i][j][k])/(f[i][j][k+1]-f[i][j][k]);
                  itab[n] = n2;
                  xs[n2] = xu + fac*dzx1;
                  ys[n2] = yu + fac*dzy1;
                  zs[n2] = zu + fac*dzz1;
                  ux = xnu + fac*(gx(i, j, k+1) - xnu);
                  uy = ynu + fac*(gy(i, j, k+1) - ynu);
                  uz = znu + fac*(gz(i, j, k+1) - znu);
                  ux = ux*exx+uy*eyx+uz*ezx;
                  uy =        uy*eyy+uz*ezy;
                  uz =               uz*ezz;
                  s = sqrt(ux*ux+uy*uy+uz*uz);
                  if(s <= 0.0) {
                     (void)fprintf(stderr, 
                        "\aWARNING: normal direction assumed\n");
                     xns[n2] = 1.0;
                     yns[n2] = 0.0;
                     zns[n2] = 0.0;
                  }
                  else {
                     xns[n2] = ux/s;
                     yns[n2] = uy/s;
                     zns[n2] = uz/s;
                  }
                  n2++;
               }
               else {
                  itab[n] = -1;
               }
            }
            else {
               itab[n] = -1;
            }
            n++;
            if(j<nym1) {
               if((isovalue-f[i][j][k] < 0.0 && isovalue-f[i][j+1][k] >= 0) ||
                  (isovalue-f[i][j][k] >= 0.0 && isovalue-f[i][j+1][k] < 0)) {
                  fac = (isovalue-f[i][j][k])/(f[i][j+1][k]-f[i][j][k]);
                  itab[n] = n2;
                  xs[n2] = xu + fac*dyx1;
                  ys[n2] = yu + fac*dyy1;
                  zs[n2] = zu;
                  ux = xnu + fac*(gx(i, j+1, k) - xnu);
                  uy = ynu + fac*(gy(i, j+1, k) - ynu);
                  uz = znu + fac*(gz(i, j+1, k) - znu);
                  ux = ux*exx+uy*eyx+uz*ezx;
                  uy =        uy*eyy+uz*ezy;
                  uz =               uz*ezz;
                  s = sqrt(ux*ux+uy*uy+uz*uz);
                  if(s <= 0.0) {
                     (void)fprintf(stderr, 
                        "\aWARNING: normal direction assumed\n");
                     xns[n2] = 1.0;
                     yns[n2] = 0.0;
                     zns[n2] = 0.0;
                  }
                  else {
                     xns[n2] = ux/s;
                     yns[n2] = uy/s;
                     zns[n2] = uz/s;
                  }
                  n2++;
               }
               else {
                  itab[n] = -1;
               }
            }
            else {
               itab[n] = -1;
            }
            n++;
         }
      }
   }
   
   if(n2 == 0) {
      (void)fprintf(stderr, "\aERROR: no isosurfaces found\n");
      return 1;
   }

   for(i=0; i<nxm1; i++) {
      for(j=0; j<nym1; j++) {
         for(k=0; k<nzm1; k++) {
            cubeindex = 0;
            if (f[i  ][j  ][k  ] < isovalue) cubeindex |=   1;
            if (f[i+1][j  ][k  ] < isovalue) cubeindex |=   2;
            if (f[i+1][j  ][k+1] < isovalue) cubeindex |=   4;
            if (f[i  ][j  ][k+1] < isovalue) cubeindex |=   8;
            if (f[i  ][j+1][k  ] < isovalue) cubeindex |=  16;
            if (f[i+1][j+1][k  ] < isovalue) cubeindex |=  32;
            if (f[i+1][j+1][k+1] < isovalue) cubeindex |=  64;
            if (f[i  ][j+1][k+1] < isovalue) cubeindex |= 128;
   
            e = edgeTable[cubeindex];

            if (e == 0) continue;
               
            n = i*nyz3+j*nz3+k*3; 

            if (e & 1   ) vertlist[0 ] = itab[n];
            if (e & 2   ) vertlist[1 ] = itab[n + nyz3 + 1];
            if (e & 4   ) vertlist[2 ] = itab[n + 3];
            if (e & 8   ) vertlist[3 ] = itab[n + 1];
            if (e & 16  ) vertlist[4 ] = itab[n + nz3];
            if (e & 32  ) vertlist[5 ] = itab[n + nyz3 + 1 + nz3];
            if (e & 64  ) vertlist[6 ] = itab[n + 3 + nz3];
            if (e & 128 ) vertlist[7 ] = itab[n + 1 + nz3];
            if (e & 256 ) vertlist[8 ] = itab[n + 2];
            if (e & 512 ) vertlist[9 ] = itab[n + nyz3 + 2];
            if (e & 1024) vertlist[10] = itab[n + nyz3 + 3 + 2];
            if (e & 2048) vertlist[11] = itab[n + 3 + 2];
            
            for (n = 0; triTable[cubeindex][n]!= -1; n += 3) {
               m1 = vertlist[triTable[cubeindex][n  ]];
               if(m1 == -1) {
                  (void)fprintf(stderr, "\aERROR m1 %d %d %d\n", i, j, k);
              return 1;
               }
               m2 = vertlist[triTable[cubeindex][n+1]];
               if(m2 == -1) {
                  (void)fprintf(stderr, "\aERROR m1 %d %d %d\n", i, j, k);
                  return 1;
               }
               m3 = vertlist[triTable[cubeindex][n+2]];
               if(m3 == -1) {
                  (void)fprintf(stderr, "\aERROR m1 %d %d %d\n", i, j, k);
              return 1;
               }
               ux = xs[m2] - xs[m1];
               uy = ys[m2] - ys[m1];
               uz = zs[m2] - zs[m1];   
               vx = xs[m3] - xs[m1];
               vy = ys[m3] - ys[m1];
               vz = zs[m3] - zs[m1];   
               wx = uy*vz-uz*vy;
               wy = uz*vx-ux*vz;
               wz = ux*vy-uy*vx;
               s = sqrt(wx*wx+wy*wy+wz*wz);
               if(s <= 0.0) continue;
               (void)fprintf(out, "3 %d %d %d\n", m1+1, m2+1, m3+1);
               n1++;
            }
         }
      }
   }
   
   if(n1 == 0) {
      (void)fprintf(stderr, "\aERROR: no isosurface found\n");
      return 1;
   }
   
   
   for(n=0; n<n2; n++) {
      (void)fprintf(out, "%g %g %g  %g %g %g\n",
         xs[n], ys[n], zs[n], xns[n], yns[n], zns[n]);
   }

   return 0;

}

int pot2surf_direct(const gchar *surf_file_to_write, const gchar *pot_file_to_read, int nsurfs_to_build, 
                const float *surf_value, const gchar **surf_name)
{
   register int i, j, k;
   char ttt[37];
   char rep[MAX_LINE_LENGTH], info[MAX_LINE_LENGTH];
   double v;
   char uuu[37];
   long int posi1, posi2;

   n1_tot = 0;
   n2_tot = 0;

   in = fopen(pot_file_to_read,"r");
   if (!in) {
      (void)fprintf(stderr, "\aERROR: cannot read potential file %s\n", pot_file_to_read);
      return 1;
   }   
   DBG_fprintf(stderr, "pot2surf : Reading potential file %s\n", pot_file_to_read);
   
   out = fopen(surf_file_to_write,"w+");
   if (!out) {
      (void)fprintf(stderr, 
         "\aERROR: cannot write+ surfaces in file %s\n", surf_file_to_write);
      return 1;
   }
   DBG_fprintf(stderr, "pot2surf : Writing surface file %s ...\n", surf_file_to_write);

   /* 1st line (comment) */
   (void)fgets(rep, MAX_LINE_LENGTH, in);
   rep[strlen(rep)-1] = 0; /* skipping \n */
   DBG_fprintf(stderr, "Pot2surf : read comment line '%s'\n", rep);
   (void)fprintf(out, "%s\n", rep);
   
   (void)fgets(rep, MAX_LINE_LENGTH, in);
   (void)sscanf(rep, "%d %d %d", &nx, &ny, &nz);
   
   f = g_malloc(nx*sizeof(double **));
   for(i=0; i<nx; i++) {
      f[i] = g_malloc(ny*sizeof(double *));
      for(j=0; j<ny; j++) {
         f[i][j] = g_malloc(nz*sizeof(double));
      }
   }
   
   nxm1 = nx - 1;
   nym1 = ny - 1;
   nzm1 = nz - 1;
   
   nz3 = 3*nz;
   nyz3 = ny*nz3;
   nxyz3 = nx*nyz3;
   
   itab = g_malloc(nxyz3*sizeof(int));
   
   xs = g_malloc(nxyz3*sizeof(double));
   ys = g_malloc(nxyz3*sizeof(double));
   zs = g_malloc(nxyz3*sizeof(double));
   xns = g_malloc(nxyz3*sizeof(double));
   yns = g_malloc(nxyz3*sizeof(double));
   zns = g_malloc(nxyz3*sizeof(double));
   
   (void)fgets(rep, MAX_LINE_LENGTH, in);
   (void)sscanf(rep, "%lf %lf %lf", &dxx, &dyx, &dyy);
   (void)fgets(rep, MAX_LINE_LENGTH, in); 
   (void)sscanf(rep, "%lf %lf %lf", &dzx, &dzy, &dzz);
   
   (void)fprintf(out, "%f %f %f\n", dxx, dyx, dyy);
   (void)fprintf(out, "%f %f %f\n", dzx, dzy, dzz);

   v = dxx*dyy*dzz;
   exx = dyy*dzz/v;
   eyx = -dyx*dzz/v;
   ezx = (dyx*dzy-dyy*dzx)/v;
   eyy = dxx*dzz/v;
   ezy = -dxx*dzy/v;
   ezz = dxx*dyy/v;   
   
   dxx1 = dxx/nx;
   dyx1 = dyx/ny;
   dyy1 = dyy/ny;
   dzx1 = dzx/nz;
   dzy1 = dzy/nz;
   dzz1 = dzz/nz;
   
   (void)fgets(rep, MAX_LINE_LENGTH, in);
   if(sscanf(rep, "%s", info) != 1) {
      (void)fprintf(stderr, "\aERROR: missing chain xyz OR zyx\n");
      return 1;
   }
   if(!strcmp(info, "xyz")) {   
      for ( k = 0; k < nz; k++ ) 
         for ( j = 0; j < ny; j++ ) 
            for ( i = 0; i < nx; i++ ) 
               (void)fscanf(in, "%lf", &f[i][j][k]);
   }
   else if(!strcmp(info, "zyx")) {   
      for ( i = 0; i < nx; i++ ) 
         for ( j = 0; j < ny; j++ ) 
            for ( k = 0; k < nz; k++ ) 
               (void)fscanf(in, "%lf", &f[i][j][k]);
   }
   else {
      (void)fprintf(stderr, 
         "\aERROR: 5th line should contain chain xyz OR zyx\n");
      return 1;
   } 
   
   (void)sprintf(ttt, "%12d%12d%12d", nsurfs_to_build, n1_tot, n2_tot);
   (void)fprintf(out, "%s\n", ttt);
   
   for(k=0; k<nsurfs_to_build; k++) {
   
     isovalue = surf_value[k];

      DBG_fprintf(stderr, "Building isosurface for isovalue = %f (%s)\n", 
              isovalue, surf_name[k]);
      fprintf(out, "%s %f\n", ISOSURFACES_FLAG_POTENTIAL, isovalue);
      (void)fprintf(out, "%s\n", surf_name[k]);
      
      n1 = 0;
      n2 = 0;
      (void)fflush(out);
      posi1 = ftell(out);
      (void)sprintf(uuu, "%12d%12d", n1, n2);
      (void)fprintf(out, "%s\n", uuu);
      
      if(calc() != 0)
      return 1;
      
      DBG_fprintf(stderr, "  Found %d facets and %d points\n", n1, n2);
      (void)fflush(out);
      posi2 = ftell(out);
      (void)fseek(out, posi1, SEEK_SET);
      (void)fflush(out);
      (void)sprintf(uuu, "%12d%12d", n1, n2);
      (void)fprintf(out, "%s\n", uuu);
      (void)fseek(out, posi2, SEEK_SET);
      (void)fflush(out);
      
      n1_tot += n1;
      n2_tot += n2;
   }
   
   (void)fflush(out);
   (void)rewind(out);
   (void)fgets(rep, MAX_LINE_LENGTH, out);
   (void)fgets(rep, MAX_LINE_LENGTH, out);
   (void)fgets(rep, MAX_LINE_LENGTH, out);
   (void)fflush(out);
   (void)sprintf(ttt, "%12d%12d%12d", nsurfs_to_build, n1_tot, n2_tot);
   (void)fprintf(out, "%s\n", ttt);
   (void)fflush(out);
               
   (void)fclose(out);
   
   DBG_fprintf(stderr, 
      "Found a total of %d facets and %d points\n", n1_tot, n2_tot);
   DBG_fprintf(stderr, "Done.\n"); 


   for(i=0; i<nx; i++) 
     {
       for(j=0; j<ny; j++) 
       g_free(f[i][j]);
       g_free(f[i]);
     }
   g_free(f);   
   
   g_free(itab);
   
   g_free(xs);
   g_free(ys);
   g_free(zs);
   g_free(xns);
   g_free(yns);
   g_free(zns);

   
   return 0;
}  
/*int main(int argc, char *argv[]) {*/
/* 0 in case of success
   n != 0 otherwise */
int pot2surf(const char *file_to_convert, GString *file_created)
{
   FILE *ins;
   register int i, j, k;
   char ttt[37];
   int nsurf;
   char rep[MAX_LINE_LENGTH], repn[MAX_LINE_LENGTH], info[MAX_LINE_LENGTH];
   double v;
   char uuu[37];
   long int posi1, posi2;

   if(file_to_convert == NULL)
     {
       fprintf(stderr, "WARNING : file to convert is NULL\n");
       return 1;
     }
   
   n1_tot = 0;
   n2_tot = 0;

   ins = fopen(file_to_convert,"r");
   if (!ins) {
      (void)fprintf(stderr, "\aERROR: cannot read in file %s\n", file_to_convert);
      return 1;
   }

   (void)fgets(rep, MAX_LINE_LENGTH, ins);
   rep[strlen(rep)-1] = 0;
   in = fopen(rep,"r");
   if (!in) {
      (void)fprintf(stderr, "\aERROR: cannot read potential file %s\n", rep);
      return 1;
   }
   
   DBG_fprintf(stderr, "Reading potential file %s\n", rep);
   
   (void)fgets(rep, MAX_LINE_LENGTH, ins);
   rep[strlen(rep)-1] = 0;
   out = fopen(rep,"w+");
   if (!out) {
      (void)fprintf(stderr, 
         "\aERROR: cannot write+ surfaces in file %s\n", rep);
      return 1;
   }
   if(file_created != NULL)
     g_string_append(file_created, rep);
   DBG_fprintf(stderr, "Writing surface file %s ...\n", rep);

   /* 1st line (comment) */
   (void)fgets(rep, MAX_LINE_LENGTH, in);
   rep[strlen(rep)-1] = 0; /* skipping \n */
   (void)fprintf(out, "%s\n", rep);
   
   (void)fgets(rep, MAX_LINE_LENGTH, in);
   (void)sscanf(rep, "%d %d %d", &nx, &ny, &nz);
   
   f = g_malloc(nx*sizeof(double **));
   for(i=0; i<nx; i++) {
      f[i] = g_malloc(ny*sizeof(double *));
      for(j=0; j<ny; j++) {
         f[i][j] = g_malloc(nz*sizeof(double));
      }
   }
   
   nxm1 = nx - 1;
   nym1 = ny - 1;
   nzm1 = nz - 1;
   
   nz3 = 3*nz;
   nyz3 = ny*nz3;
   nxyz3 = nx*nyz3;
   
   itab = g_malloc(nxyz3*sizeof(int));
   
   xs = g_malloc(nxyz3*sizeof(double));
   ys = g_malloc(nxyz3*sizeof(double));
   zs = g_malloc(nxyz3*sizeof(double));
   xns = g_malloc(nxyz3*sizeof(double));
   yns = g_malloc(nxyz3*sizeof(double));
   zns = g_malloc(nxyz3*sizeof(double));
   
   (void)fgets(rep, MAX_LINE_LENGTH, in);
   (void)sscanf(rep, "%lf %lf %lf", &dxx, &dyx, &dyy);
   (void)fgets(rep, MAX_LINE_LENGTH, in); 
   (void)sscanf(rep, "%lf %lf %lf", &dzx, &dzy, &dzz);
   
   (void)fprintf(out, "%f %f %f\n", dxx, dyx, dyy);
   (void)fprintf(out, "%f %f %f\n", dzx, dzy, dzz);

   v = dxx*dyy*dzz;
   exx = dyy*dzz/v;
   eyx = -dyx*dzz/v;
   ezx = (dyx*dzy-dyy*dzx)/v;
   eyy = dxx*dzz/v;
   ezy = -dxx*dzy/v;
   ezz = dxx*dyy/v;   
   
   dxx1 = dxx/nx;
   dyx1 = dyx/ny;
   dyy1 = dyy/ny;
   dzx1 = dzx/nz;
   dzy1 = dzy/nz;
   dzz1 = dzz/nz;
   
   (void)fgets(rep, MAX_LINE_LENGTH, in);
   if(sscanf(rep, "%s", info) != 1) {
      (void)fprintf(stderr, "\aERROR: missing chain xyz OR zyx\n");
      return 1;
   }
   if(!strcmp(info, "xyz")) {   
      for ( k = 0; k < nz; k++ ) 
         for ( j = 0; j < ny; j++ ) 
            for ( i = 0; i < nx; i++ ) 
               (void)fscanf(in, "%lf", &f[i][j][k]);
   }
   else if(!strcmp(info, "zyx")) {   
      for ( i = 0; i < nx; i++ ) 
         for ( j = 0; j < ny; j++ ) 
            for ( k = 0; k < nz; k++ ) 
               (void)fscanf(in, "%lf", &f[i][j][k]);
   }
   else {
      (void)fprintf(stderr, 
         "\aERROR: 5th line should contain chain xyz OR zyx\n");
      return 1;
   } 
   
   (void)fgets(rep, MAX_LINE_LENGTH, ins);
   if(sscanf(rep, "%d", &nsurf) != 1) {
      (void)fprintf(stderr, 
         "\aERROR: cannot read in %s nbr of isosurfaces to build\n", file_to_convert);
      return 1;
   }
   
   (void)sprintf(ttt, "%12d%12d%12d", nsurf, n1_tot, n2_tot);
   (void)fprintf(out, "%s\n", ttt);
   
   for(k=0; k<nsurf; k++) {
   
      (void)fgets(rep, MAX_LINE_LENGTH, ins);
      if(sscanf(rep, "%lf %s", &isovalue, repn) != 2) {
         (void)fprintf(stderr, 
            "\aERROR: cannot read in %s #%d isovalue and name\n", file_to_convert,
            k+1);
         return 1;
      }

      DBG_fprintf(stderr, "Building isosurface for isovalue = %f (%s)\n", 
         isovalue, repn);
      fprintf(out, "%s %f\n", ISOSURFACES_FLAG_POTENTIAL, isovalue);
      (void)fprintf(out, "%s\n", repn);
      
      n1 = 0;
      n2 = 0;
      (void)fflush(out);
      posi1 = ftell(out);
      (void)sprintf(uuu, "%12d%12d", n1, n2);
      (void)fprintf(out, "%s\n", uuu);
      
      if(calc() != 0)
      return 1;
      
      (void)fprintf(stderr, "  Found %d facets and %d points\n", n1, n2);
      (void)fflush(out);
      posi2 = ftell(out);
      (void)fseek(out, posi1, SEEK_SET);
      (void)fflush(out);
      (void)sprintf(uuu, "%12d%12d", n1, n2);
      (void)fprintf(out, "%s\n", uuu);
      (void)fseek(out, posi2, SEEK_SET);
      (void)fflush(out);
     
      n1_tot += n1;
      n2_tot += n2;
   }
   
   (void)fflush(out);
   (void)rewind(out);
   (void)fgets(rep, MAX_LINE_LENGTH, out);
   (void)fgets(rep, MAX_LINE_LENGTH, out);
   (void)fgets(rep, MAX_LINE_LENGTH, out);
   (void)fflush(out);
   (void)sprintf(ttt, "%12d%12d%12d", nsurf, n1_tot, n2_tot);
   (void)fprintf(out, "%s\n", ttt);
   (void)fflush(out);
               
   (void)fclose(out);
   
   DBG_fprintf(stderr, 
      "Found a total of %d facets and %d points\n", n1_tot, n2_tot);
   DBG_fprintf(stderr, "Done.\n"); 


   for(i=0; i<nx; i++) 
     {
       for(j=0; j<ny; j++) 
       g_free(f[i][j]);
       g_free(f[i]);
     }
   g_free(f);   
   
   g_free(itab);
   
   g_free(xs);
   g_free(ys);
   g_free(zs);
   g_free(xns);
   g_free(yns);
   g_free(zns);

   
   return 0;
}  
   
/******************************************************************************/
/******************************************************************************/


/* DEVELOPMENT AREA
   Methods are copied to be callable without all static stuff everywhere. */

static double gradX(int i, int j, int k, double ***data, int n) {
   if(i > 0 && i < n) {
      return (data[i+1][j][k]-data[i-1][j][k]) / 2.0;
   }
   else if(i == 0) {
      return (data[i+1][j][k]-data[i][j][k]);
   }
   else {
      return (data[i][j][k]-data[i-1][j][k]);
   }
}

static double gradY(int i, int j, int k, double ***data, int n) {
   if(j > 0 && j < n) {
      return (data[i][j+1][k]-data[i][j-1][k]) / 2.0;
   }
   else if(j == 0) {
      return (data[i][j+1][k]-data[i][j][k]);
   }
   else {
      return (data[i][j][k]-data[i][j-1][k]);
   }
}

static double gradZ(int i, int j, int k, double ***data, int n) {
   if(k > 0 && k < n) {
      return (data[i][j][k+1]-data[i][j][k-1]) / 2.0;
   }
   else if(k == 0) {
      return (data[i][j][k+1]-data[i][j][k]);
   }
   else {
      return (data[i][j][k]-data[i][j][k-1]);
   }   
}

gboolean pot2surfCreate(Surfaces **surf, ScalarField *field,
                  double isoValue, int id, gchar *name)
{

/*  I label local cell vertices as in ref. above 
    and assigns x y z names like this :
    
                 y
                 |

                 4----------5
                /.         /|
               7----------6 |
               | .        | |
               | .        | |
               | 0........|.1  --x
               |/         |/
               3----------2
             
             /
            z
            
*/         
  register int i, j, k;
  int n;
  double fac;
  double xu, yu, zu, xui, xuij, yuj;
  int cubeindex;
  int e;
  int vertlist[12];
  int m1, m2, m3;
  double ux, uy, uz, vx, vy, vz, wx, wy, wz, s;
  double xnu, ynu, znu, v;
  int sizem1[3];
  int nelez3, neleyz3, nelexyz3;
  double orthoBox[6];
  double dBox[6];
  int nPolys, nPoints, nPolys_old, nPoints_old;
  double *xsurf, *ysurf, *zsurf, *xnsurf, *ynsurf, *znsurf;
  int *iTab, *vTab;
  gchar *nameSurf;
  float* densityData;
  double ***scalarData;
  float scalarBox[6];
  int scalarGrid[3];
  gboolean new;

  g_return_val_if_fail(surf, FALSE);

  /* edges are labelled from each corner:
     x+ then z+ then y+
     for k varying fastest
     then j varying fastest
     then i varying fastest
  */

  /* Get data from the field. */
  scalarData = scalarFieldGet_data(field);
  scalarFieldGet_box(field, scalarBox);
  scalarFieldGet_gridSize(field, scalarGrid);
  
  DBG_fprintf(stderr, "Pot2surf : compute isosurface at value %g from"
            " field %p.\n", isoValue, (gpointer)field);
  DBG_fprintf(stderr, " | size %d %d %d.\n", scalarGrid[0],
            scalarGrid[1], scalarGrid[2]);
  DBG_fprintf(stderr, " | box %g %g %g\n", scalarBox[0],
            scalarBox[1], scalarBox[2]);
  DBG_fprintf(stderr, " |     %g %g %g.\n", scalarBox[3],
            scalarBox[4], scalarBox[5]);

  /* Compute some local values. */
  for (i = 0; i < 3; i++)
    sizem1[i] = scalarGrid[i] - 1;
  nelez3 = 3 * scalarGrid[2];
  neleyz3 = scalarGrid[1] * nelez3;
  nelexyz3 = scalarGrid[0] * neleyz3;
  v   = scalarBox[0] * scalarBox[2] * scalarBox[5];
  orthoBox[0] = scalarBox[2] * scalarBox[5] / v;
  orthoBox[1] = -scalarBox[1] * scalarBox[5] / v;
  orthoBox[2] = scalarBox[0] * scalarBox[5] / v;
  orthoBox[3] = (scalarBox[1] * scalarBox[4] - scalarBox[2] * scalarBox[3]) / v;
  orthoBox[4] = -scalarBox[0] * scalarBox[4] / v;
  orthoBox[5] = scalarBox[0] * scalarBox[2] / v;   
  dBox[0] = scalarBox[0] / scalarGrid[0];
  dBox[1] = scalarBox[1] / scalarGrid[1];
  dBox[2] = scalarBox[2] / scalarGrid[1];
  dBox[3] = scalarBox[3] / scalarGrid[2];
  dBox[4] = scalarBox[4] / scalarGrid[2];
  dBox[5] = scalarBox[5] / scalarGrid[2];
  nPolys = 0;
  nPoints = 0;

  /* Allocate buffers. */
  /* Point values. */
  xsurf = g_malloc(nelexyz3 * sizeof(double));
  ysurf = g_malloc(nelexyz3 * sizeof(double));
  zsurf = g_malloc(nelexyz3 * sizeof(double));
  /* Normal values. */
  xnsurf = g_malloc(nelexyz3 * sizeof(double));
  ynsurf = g_malloc(nelexyz3 * sizeof(double));
  znsurf = g_malloc(nelexyz3 * sizeof(double));
  /* Indexes. */
  iTab = g_malloc(nelexyz3 * sizeof(int));
  vTab = g_malloc(nelexyz3 * sizeof(int));

   n = 0;
   for(i=0; i<scalarGrid[0]; i++) {
      xui = i*dBox[0];
      for(j=0; j<scalarGrid[1]; j++) {
         xuij = xui+j*dBox[1];
         yuj = j*dBox[2];
         for(k=0; k<scalarGrid[2]; k++) {
            xu = xuij+k*dBox[3];
            yu = yuj+k*dBox[4];
            zu = k*dBox[5];
            xnu = gradX(i, j, k, scalarData, sizem1[0]);
            ynu = gradY(i, j, k, scalarData, sizem1[1]);
            znu = gradZ(i, j, k, scalarData, sizem1[2]);
            if(i<sizem1[0]) {
               if((isoValue-scalarData[i][j][k] < 0.0 && isoValue-scalarData[i+1][j][k] >= 0) ||
                  (isoValue-scalarData[i][j][k] >= 0.0 && isoValue-scalarData[i+1][j][k] < 0)) {
                  fac = (isoValue-scalarData[i][j][k])/(scalarData[i+1][j][k]-scalarData[i][j][k]);
                  iTab[n] = nPoints;
                  xsurf[nPoints] = xu + fac*dBox[0];
                  ysurf[nPoints] = yu;
                  zsurf[nPoints] = zu;
                  ux = xnu + fac*(gradX(i+1, j, k, scalarData, sizem1[0]) - xnu);
                  uy = ynu + fac*(gradY(i+1, j, k, scalarData, sizem1[1]) - ynu);
                  uz = znu + fac*(gradZ(i+1, j, k, scalarData, sizem1[2]) - znu);
                  ux = ux*orthoBox[0]+uy*orthoBox[1]+uz*orthoBox[3];
                  uy =        uy*orthoBox[2]+uz*orthoBox[4];
                  uz =               uz*orthoBox[5];
                  s = sqrt(ux*ux+uy*uy+uz*uz);
                  if(s <= 0.0) {
                     (void)fprintf(stderr, 
                        "\aWARNING: normal direction assumed\n");
                     xnsurf[nPoints] = 1.0;
                     ynsurf[nPoints] = 0.0;
                     znsurf[nPoints] = 0.0;
                  }
                  else {
                     xnsurf[nPoints] = ux/s;
                     ynsurf[nPoints] = uy/s;
                     znsurf[nPoints] = uz/s;
                  }
                  nPoints++;
               }
               else {
                  iTab[n] = -1;
               }
            }
            else {
               iTab[n] = -1;
            }
            n++;
            if(k<sizem1[2]) {
               if((isoValue-scalarData[i][j][k] < 0.0 && isoValue-scalarData[i][j][k+1] >= 0) ||
                  (isoValue-scalarData[i][j][k] >= 0.0 && isoValue-scalarData[i][j][k+1] < 0)) {
                  fac = (isoValue-scalarData[i][j][k])/(scalarData[i][j][k+1]-scalarData[i][j][k]);
                  iTab[n] = nPoints;
                  xsurf[nPoints] = xu + fac*dBox[3];
                  ysurf[nPoints] = yu + fac*dBox[4];
                  zsurf[nPoints] = zu + fac*dBox[5];
                  ux = xnu + fac*(gradX(i, j, k+1, scalarData, sizem1[0]) - xnu);
                  uy = ynu + fac*(gradY(i, j, k+1, scalarData, sizem1[1]) - ynu);
                  uz = znu + fac*(gradZ(i, j, k+1, scalarData, sizem1[2]) - znu);
                  ux = ux*orthoBox[0]+uy*orthoBox[1]+uz*orthoBox[3];
                  uy =        uy*orthoBox[2]+uz*orthoBox[4];
                  uz =               uz*orthoBox[5];
                  s = sqrt(ux*ux+uy*uy+uz*uz);
                  if(s <= 0.0) {
                     (void)fprintf(stderr, 
                        "\aWARNING: normal direction assumed\n");
                     xnsurf[nPoints] = 1.0;
                     ynsurf[nPoints] = 0.0;
                     znsurf[nPoints] = 0.0;
                  }
                  else {
                     xnsurf[nPoints] = ux/s;
                     ynsurf[nPoints] = uy/s;
                     znsurf[nPoints] = uz/s;
                  }
                  nPoints++;
               }
               else {
                  iTab[n] = -1;
               }
            }
            else {
               iTab[n] = -1;
            }
            n++;
            if(j<sizem1[1]) {
               if((isoValue-scalarData[i][j][k] < 0.0 && isoValue-scalarData[i][j+1][k] >= 0) ||
                  (isoValue-scalarData[i][j][k] >= 0.0 && isoValue-scalarData[i][j+1][k] < 0)) {
                  fac = (isoValue-scalarData[i][j][k])/(scalarData[i][j+1][k]-scalarData[i][j][k]);
                  iTab[n] = nPoints;
                  xsurf[nPoints] = xu + fac*dBox[1];
                  ysurf[nPoints] = yu + fac*dBox[2];
                  zsurf[nPoints] = zu;
                  ux = xnu + fac*(gradX(i, j+1, k, scalarData, sizem1[0]) - xnu);
                  uy = ynu + fac*(gradY(i, j+1, k, scalarData, sizem1[1]) - ynu);
                  uz = znu + fac*(gradZ(i, j+1, k, scalarData, sizem1[2]) - znu);
                  ux = ux*orthoBox[0]+uy*orthoBox[1]+uz*orthoBox[3];
                  uy =        uy*orthoBox[2]+uz*orthoBox[4];
                  uz =               uz*orthoBox[5];
                  s = sqrt(ux*ux+uy*uy+uz*uz);
                  if(s <= 0.0) {
                     (void)fprintf(stderr, 
                        "\aWARNING: normal direction assumed\n");
                     xnsurf[nPoints] = 1.0;
                     ynsurf[nPoints] = 0.0;
                     znsurf[nPoints] = 0.0;
                  }
                  else {
                     xnsurf[nPoints] = ux/s;
                     ynsurf[nPoints] = uy/s;
                     znsurf[nPoints] = uz/s;
                  }
                  nPoints++;
               }
               else {
                  iTab[n] = -1;
               }
            }
            else {
               iTab[n] = -1;
            }
            n++;
         }
      }
   }
  DBG_fprintf(stderr, " | found %d points.\n", nPoints);
   
   if(nPoints == 0)
     {
       g_warning("no isosurface found.");
       g_free(xsurf);
       g_free(ysurf);
       g_free(zsurf);
       g_free(xnsurf);
       g_free(ynsurf);
       g_free(znsurf);
       g_free(iTab);
       g_free(vTab);
       return FALSE;
     }

   for(i=0; i<sizem1[0]; i++) {
      for(j=0; j<sizem1[1]; j++) {
         for(k=0; k<sizem1[2]; k++) {
            cubeindex = 0;
            if (scalarData[i  ][j  ][k  ] < isoValue) cubeindex |=   1;
            if (scalarData[i+1][j  ][k  ] < isoValue) cubeindex |=   2;
            if (scalarData[i+1][j  ][k+1] < isoValue) cubeindex |=   4;
            if (scalarData[i  ][j  ][k+1] < isoValue) cubeindex |=   8;
            if (scalarData[i  ][j+1][k  ] < isoValue) cubeindex |=  16;
            if (scalarData[i+1][j+1][k  ] < isoValue) cubeindex |=  32;
            if (scalarData[i+1][j+1][k+1] < isoValue) cubeindex |=  64;
            if (scalarData[i  ][j+1][k+1] < isoValue) cubeindex |= 128;
   
            e = edgeTable[cubeindex];

            if (e == 0) continue;
               
            n = i*neleyz3+j*nelez3+k*3; 

            if (e & 1   ) vertlist[0 ] = iTab[n];
            if (e & 2   ) vertlist[1 ] = iTab[n + neleyz3 + 1];
            if (e & 4   ) vertlist[2 ] = iTab[n + 3];
            if (e & 8   ) vertlist[3 ] = iTab[n + 1];
            if (e & 16  ) vertlist[4 ] = iTab[n + nelez3];
            if (e & 32  ) vertlist[5 ] = iTab[n + neleyz3 + 1 + nelez3];
            if (e & 64  ) vertlist[6 ] = iTab[n + 3 + nelez3];
            if (e & 128 ) vertlist[7 ] = iTab[n + 1 + nelez3];
            if (e & 256 ) vertlist[8 ] = iTab[n + 2];
            if (e & 512 ) vertlist[9 ] = iTab[n + neleyz3 + 2];
            if (e & 1024) vertlist[10] = iTab[n + neleyz3 + 3 + 2];
            if (e & 2048) vertlist[11] = iTab[n + 3 + 2];
            
            for (n = 0; triTable[cubeindex][n]!= -1; n += 3) {
               m1 = vertlist[triTable[cubeindex][n  ]];
               if(m1 == -1)
             {
               g_warning("m1 %d %d %d.", i, j, k);
               return FALSE;
             }
               m2 = vertlist[triTable[cubeindex][n+1]];
               if(m2 == -1)
             {
               g_warning("m2 %d %d %d.", i, j, k);
               return FALSE;
             }
               m3 = vertlist[triTable[cubeindex][n+2]];
               if(m3 == -1)
             {
               g_warning("m3 %d %d %d.", i, j, k);
               return FALSE;
             }
               ux = xsurf[m2] - xsurf[m1];
               uy = ysurf[m2] - ysurf[m1];
               uz = zsurf[m2] - zsurf[m1];   
               vx = xsurf[m3] - xsurf[m1];
               vy = ysurf[m3] - ysurf[m1];
               vz = zsurf[m3] - zsurf[m1];   
               wx = uy*vz-uz*vy;
               wy = uz*vx-ux*vz;
               wz = ux*vy-uy*vx;
               s = sqrt(wx*wx+wy*wy+wz*wz);
               if(s <= 0.0) continue;
             vTab[3 * nPolys + 0] = m1;
             vTab[3 * nPolys + 1] = m2;
             vTab[3 * nPolys + 2] = m3;
               nPolys++;
            }
         }
      }
   }
  DBG_fprintf(stderr, " | found %d polygons.\n", nPolys);
   
   if(nPolys == 0)
     {
       g_warning("no isosurface found.");
       g_free(xsurf);
       g_free(ysurf);
       g_free(zsurf);
       g_free(xnsurf);
       g_free(ynsurf);
       g_free(znsurf);
       g_free(iTab);
       g_free(vTab);
       return FALSE;
     }

   /* Ok, try to create a Surfaces object from here. */
   if (*surf == (Surfaces*)0)
     {
       DBG_fprintf(stderr, " | create new Surfaces.\n");
       *surf = isosurfacesNew();
       (*surf)->nsurf = 1;
       isosurfacesAllocate(*surf, (*surf)->nsurf, nPolys, nPoints);
       nPolys_old = 0;
       nPoints_old = 0;
       for (i = 0; i < 6; i++)
       (*surf)->local_box[i] = scalarBox[i];
       densityData = isosurfacesAdd_floatProperty(*surf, ISOSURFACES_PROPERTY_POTENTIAL);
     }
   else
     {
       DBG_fprintf(stderr, " | expand existing Surfaces.\n");
       nPolys_old = (*surf)->num_polys;
       nPoints_old = (*surf)->num_points;
       (*surf)->nsurf += 1;
       isosurfacesReallocate(surf, (*surf)->nsurf, nPolys_old + nPolys,
                       nPoints_old + nPoints);
       for (i = 0; i < 6; i++)
       if ((*surf)->local_box[i] != scalarBox[i])
         g_warning("box inconstitency between field and surface.");
       densityData = isosurfacesGet_floatProperty(*surf, ISOSURFACES_PROPERTY_POTENTIAL);
     }
   
  DBG_fprintf(stderr, " | copy vertices.\n");
   for(n = 0; n < nPoints; n++)
     {
       (*surf)->poly_points[nPoints_old + n][0] = xsurf[n];
       (*surf)->poly_points[nPoints_old + n][1] = ysurf[n];
       (*surf)->poly_points[nPoints_old + n][2] = zsurf[n];
       (*surf)->poly_normals[nPoints_old + n][0] = -xnsurf[n];
       (*surf)->poly_normals[nPoints_old + n][1] = -ynsurf[n];
       (*surf)->poly_normals[nPoints_old + n][2] = -znsurf[n];
     }
   for (n = 0; n < nPolys; n++)
     {
       (*surf)->poly_surf_index[nPolys_old + n] = (*surf)->nsurf - 1;
       (*surf)->poly_num_vertices[nPolys_old + n] = 3;
       (*surf)->poly_vertices[nPolys_old + n] =
       g_malloc(sizeof(int) * 3);
       for (i = 0; i < 3; i++)
       (*surf)->poly_vertices[nPolys_old + n][i] =
         vTab[3 * n + i] + nPoints_old;
     }

   /* Additional stuff for the given surface. */
   DBG_fprintf(stderr, " | finishing surfaces.\n");
   if (name)
     nameSurf = g_strdup(name);
   else
     nameSurf = g_strdup_printf(_("Surface %d"), id + 1);
   (*surf)->resources[(*surf)->nsurf - 1] =
     isosurfacesGet_resourceFromName(nameSurf, &new);
   if (new)
     (*surf)->resources[(*surf)->nsurf - 1]->rendered = TRUE;
   densityData[(*surf)->nsurf - 1] = isoValue;
   (*surf)->ids[(*surf)->nsurf - 1] = id;

#if DEBUG == 1
   isosurfacesCheck_consistency(*surf);
#endif

   g_free(xsurf);
   g_free(ysurf);
   g_free(zsurf);
   g_free(xnsurf);
   g_free(ynsurf);
   g_free(znsurf);
   g_free(iTab);
   g_free(vTab);
   return TRUE;
}

Generated by  Doxygen 1.6.0   Back to index