/* ReSCAL - genesis entry
 *
 * Copyright (C) 2011-2012
 *
 * Author: Olivier Rozier <rozier@ipgp.fr>
 *
 * This file is part of ReSCAL.
 *
 * ReSCAL is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */

#define _MAIN_

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include "defs.h"
#include "macros.h"
#include "genesis.h"
#include "cells.h"
#include "param.h"
#include "format.h"


int prog=PROG_TOOL;
int  H, L, P, HL, HLP;        // size of cellular space
Cell *TE=NULL;          // cellular space
//int  node = -1;            // node number (with PARALLEL option)
char *model_name = NULL; //name of model
char *bin_filename = NULL; //name of BIN file
char *csp_filename = NULL; //name of CSP file
char opt_bin = 0;   //raw binary format

char *csp_template_str = NULL; //value of CSP_template parameter
CSP_Template csp_template; //CSP template
CSP_Template t_templates[MAX_TEMPLATES]; //available templates
int nb_templates = 0; //number of templates

char *boundary_str = NULL;
int boundary = BC_PERIODIC;  //boundary conditions

double drand48();
void  srand48();

float* load_surface(char *, int, int);

void init_template(int type, char *name, char *desc, int nb_args, ...)
{
  va_list vl;
  int i;

  t_templates[nb_templates].type = type; 
  t_templates[nb_templates].name = name;
  t_templates[nb_templates].desc = desc;
  t_templates[nb_templates].nb_args = nb_args;
  
  AllocMemory(t_templates[nb_templates].args, float, nb_args);
  va_start(vl, nb_args);
  for (i=0; i<nb_args; i++){
    t_templates[nb_templates].args[i] = (float) va_arg(vl, double);
    //LogPrintf("t_templates[%d].args[%d] = %f\n", nb_templates, i, t_templates[nb_templates].args[i]);
  }
  va_end(vl);
  nb_templates++;
}

// TEMPLATES USER'S GUIDE
// ----------------------
//
// A CSP template is a predefined cellular space configuration, with optional arguments.
// The user gives to the parameter Csp_template a value of the form:
// TEMPLATE_NAME, or TEMPLATE_NAME(arg1, ...).
//
// Example: LAYER(5) will generate a sand layer of 5 cells height
//
// Here we explain how to create a new CSP template:
// 1) add a new type in CSP_TEMPLATES.
// 2) define the new template in init_template_list() by adding a line of the form
//    init_template(template type, "TEMPLATE NAME", "Template usage description", number of optional arguments, default value for argument 1, default value for argument 2, ...);
//    Note: default values of arguments MUST be floats (in case of an integer, just add a zero fractional part).
// 3) in genesis() function, write the implementation for the new CSP template.
//

//available CSP template types
#ifdef MODELE_DUN
enum CSP_TEMPLATES {CSP_CUSTOM, CSP_LAYER, CSP_BLOCK, CSP_CYLINDER, CSP_CONE, CSP_RCONE, CSP_CONE3, CSP_CONE5, CSP_RCONE5, CSP_RWALL, CSP_WAVES_2D, CSP_SMILEY};
#else
enum CSP_TEMPLATES {CSP_CUSTOM};
#endif

//initialization of available templates
void init_template_list()
{
  init_template(CSP_CUSTOM, "CUSTOM", "CUSTOM: no template", 0);
#ifdef MODELE_DUN
  init_template(CSP_LAYER, "LAYER", "LAYER(h=1.0): sand layer of height <h>", 1, 1.0);
  init_template(CSP_BLOCK, "BLOCK", "BLOCK(xmin, xmax, ymin, ymax, h): sand block of height <h>", 5, 30.0, 50.0, 50.0, 150.0, 20.0);
  init_template(CSP_CYLINDER, "CYLINDER", "CYLINDER(h=20, w=100): cylinder of height <h> and width <w>", 2, 20.0, 100.0);
  init_template(CSP_CONE, "CONE", "CONE(h=0, w=100): cone of height <h> and width <w>", 2, 30.0, 100.0);
  init_template(CSP_RCONE, "RCONE", "RCONE(rh=0.66, rw=0.4): cone of relative height <rh>*H and relative width <rw>*Min(L,D)", 2, 0.66, 0.4);
  init_template(CSP_CONE3, "CONE3", "CONE3(h=40, w=126): 3 cones of height <h> and width <w>", 2, 40.0, 126.0);
  init_template(CSP_CONE5, "CONE5", "CONE5(h=40, w=126): 5 cones of height <h> and width <w>", 2, 40.0, 126.0);
  init_template(CSP_RCONE5, "RCONE5", "RCONE5(rh=0.66, rw=0.2): 5 cones of relative height <rh>*H and relative width <rw>*Min(L,D)", 2, 0.66, 0.2);
  init_template(CSP_RWALL, "RWALL", "RWALL(rh=0.5, rx=0.3): wall of relative height <rh>*H and relative position <rx>*L", 2, 0.5, 0.3);
  init_template(CSP_WAVES_2D, "WAVES_2D", "WAVES_2D: 2d sand waves with various periods for stability analysis", 2, 0.5, 0.3);
  init_template(CSP_SMILEY, "SMILEY", "SMILEY: :-)", 0);
#endif
  LogPrintf("init_template_list: %d templates\n", nb_templates);
}

int parse_template()
{
  char *ptr;
  int i;

  //read template name
  ptr = strtok(csp_template_str,"(");
  LogPrintf("csp_template.name = %s\n", ptr);

  //find and copy the full template
  for (i=0; (i < nb_templates) && strcmp(ptr,t_templates[i].name); i++);
  if (i < nb_templates){
    memcpy(&csp_template, &t_templates[i], sizeof(CSP_Template));
    LogPrintf("csp_template.type = %d\n", csp_template.type);
    AllocMemory(csp_template.args, float, csp_template.nb_args);
    memcpy(csp_template.args, t_templates[i].args, csp_template.nb_args*sizeof(float));
  }
  else{
    ErrPrintf("ERROR: Incorrect value for Csp_template, %s\n", csp_template_str);
    exit(-2);
  }

  //read template optional arguments, if any
  ptr = strtok(NULL,",)");
  i = 0;
  while (ptr && (i < csp_template.nb_args)){
    //LogPrintf("csp_template.args[%d] = %s\n", i, ptr);
    csp_template.args[i++] = read_float(ptr);
    ptr = strtok(NULL,",)");
  }
  if (ptr){
    ErrPrintf("ERROR: Too many arguments in Csp_template (%s)\n", ptr);
    exit(-2);
  }
  for (i=0; i < csp_template.nb_args; i++){
    LogPrintf("csp_template.args[%d] = %f\n", i, csp_template.args[i]);
  }

  return i;
}


void genesis_parse()
{
  if (model_name && strcmp(model_name, MOD_NAME)){
    WarnPrintf("ERROR: Wrong model name %s in parameter file. Check model in defs.h and rebuild all.\n", model_name);
    exit(-2);
  }
  
  if (csp_template_str){
    parse_template();
  }

  if (boundary_str){
    if (!strcmp(boundary_str, "PERIODIC")){
#ifdef CYCLAGE_HOR
      boundary = BC_PERIODIC;
#else
      ErrPrintf("Incorrect value for Boundary : %s, please define CYCLAGE_HOR in defs.h and recompile\n", boundary_str);
      exit(-2);
#endif
    }
    else if (!strcmp(boundary_str, "OPEN")){
      boundary = BC_OPEN;
    }
    else if (!strcmp(boundary_str, "REINJECTION")){
      boundary = BC_REINJECTION;
    }
    else{
      ErrPrintf("Incorrect value for Boundary : %s\n", boundary_str);
      exit(-2);
    }
    //pbc_mode = (boundary == BC_PERIODIC);
  }
}

void genesis()
{
  int i, j, k, Hd2, Ld2, Pd2, Pd4, Hd4, Hd3;
  Cell *aux;
  
  HL = H*L;
  HLP = HL*P;

  // allocation of the cellular space
  AllocMemoryPrint("TE",TE, Cell, HLP);
  ResetMemory(TE, Cell, HLP);

  // center of the cellular space
  Hd2 = (int) H/2;
  Ld2 = (int) L/2;
  Pd2 = (int) P/2;
  
  Pd4 = (int) P/4;
  Hd4 = (int) H/4;
  Hd3 = (int) H/3;

  aux = TE;

#ifdef MODELE_DUN
  float di, dj, dk;
  int Ld3, Ld34, Ld4, Ld8;
  float Ldx, Ldy, Ldz;
  float pente = 2.0;
  float h, w, rc;
  int xmin, xmax, ymin, ymax;

  LogPrintf("modele DUNES\n");

//// normal mode ////

  Ld3 = (int) L/3;
  Ld4 = (int) L/4;
  Ld8 = (int) L/8;
  Ld34= (int) 3*L/4;

  //Ldx= Ld2 + 0.5; //Ld4; //10+Hd2;
  Ldx= Ld2-0.5; //Ld2-0.5; //Ld3;//Ld4;
  Ldy= /*0.66**/Pd2-0.5; //Ld4; //(int) L/2;
  Ldz= Hd2;

  int ll1 = (int) L*0.2; //L*0.15; //L*0.5;//L*0.1; //10; //0.1 //0.22; //Ld3
  int ll2 = (int) L*0.2; //0.22; //Ld3
  int l0 = (int) L*0.1; //L*0.1; //(L-ll1-ll2)*0.4;//25; //0.45
  int l1 = (int) l0+ll1*1.3; //pour 2 triangles
  int l2 = (int)(L-ll1-ll2)*0.25; //0.45
 //int l1 = (int)(L-l0-l2-ll1-ll2);
  int hh = (int) (H/2);//(100/3); //8+(100/3); //(H/3);//(60/2); //H/6; //H/2;
  int H0 = H-2; //H-8
  int hh1 = hh/1.5;
  int hh2 = hh;
  int lc = P; //1; //400; //Ld4; //Ld2; //L/10; //Ld4; //largeur du couloir
  int lcone = lc*0.6; //lc*0.6; //lc/2; //lc*0.6; //lc*2/3;  //Ld2*0.8; //largeur du cone
  int jcone = H0-lcone/2; //sommet du cone
  int hcone = 2*Hd3; //Hd2;//2*Hd3; //H-10; //hauteur du cone
 //int h0 = (int) (H-hh)/2;
 //int h2 = H-h0;
  int rcyl = Ld2-10; //rayon cylindre
  int hgrain = Hd2; //epaisseur de grains
  float periode = 10.0; //periode des ondulations
  float kk, kkk;
  // initialization of the cellular space
  for(k=0; k < P; k++){ // profondeur
 //   periode = 10.0+(k-(Pd2-lc/2))/4;
    periode = 10.0+k/2;
    kkk = k+9;
    if (k > 0 & k < P/10) kk =  L/46;
    else if (k > P/10 & k < P/10*2 ) kk =  L/42+1;
    else if (k > P/10*2 & k < P/10*3) kk =  L/38;
    else if (k > P/10*3 & k < P/10*4) kk =  L/34;
    else if (k > P/10*4 & k < P/10*5) kk =  L/30;
    else if (k > P/10*5 & k < P/10*6) kk =  L/26;
    else if (k > P/10*6 & k < P/10*7) kk =  L/22;
    else if (k > P/10*7 & k < P/10*8) kk =  L/18;
    else if (k > P/10*8 & k < P/10*9) kk =  L/14;
    else  kk =  L/10;

    //periode += 0.5*drand48();
    for(j=0; j < H; j++){ // hauteur
      for(i=0; i < L; i++, aux++){ //largeur
        di = i - Ldx;
        dj = j - Ldz;
        dk = k - Ldy;

        /*** EAUC cells ***/
        aux->celltype = EAUC; //default state
        
        /*** GR cells ***/
        if (csp_template.type == CSP_LAYER){
          //CSP_LAYER: sand layer
          //format: LAYER(h)
          hh = (int)csp_template.args[0];
          float frac = fmodf(csp_template.args[0],1.0);
          if (j > H-2-hh) aux->celltype = GR;
          if ((j == H-2-hh) && frac){ //h is not an integer
            if (drand48() < frac) aux->celltype = GR;
          }
        }
        if (csp_template.type == CSP_BLOCK){
          //CSP_BLOCK: sand block
          //format: BLOCK(xmin, xmax, ymin, ymax, h)
          xmin = (int)csp_template.args[0];
          xmax = (int)csp_template.args[1];
          ymin = (int)csp_template.args[2];
          ymax = (int)csp_template.args[3];
          hh = (int)csp_template.args[4];
          if ((i >= xmin) && (i < xmax) && (k >= ymin) && (k < ymax) && (j > H-2-hh)){
            aux->celltype = GR;
          }
        }
        else if (csp_template.type == CSP_CYLINDER){
          //CSP_CONE: sand cylinder
          //format: CYLINDER(h, w)
          h = csp_template.args[0];
          w = csp_template.args[1];
          if ((j >= (H-h-1)) && (di*di + dk*dk < 0.25*w*w)) aux->celltype = GR;
        }
        else if (csp_template.type == CSP_CONE){
          //CSP_CONE: sand cone
          //format: CONE(h, w)
          hcone = csp_template.args[0];
          lcone = csp_template.args[1];
          rc = 0.5*lcone*(j-(H-hcone))/(hcone-1);
          if ((j > (H-hcone)) && (di*di + dk*dk < rc*rc)) aux->celltype = GR;
        }
        else if (csp_template.type == CSP_RCONE){
          //CSP_RCONE: sand cone with relative dimensions
          //format: RCONE(rh, rw)
          hcone = csp_template.args[0]*H;
          lcone = csp_template.args[1]*Min(L,P);
          rc = 0.5*lcone*(j-(H-hcone))/(hcone-1);
          if ((j > (H-hcone)) && (di*di + dk*dk < rc*rc)) aux->celltype = GR;
          //if (j > H-3) aux->celltype = GR; //uniform layer
        }
        else if (csp_template.type == CSP_CONE3){
          //CSP_CONE3: 3 sand cones with fixed dimensions
          //format: RCONE3(h, w)
          hcone = csp_template.args[0];
          lcone = csp_template.args[1];
          rc = 0.5*lcone*(j-(H-hcone))/(hcone-1);
          float di0 = i - Ldx - 0.5*Ldx*cos(0);
          float dk0 = k - Ldy - 0.5*Ldy*sin(0);
          float di2 = i - Ldx - 0.5*Ldx*cos(2*PI/3);
          float dk2 = k - Ldy - 0.5*Ldy*sin(2*PI/3);
          float di4 = i - Ldx - 0.5*Ldx*cos(4*PI/3);
          float dk4 = k - Ldy - 0.5*Ldy*sin(4*PI/3);
          if ((j > (H-hcone)) && (di0*di0 + dk0*dk0 < rc*rc)) aux->celltype = GR;
          if ((j > (H-hcone)) && (di2*di2 + dk2*dk2 < rc*rc)) aux->celltype = GR;
          if ((j > (H-hcone)) && (di4*di4 + dk4*dk4 < rc*rc)) aux->celltype = GR;
        }
        else if ((csp_template.type == CSP_CONE5) || (csp_template.type == CSP_RCONE5)){
          if (csp_template.type == CSP_CONE5){
            //CSP_CONE5: 5 sand cones with fixed dimensions
            //format: RCONE5(h, w)
            hcone = csp_template.args[0];
            lcone = csp_template.args[1];
          }
          else{
            //CSP_RCONE5: 5 sand cones with relative dimensions
            //format: RCONE5(rh, rw)
            hcone = csp_template.args[0]*H;
            lcone = csp_template.args[1]*Min(L,P);
          }
          rc = 0.5*lcone*(j-(H-hcone))/(hcone-1);
          float di1 = i - Ldx - 0.6*Ldx*cos(PI/5);
          float dk1 = k - Ldy - 0.6*Ldy*sin(PI/5);
          float di3 = i - Ldx - 0.6*Ldx*cos(3*PI/5);
          float dk3 = k - Ldy - 0.6*Ldy*sin(3*PI/5);
          float di5 = i - Ldx - 0.6*Ldx*cos(5*PI/5);
          float dk5 = k - Ldy - 0.6*Ldy*sin(5*PI/5);
          float di7 = i - Ldx - 0.6*Ldx*cos(7*PI/5);
          float dk7 = k - Ldy - 0.6*Ldy*sin(7*PI/5);
          float di9 = i - Ldx - 0.6*Ldx*cos(9*PI/5);
          float dk9 = k - Ldy - 0.6*Ldy*sin(9*PI/5);
          if ((j > (H-hcone)) && (di1*di1 + dk1*dk1 < rc*rc)) aux->celltype = GR;
          if ((j > (H-hcone)) && (di3*di3 + dk3*dk3 < rc*rc)) aux->celltype = GR;
          if ((j > (H-hcone)) && (di5*di5 + dk5*dk5 < rc*rc)) aux->celltype = GR;
          if ((j > (H-hcone)) && (di7*di7 + dk7*dk7 < rc*rc)) aux->celltype = GR;
          if ((j > (H-hcone)) && (di9*di9 + dk9*dk9 < rc*rc)) aux->celltype = GR;
          //if ((j > (H-hcone)) && (di*di + dk*dk < rc*rc)) aux->celltype = GR; //central cone
          //if (j > H-3) aux->celltype = GR; //uniform layer
        }
        else if (csp_template.type == CSP_RWALL){
          //CSP_RWALL(rh,rx): north-south wall or DUM cells
          int h = csp_template.args[0]*H;
          int x = csp_template.args[1]*L;
          if ((i>=x-1) && (i<=x+1) && (j>=H-h)) aux->celltype = DUM;
        }
        else if (csp_template.type == CSP_WAVES_2D){
          //CSP_WAVES_2D: 2d sand waves with various periods for stability analysis
         // if (j > Hd2-1.5*sin(kk*i*2*PI/L)) aux->celltype = GR; //variable ondulations
          if (j > Hd2-2*sin(i*2*PI/periode)) aux->celltype = GR; //variable ondulations

         // if ((k%(P/10) == 0) || (k > 0 & k <1.1) || (k > P-1.01 & k <P) || (k==P) || (j > H0) || (j==0)) aux->celltype = DUM;
        // if ((k & 1) || (j > H0) || (j==0)) aux->celltype = DUM; //ground + ceiling + walls east-west
         if ((k % 2==0) || (j > H0) || (j==0)) aux->celltype = DUM; //ground + ceiling + walls east-west

        }
        else if (csp_template.type == CSP_SMILEY){
          //CSP_SMILEY: :-)
          if (j<H/2 && ((cos(PI*(i-L/2)/L/2)*cos(PI*(k-P/2)/P/2)<0.8 && cos(PI*(i-L/2)/L/2)*cos(PI*(k-P/2)/P/2)>0.75)
             || (cos(2*PI*(i-L/2)/L)<0.4 && sin(PI/4+2*PI*(k-P/2)/P)<0.4 && cos(2*PI*(i-L/2)/L)>0.2 && sin(PI/4+2*PI*(k-P/2)/P)>0.2 && k<L/2)
             || (cos(PI*(i-L/2)/L)*cos(PI*(k-P/2)/P)<0.8 && cos(PI*(i-L/2)/L)*cos(PI*(k-P/2)/P)>0.75&&k>3*P/5)))
               aux->celltype = GR;
        }
        else{
          //CSP_CUSTOM
          //if (0) //pas de grain
            // if ((j > H-0.4*H) && (i <= 0.2*L) && (i >= 0.025*L))
           //  if ((i >= 0.15*L) && (j > H-0.25*H))
          //if ((di*di + dk*dk < (j-(H-hcone))*(j-(H-hcone))) && (j > (H-hcone))) //cone de hauteur hcone
          //if ((di*di + dk*dk < 4*(j-2*Hd3)*(j-2*Hd3)) && (j > 2*Hd3)) //cone de hauteur Hd3 et de rayon 2*Hd3
          //if ((di*di + dk*dk < 2.5*(j-Hd3)*(j-Hd3)) && (j > Hd3)) //cone de hauteur 2*Hd3 pour les dunes etoiles
          //if ((di*di + dk*dk < (j-jcone)*(j-jcone)) && (j > Hd2)) //cone tronque
          //if (((di*di + dk*dk < (j-jcone)*(j-jcone)) && (j > jcone)) || (j > H-10)) //cone + couche uniforme
          //if ((j > 12) && (di*di + dk*dk < 400)) // disque
          //if ((j > H*0.65) && (fabs(di) < Ld8) && (fabs(dk) < Ld8)) // parallelepipede rectangle
          //if ((j > Hd2) && ((di+dk)*(di+dk) < 400) && (dk*dk < 400)) // parallelepipede 'tordu'
          //if (j + (i/20)> H-10) //pente EO
          //if (j + k*pente > H-20) //pente NS
          //if (j + i > H) //pente EO
         // if (j > H-3) //couche uniforme
          //if (j>H-8-(drand48()-0.90)) //couche aleaoire
          //if ((j > H-4) && (k>=Pd4) && (k<P-Pd4)) //couche uniforme est-ouest
          //if ((j > H-5) && (di*di + dk*dk < rcyl*rcyl)) //couche uniforme circulaire
          //if (/*(i>1) &&*/ (j > H-Hd2*drand48())) //couche non-uniforme
          //if ((j > H-2) || ((j == H-2) && i < L*0.5 && (drand48()>0.75/*0.5*/))) //couche non-uniforme
//if ((j > H-0.05*H) && (i <= 0.9*L) && (i >= 0.1*L))
     // if (j >H-0.5*H)(j > H0) || (j==0) 
        //if ((dj*dj + dk*dk < 0.42*H*0.42*H) && (dj*dj + dk*dk > 0.40*H*0.40*H) || (j > H0) || (j==0)) 
          if ((dj*dj + dk*dk > 0.48*H*0.48*H)) 
 //     if (i+0.05*H-(H-0.95*j) > 0 && j >H-0.5*H)
      //if (j > Hd2-2*sin(i*2*PI/20)) //ondulations
     //  if (((di+L*0.3)*(di+L*0.3) + dk*dk < 2*(j-1.4*Hd3)*(j-1.4*Hd3)) && (j > 1.4*Hd3)) //barchan dunes
 //  if (j > H-0.5*H && i > 0.05*L && i < 0.4*L) 
          //if (j > Hd2-2*sin(i*2*PI/periode)) //ondulations variables
          //if (((i>=l0) && (i<l0+1/*l0*/) && (j>H-hh)) || (j>=H-1)) //rectangle
          //if (((i>=l0) && (i<l0+1) && (j>H-hh)) || (j< (H-hh-10)*(Ld2-k)/lc) || (j>=H-1)) //rectangle avec toit incline
         // if (((i>=l0) && (i<l0+ll1) && ((i-l0)-((H0-j/*+drand48()*5*/)*ll1/hh) > 0))
           //    || (j>=H-40)) //triangle
      //    if (((i>=l0) && (i<l0+ll1) && ((i-l0)-((H0-j/*+drand48()*5*/)*ll1/hh) > 0))
        //            || ((i>=l0+ll1) && ((i-l0-ll1)-(hh-(H0-j))/3 < 0))
          //        || (j>=H0-35)) //triangle + triangle inverse
          //if (((i>=l0) && (i<l0+ll1) && ((i-l0)-((H-j)*ll1/hh1) > 0))
          //      || ((i>=l1) && (i<l1+ll2) && ((i-l1)-((H-j)*ll2/hh2) > 0))
          //      || (j>=H-1)) //2 triangles
          /*if (((i>=l0) && (i<=l0+ll1) && ((i-l0)-((H-j)*ll1/hh) > 0)) || (j< H*((k<Ld2) ? Ld2-k : k-Ld2)/Ld2)
            || (j>=H-1)) //triangle avec toit incline*/
                          /*else if (((i>=l0) && (i<=l0+ll1) && ((i-l0)-((H-j)*ll1/hh) > 0)) || (j< Hd2*((k<Ld2) ? Ld2-k : k-Ld2)/ll)
            || (j>=H-1)) //triangle avec toit en biseau*/
          //if (((i>=l0) && (i<=l0+ll1) && ((i-l0)-((H-j)*ll1/hh) > 0)) || (j< (H-hh-10)*(Ld2-k)/lc) || (j>=H-1)) //triangle avec toit incline
          /*if (((i>=l0) && (i<l0+ll1) && ((i-l0)-((H-j)*ll1/hh) > 0))
            || ((i>=l0+ll1) && (i<L-ll2-l2) && (j>=hh))
            || ((i>=L-ll2-l2) && (i<L-l2) && ((L-l2-i)-((H-j)*ll2/hh) > 0))
            || (j>=H-1)) //trapeze*/
          //if (di*di + dj*dj + dk*dk < Hd4*Hd4) // sphere
          //if ((di*di + dj*dj + dk*dk < Hd4*Hd4) || ((di+50)*(di+50) + (dj+30)*(dj+30) + dk*dk < Hd4*Hd4/4)) // 2 spheres
          //if (j >H0-3*(1+cos(10*i*k*PI/L))) //test mh
          //if (j >H0-3*(1+cos(10*i*k*PI/L/P))) //test mh
          //if (j >H0-3*(0+cos(10*i*PI/L)-cos(10*k*PI/P))) // champs de bosses
    

     {
        aux->celltype = DUM;
       // aux->celltype   = GR;
         }
        }

        /*** DUM cells ***/
  	    //if (0)
        //if (j > H-2) //plateau
     //if ((j==0) || (j > H0)) //plateau + plafond
           if ((j > H-0.05*H) && (dj*dj + dk*dk <=0.48*H*0.48*H))
	// if ((j==0) || (j > H-0.5*H)) 
        //if ((j > H0) || (j==0) || ((i==Ld3) && (j>=Hd2))) //plateau + plafond + mur
        //if ((j > H0) || (j==0) || (k>=P-5) || (k<5)) //plateau + plafond + couloir
        //if ((j > H0) || (j<4)) //plateau + plafond epais (4)
        //if ((k < Pd2-lc/2) || (k > Pd2+lc/2) || (k & 1) || (j > H0) || (j==0)) //plateau + plafond + couloir + peigne
        //if ((j > H0-5*(1+cos(10*i*PI/L))) || (j==0)) //plateau ondule + plafond
        //if ((j > H0-5*(1+cos(10*i*PI/L))) || (j < 4)) //plateau ondule + plafond epais (4)
        //if (((i == L-1) && (j>H-L))) //couloir + mur
        //if ((j > H-2) && (di*di + dk*dk < Ld2*Ld2)) // disque + plateau + plafond
        //if (((j > hgrain) && (di*di + dk*dk >= rcyl*rcyl)) || (j > H0) || (j==0)) // plateau + trou circulaire + plafond
        //if (j + (i/20)> H-2) //pente EO
        //if ((j + (i/20)> H-2) || (j==0)) //pente EO + plafond
        //if (j + k*pente > H+20) //pente NS
          aux->celltype   = GR;
         // aux->celltype = DUM;

#ifdef CELL_COLOR      
        /*** color of cells ***/
       if (aux->celltype == GR){
          //aux->color = 1; //unformly colored
          aux->color = (drand48() < 0.2) ? 1 : 0; //randomly colored
          //aux->color = (k > Pd2) ? 1 : 0; //colored area to the north
          //aux->color = (k == Pd2+1) && (j == H-20) && (i == 40) ? 1 : 0; //one colored cell for TRACE_PAR_COL
          //aux->color = (k == Pd2+1) && (i == 40) ? H-1-j : 0; //column of colored cells for TRACE_PAR_COL
        }
#endif

        /*** IN cells ***/
       // if (boundary == BC_REINJECTION){
        // if ((j > H-0.15*H && j < H-0.1*H && i < 0.009*L && i > 0.003*L) ||(j==1 && i >0.01*L))
           //  if (j==1 && i >0.3*L)
          // if ((i > 0.015*L) &&(i < 0.025*L) && (j > 0) && (j < 2)) 
      //     if ((0.015*L<i) && (i < 0.025*L) && (j < H-0.35*H) && (j > H-0.4*H))
         //  if ((k == Ld2) && (i == Ld2) && (j == 0)) //source de grains en haut au centre
//          if (/*(node == 0) &&*/ (i == 5) && (j == 0)) //source de grains en haut a gauche
          //if ((j == 0) && (di*di + dk*dk < 5*5)) ////source de grains en disque
          //if ((j == 0) && ((di*di + dk*dk >= rcyl*rcyl) && ((di+1)*(di+1) + dk*dk < rcyl*rcyl))) ////source de grains en demi-cercle
          //if ((j == 0) && (dk > -5) && (dk < 5) && ((di*di + dk*dk >= rcyl*rcyl) && ((di+1)*(di+1) + dk*dk < rcyl*rcyl))) ////source de grains en arc de cercle
          //if ((j == 0) && (((k >= i) && (i+k <= L) && (di*di + dk*dk >= rcyl*rcyl) && ((di+1)*(di+1) + dk*dk < rcyl*rcyl))
          //     || (((k <= i) || (i+k >= L)) && (i+1==round(Ld2-rcyl/sqrt(2)))))) //source de grains en quart de cercle + segments de droite
          //  aux->celltype = IN;
       // }

        /*** OUT cells ***/
      //if ((i > L-7) && (j > H0))

        //if ((j-1 == hgrain) && (di*di + dk*dk >= rcyl*rcyl)) ////sortie de grains en cercle
        //if ((j-1 == hgrain) && ((di>0) && ((di-1)*(di-1) + (dk-sgn(dk))*(dk-sgn(dk)) >= rcyl*rcyl))) ////sortie de grains en demi-cercle
        // aux->celltype = OUT;

      
      }
    }
  }


//// precomputed mode (height map) ////
/*
  float* surf;
  int lx = 800; //799;
  int ly = 400; //399;
  int cx = L/2;  //lx/2;
  int cy = P/2; //L*0.4;  //L/3; //L/2;
  int ly1 = (L-ly)/2;
  int ly2 = ly1+ly;
  float alpha = 0; //-PI/6.0;

  surf = load_surface("data/scale_dune_4.txt", lx, ly);
  //surf = load_surface("data/trunc_dune_2.txt", lx, ly);

  for(k=0; k < P; k++){ // profondeur
    for(j=0; j < H; j++){ // hauteur
      for(i=0; i < L; i++, aux++){ //largeur
        // rotation
        int i0 = (i-cx)*cos(alpha) + (k-cy)*sin(alpha) + cx;
        int k0 = (i-cx)*sin(alpha) - (k-cy)*cos(alpha) + cy;
        if (0)
        //if (k>L*0.8)
        //if ((k<L*0.4) || (k>L*0.6))
        //if ((k<Ld2) || (k>Ld2))
          aux->celltype = DUM;
        else if ((k0>=ly1) && (k0<ly2) && (i0>=0) && (i0<lx) && (j > H-surf[i0+(k0-ly1)*lx]))
          aux->celltype = GR;
        else
          //aux->celltype = EAUC;
          //aux->celltype = (drand48() < 0.92) ? EAUC : EAUT;
          //aux->celltype = (j>Hd2) ? EAUC : EAUT;
          aux->celltype = (j>0)||(i>2) ? EAUC : EAUT;
       }
    }
  }
*/
#endif

#ifdef MODELE_AVA
  int Ld34, Ld4, Ld3;
  int hc = H*0.4;
  float Ldx, Ldy;
  float di, dk;
  float tg = tan(PI/6);

  LogPrintf("modele AVA\n");

  Ld4 = (int) L/4;
  Ld3 = (int) L/3;
  Ld34= (int) 3*L/4;

  Ldx= Ld2;//5 + 0.5; //Ld4; //10+Hd2;
  Ldy= Pd2;//5 + 0.5; //Ld4; //(int) L/2;

  for(k=0; k < P; k++){ // profondeur
    for(j=0; j < H; j++){ // hauteur
      for(i=0; i < L; i++, aux++){ //largeur
        di = i - Ldx;
        dk = k - Ldy;
        if (di<0) di = -di;
        if (dk<0) dk = -dk;

        //if (0)
        if (j > H-2) //plateau
        //if ((j > H-2) || ((j >= H - 20) && (di*di + dk*dk < 200))) //plateau et socle
           aux->celltype   = DUM;
        //else if ((i==0) && (k==0) && (j==0)) //source de grains
        //else if ((i==Ld2) && (k==Ld2) && (j==0)) //source ponctuelle de grains
        else if ((di*di + dk*dk <= 25) && (j==0)) //source etendue de grains
        //else if ((((i-Ld3)*(i-Ld3) + dk*dk <= 25) || ((i-2*Ld3)*(i-2*Ld3) + dk*dk <= 25)) && (j==0)) //2 sources etendues de grains
        //else if ((j==0) && (di*di + dk*dk < 55*55) && ((di-1)*(di-1) + (dk-1)*(dk-1) >= 45*45)) //source circulaire de grains
          aux->celltype = IN;
        //else if ((di*di + dk*dk <= 25) && (j==H-2)) //sortie circulaire de grains
        //  aux->celltype = OUT;
        //else if ((j>H/4) && (di*di + dk*dk <= (Ld2-5)*(Ld2-5))) //cylindre
        //else if ((di*di + dk*dk < 2*(j-Hd3)*(j-Hd3)) && (j > Hd3)) //cone de hauteur 2*Hd3
        //else if ((di*di + dk*dk < (j-H+hc)*(j-H+hc)/(tg*tg)) && (j > H-hc)) //cone de hauteur hc et de pente tg
        //    aux->celltype   = DUM;
        else{
          aux->celltype   = AIR;
          //aux->celltype   = GR;
#ifdef CELL_COLOR
          //aux->color = (drand48()<0.5) ? 1 : 0;
#endif
        //if ((j == H-1) && (di*di + dk*dk >= (Ld2-5)*(Ld2-5))) ////sortie de grains en cercle
        //  aux->celltype = OUT;
        }
      }
    }
  }
#endif

#ifdef MODELE_CMB
  LogPrintf("modele CMB\n");
  for(k=0; k < P; k++) // profondeur
    for(j=0; j < H; j++) // hauteur
      for(i=0; i < L; i++, aux++){ //largeur
 //Modele classique
        if (j < Hd2)
        //else if ((j < Hd2) || ((j < Hd2+20) && (i > Ld2-10) && (i < Ld2+10)))
          aux->celltype   = PLUS;
        else
          aux->celltype   = MOINS;
        //aux->celltype   = (int) floor(drand48()*3);

/* Modele experimental 1 (2 plaques)*/
/*
        int L1 = L/3;
        int L2 = L*2/3;
        int H1 = H/3;
        int H2 = H*2/3;
        if ((i>L1) && (i<L2) && (k>L1) && (k<L2) && ((j<=H1) || (j>=H2)))
          aux->celltype = PLUS;
        else if ((i==L1) && (k>=L1) && (k<=L2) && ((j<=H1+5) || (j>=H2-5)))
          aux->celltype = DUM;
        else if ((i==L2) && (k>=L1) && (k<=L2) && ((j<=H1) || (j>=H2)))
          aux->celltype = DUM;
        else if ((k==L1) && (i>=L1) && (i<=L2) && ((j<=H1+5) || (j>=H2-5)))
          aux->celltype = DUM;
        else if ((k==L2) && (i>=L1) && (i<=L2) && ((j<=H1) || (j>=H2)))
          aux->celltype = DUM;
        else if (drand48() < 0.5)
          aux->celltype = ZERO;
        else
          aux->celltype = MOINS;
*/
// Modele experimental 2 (sillon diagonal)
 /*
        int D = H/6;
        int L1 = (L/4)-D;
        int L2 = (L/4)+D;
        int L3 = (3*L/4)-D;
        int L4 = (3*L/4)+D;
        int H1 = Hd2-D;
        int H2 = Hd2+D;
        int ii = i+k-Ld2; //fenetre courante
        if (ii<0) ii += L;
        if (ii>=L) ii -= L;
        if (ii<L1)
          aux->celltype = (j<H2) ? PLUS : MOINS;
        else if (ii<L2)
          aux->celltype = (j<H2-(ii-L1)) ? PLUS : MOINS;
        else if (ii<L3)
          aux->celltype = (j<H1) ? PLUS : MOINS;
        else if (ii<L4)
          aux->celltype = (j<H1+(ii-L3)) ? PLUS : MOINS;
        else
          aux->celltype = (j<H2) ? PLUS : MOINS;

  */
/* Modele experimental 3 (1 faille oblique)*/
/*
        int L1 = (L-Hd2)*0.9;
        int L2 = (L-Hd2)*1.1;
        int ii = i+j; // fenetre courante
        if ((ii>L1) && (ii<L2))
          aux->celltype = MOINS;
        else
          aux->celltype = PLUS;
*/

/* Modele experimental 4 (2 fractures auto-affines)*/
/*
        //premiere fracture oblique
        static float* surf1 = NULL;
        static float* surf2 = NULL;
        if (!surf1) surf1 = load_surface("surf256.dat", 256, 256);
        if (!surf2) surf2 = load_surface("surf256b.dat", 256, 256);
        float alpha = 70*PI/180;
        float amp = 1.0;

        int H1 = Hd2-5;
        int H2 = Hd2+5;

        // rotation
        int i0 = (i-Ld2)*cos(alpha) + (j-Hd2)*sin(alpha) + Ld2;
        int j0 = (i-Ld2)*sin(alpha) - (j-Hd2)*cos(alpha) + Hd2;
        if (i0 < 0) i0 += 256;
        if ((j0 > H1+amp*surf1[i0*256+k]) && (j0 < H2-amp*surf2[i0*256+k]))
          aux->celltype = ZERO;
        else
          aux->celltype = PIERRE;

        //deuxieme fracture verticale
        static float* surf3 = NULL;
        static float* surf4 = NULL;
        if (!surf3) surf3 = load_surface("surf256c.dat", 256, 256);
        if (!surf4) surf4 = load_surface("surf256d.dat", 256, 256);
        amp = 1.0;

        int L1 = (L/3)-5;
        int L2 = (L/3)+5;

        if ((k > L1+amp*surf3[i+j*256]) && (k < L2-amp*surf4[i+j*256]))
          aux->celltype = ZERO;
*/
      }
#endif

#ifdef MODELE_ICB
  LogPrintf("modele ICB\n");
  for(k=0; k < P; k++) // profondeur
    for(j=0; j < H; j++) // hauteur
      for(i=0; i < L; i++, aux++){ //largeur
        if (j < H*2/3)
          aux->celltype   = (drand48() < 0.5)?ZERO:MOINS;
        else
          aux->celltype   = PLUS;
        //aux->celltype   = (int) floor(drand48()*3);
      }
#endif

#ifdef MODELE_RIV
  LogPrintf("modele RIVIERE\n");

// sillon diagonal unique parallele au flux, en haut
  for(k=0; k < P; k++) // profondeur
    for(j=0; j < H; j++) // hauteur
      for(i=0; i < L; i++, aux++){ //largeur
        int L1 = L*0.3;
        int L2 = L*0.7;
        int H1 = H*0.25;
        int H2 = H*0.4;
        int ii = i-k+Ld2; //fenetre courante
        if (ii<0) ii += L;
        if (ii>=L) ii -= L;
        if ((ii<L1) || (ii>L2))
          aux->celltype = (j<H1) ? EAU : TERRE;
        else
          aux->celltype = (j<H2) ? EAU : TERRE;
        if ((j==0) && (i==0) && (k==0)) aux->celltype = IN;  //source ponctuelle
     }


//modele precalcule (altimap)
#if 0
     float* surf;
     //int lx = 501;
     //int ly = 529;
     //int ly1 = (L-ly)/2;
     //int ly2 = ly1+ly;
     int lx = 512; //128;
     int ly = 512; //128;
     int ly1 = 0;
     int ly2 = ly-1; //ly*0.48;
     int ll = Ld2;
     float niv = H*0.9;
     //float niv = H-1;
     float amp = 30.0;

     //surf = load_surface("data/guade.txt", lx, ly);
     //surf = load_surface("../../surface/surf128.dat", lx, ly);
     //surf = load_surface("surf512_05.txt", lx, ly);

     for(k=0; k < P; k++){ // profondeur
       for(j=0; j < H; j++){ // hauteur
         for(i=0; i < L; i++, aux++){ //largeur
           //if ((k>=ly1) && (k<ly2) && (i<lx) && (j > H-surf[i+(k-ly1)*lx]/10.0))
           //if ((k>=ly1) && (k<ly2) && (i<lx) && (j > H-surf[i+(k-ly1)*lx]*amp))
           //if ((k>=ly1) && (k<ly2) && (i<lx) && (j > (H*0.4)-surf[i+(k-ly1+(L/2))*lx]*2.0))
           //if ((k>=ly1) && (k<ly2) && (i<lx) && (j > niv-((L-1-k)/2)-surf[i+(k-ly1)*lx]*2.0))
           //if ((k>=ly1) && (k<ly2) && (i<lx) && (j > niv-((L-1-k)/2)-surf[i+(k-ly1)*lx]/2.0)) //pente NS avec rugosite
           //if ((k>=ly1) && (k<ly2) && (i<lx) && (j > niv-((L-1-k)/2)-((i<Ld2)?(Ld2-i):(i-Ld2)))) //pente NS avec sillon
           //  aux->celltype = PLUS; //(j<H-1) ? PLUS : DUM;
           /*else if ((k>=ly2) && (i<lx) && (j > niv)) //plaine
           aux->celltype = PLUS;*/
           if (j > niv-((L-1-k)*0.9)) //pente de 0.9
           //if (j > niv-((L-1-k)*0.8)+H*sin(5.0*i*6.28/L)/20) //pente de 0.8 + sinusoide
           //if ((j > niv-(k*0.9)) && (j > niv-((L-1-k)*0.9))) //triangle de pente 0.9
             aux->celltype = TERRE;
             //aux->celltype = (((int) floor(drand48()*100) == 0)) ? PIERRE : TERRE;
           else
             aux->celltype = EAU;
           //if ((i==0) || (k==0) || (i==L-1) || (k==L-1)) aux->celltype = DUM;  //barriere infranchissable tout autour (modele ferme)
           if ((i==0) || (k==0) || (i==L-1) || ((k==P-1) && (j>=niv))) aux->celltype = DUM;  //barriere infranchissable tout autour et jusqu'a une certaine hauteur sur le bord sud
           //if ((i==0) || (i==L-1) || (((k==0) || (k==L-1)) && (j>=niv))) aux->celltype = DUM;  //barriere infranchissable tout autour et jusqu'a une certaine hauteur sur les bords sud et nord
           //if ((i==0) || (i==L-1) || (((k==0) || (k==L-1)) && ((j>=niv) || (i!=Ld2)))) aux->celltype = DUM;  //barriere infranchissable tout autour et jusqu'a une certaine hauteur au centre des bords sud et nord
           //if ((i<Ld2-(ll/2)) || (i>Ld2+(ll/2))) aux->celltype = DUM; //couloir
           //if ((i==Ld2) && (k==L/10) && (j==0)) aux->celltype = IN; //source
         }
       }
     }
#endif
#endif

#ifdef MODELE_RIV2
  LogPrintf("modele RIVIERE-2\n");
//#define WAVES
#ifdef WAVES
#define NB_WAV 6
  float offsetx[NB_WAV];
  float offsetz[NB_WAV];
  float alphaw[NB_WAV];
  float Hwave[NB_WAV];
  float Fwave[NB_WAV];
  int iw;
  for (iw=0; iw<NB_WAV; iw++){
    offsetx[iw] = drand48()*10;
    offsetz[iw] = drand48()*10;
    alphaw[iw] = drand48();
    LogPrintf("alphaw[%d] = %f\n", iw, alphaw[iw]);
    //Hwave[iw] = (iw==0)? (float)H*0.2 : Hwave[iw-1]/1.2;
    //Fwave[iw] = (iw==0)? (float)6.28/L : Fwave[iw-1]*1.5;
    Hwave[iw] = (iw < 2)? (float)H*0.2 : Hwave[iw-2]/1.5;
    Fwave[iw] = (iw < 2)? (float)6.28/L : Fwave[iw-2]*2.0;
  }
#endif
  for(k=0; k < P; k++){ // profondeur
    for(i=0; i < L; i++){ //largeur
      //float alpha =  1.0 - (float)(i+k)/(L+L); //pente diagonale
      float alpha =  1.0 - (float)(k)/(P);  //pente nord-sud
      int H1 = 1; //H*0.7; //H*0.1
      int H2 = H*0.8; //H*0.4
      //int H1 = Hd2*1.1;
      //int H2 = Hd2*1.7;
      //float hh = H1*alpha + H2*(1-alpha); //sol en pente
      float hh = H2;//H*0.3;
      //float hh = (i<Ld2)? H*0.5 : H*0.9;//H*0.3; //falaise
      /*int x1 = L/6;
      int x2 = L*5/6;
      int z1 = L/6;
      int z2 = L*5/6;
      if ((i>=x1) && (i<x2) && (k>=z1) && (k<z2)) hh = H*0.1; //rectangle sureleve*/
      //float hh = H*0.2; //H*0.4;
      /*float dd = ((L/2) - i)*((L/2) - i) + ((L/2) - k)*((L/2) - k);
      if (dd < 3000) // bosse centrale;
        hh -= (3000-dd)/120;*/
      /*int ii = i-k+Ld2; //fenetre pour le sillon
      int L1 = L*0.1;
      int L2 = L*0.5;
      if ((ii>L1) && (ii<L2)) //sillon diagonal
        hh += 3;*/
#ifdef WAVES
      hh = H*0.5;
      for (iw=1; iw<NB_WAV; iw++){
        hh += Hwave[iw] * cos(offsetx[iw]+(float)(alphaw[iw]*i+(1-alphaw[iw])*k)*Fwave[iw]) * cos(offsetz[iw]+(float)((1-alphaw[iw])*i-alphaw[iw]*k)*Fwave[iw]); //ondulations de surface
      }
#endif
      //hh += H*0.1*drand48();  //bruit blanc
      //int L34 = L*0.75;
      //if ((i-L34)*(i-L34)+(k-L34)*(k-L34)<100) //encoche circulaire
      //  hh=H;
      //if ((i+k<L/2) || ((L-i+k)<L/2) || ((L-k+i)<L/2))  //carre en diagonal
      //  aux->celltype = DUM;
      //else if (i+k>3*L/2)
      //  aux->celltype = MOINS;
      //else
      float niv = H*0.9;
      for(j=0; j < H; j++){ // hauteur
        aux = TE + k*HL + j*L + i;
        //aux->celltype = (j<hh) ? MOINS : PLUS;  //(j<H-1) ? PLUS : DUM;
        aux->celltype = ((k<P*0.48) && (j > niv-((P-1-k)/2))) ? PLUS : DUM;
        //aux->celltype = (j<hh) ? MOINS : ((int) floor(drand48()*20) > 0) ? PLUS : PIERRE;
        if ((j==0) && (i==Ld2) && (k==1)) aux->celltype = IN;  //source ponctuelle
        //if (((i==0) || (k==0)) && (j>=niv)) aux->celltype = DUM;  //barriere infranchissable jusqu'a une certaine hauteur sur le bord
        if (((i==L-1) || (k==P-1)) && (j>=niv)) aux->celltype = DUM;  //barriere infranchissable jusqu'a une certaine hauteur sur le bord
        if ((i==0) || (k==0) || (i==L-1) || ((k==P-1) && (j>=niv))) aux->celltype = DUM;  //barriere infranchissable tout autour et jusqu'a une certaine hauteur sur le bord sud
        //if (((i==L-1) || (k==L-1)) && (i+k<L+L-10)) aux->celltype = DUM;  //petite ouverture en coin
        //if ((i==0) || (k==0) || (i==L-1) || ((k==L-1) && ((i<Ld2-2) || (i>Ld2+2)))) aux->celltype = DUM;  //petite ouverture au milieu
      }
    }
  }

#endif

#ifdef MODELE_CHE
  LogPrintf("modele CHENAL\n");
// pente diagonale
  for(k=0; k < P; k++) // profondeur
    for(j=0; j < H; j++) // hauteur
      for(i=0; i < L; i++, aux++){ //largeur
        float alpha =  1.0 - (float)(i+k)/(L+P); //pente diagonale
        //float alpha =  1.0 - (float)(k)/(L);  //pente nord-sud
        int H1 = H*0.2;
        int H2 = H*0.8;
        int L2 = L*0.5;
        //float hh = H1*alpha + H2*(1-alpha); //sol en pente
        float hh = H*0.5; //sol plat
        //float dd = ((L/2) - i)*((L/2) - i) + ((L/2) - k)*((L/2) - k);

        //if (dd < 400) // bosse centrale;
        //  hh -= (400-dd)/40;

        //if (fabs((float)(k-i)) < 10) // chenal;
        //  hh += 2;

        //hh -= (float)(H/10)*cos((float)i/10)*cos((float)k/10); //ondulations de surface

        if ((j==H1-2) && /*(i>=L2-1) && (i<=L2+1) && (k>=L2-1) && (k<=L2+1)*/ (i<1) && (k<1))
        //if (j==0) //pluie
          aux->celltype = IN;
        //else if ((i>=L-10) || (k>=L-10)) //fosse !!
        //  aux->celltype = AIR;
        else //if (i+j+k>=H1)
          aux->celltype = (j<hh) ? AIR : PLUS;
          //aux->celltype = (j<hh) ? AIR : (((int) floor(drand48()*5) == 0)) ? PIERRE : PLUS;
        //else
        //  aux->celltype = MOINS;
      }
#endif

#ifdef MODELE_DUM
  LogPrintf("modele DUMMY\n");
// pour les tests ...
  for(k=0; k < P; k++) // profondeur
    for(j=0; j < H; j++) // hauteur
      for(i=0; i < L; i++, aux++){ //largeur
        aux->celltype = (i+k & 1) ? MOINS : PLUS;
      }
#endif

#ifdef MODELE_CRI
  LogPrintf("modele CRISTALLISATION\n");
  // sol plat et arrivee du gaz
  for(k=0; k < P; k++) // profondeur
    for(j=0; j < H; j++) // hauteur
      for(i=0; i < L; i++, aux++){ //largeur
        int H1 = H*0.4;
        int H2 = H*0.8;
        if ((j==H1) && (i<1) && (k<1))
          aux->celltype = IN;
        else
          aux->celltype = (j<H2) ? AIR : PLUS;
      }
#endif
}


float* load_surface(char *nom, int lx, int ly)
{
  FILE *f;
  float *altimap, *pt_al;
  int i, j;

  LogPrintf("lecture surface %s\n", nom);

  altimap = malloc(lx*ly*sizeof(float));

  //read ascii file
  f = fopen(nom,"r");
  if (!f){
    ErrPrintf("ERROR : can't open file %s\n", nom);
    exit(-1);
  }
  pt_al = altimap;
  for(j=0; j < ly; j++)
    for(i=0; i < lx; i++, pt_al++)
      fscanf(f, "%f \n", pt_al);

  fclose(f);

  return altimap;
}


void dump_terre(int parx, int pary)
{
  FILE *fp;
  char str[256], *filename;
  Cell *pt;
  int i, j;
  char *buf, *aux, *ext;
  int dump_type;
  
  
  if ((parx == 1) && (pary == 1)){
    // check extension
    if (csp_filename){
      filename = csp_filename;
      dump_type = DUMP_CSP;
      ext = filename + strlen(filename) - 4;
      if (strcmp(ext,".csp") && strcmp(ext,".CSP")){
        ErrPrintf("ERROR: invalid extension %s\n", filename);
        exit(-1);
      }
    }
    else if (bin_filename){
      filename = bin_filename;
      dump_type = DUMP_BIN;
      ext = filename + strlen(filename) - 4;
      if (strcmp(ext,".bin") && strcmp(ext,".BIN")){
        ErrPrintf("ERROR: invalid extension %s\n", filename);
        exit(-1);
      }
    }
    else{
      ext = opt_bin ? "bin" : "csp";
      dump_type = opt_bin ? DUMP_BIN : DUMP_CSP;
      sprintf(str,"%s_%d_%d_%d.%s",MOD_NAME,H,L,P,ext);
      filename = str;
    }
    LogPrintf("format %s\n", ext+1);
    csp_set_bounds(0);
    write_csp(dump_type, filename);
  }
  else{ //mode parallele (raw binary format) //TODO: csp format
    int LX = L / parx;
    int LY = P / pary;
    int tx, ty, tn, k;
    LogPrintf("LX = %d\nLY = %d\n", LX, LY);
    buf = (char*)malloc(sizeof(char)*LX*H);
    for(tn=0, ty=0; ty<pary; ty++){
      for(tx=0; tx<parx; tx++, tn++){
        sprintf(str, "%s_%d_%d_%d_%d.bin", MOD_NAME, H, LX, LY, tn);
        filename = str;
        fp = fopen(filename,"w");
        if (!fp){
          ErrPrintf("erreur creation fichier %s\n", filename);
          exit(-4);
        }
        LogPrintf("creation fichier %s\n", filename);

        for(k=0; k<LY; k++){
          aux = buf;
          pt = TE + (ty*LY+k)*HL + tx*LX;
          for(j=0; j<H; j++, pt+=L){
            for(i=0; i<LX; i++){
              *aux++ = (pt+i)->celltype;
            }
          }
          fwrite(buf, sizeof(char), LX*H, fp);
        }
        fclose(fp);
        //LogPrintf("Data size: %d\n", H*LX*LY);
      }
    }
    free(buf);
  }
}

void dump_par_conf(int parx, int pary)
{
  FILE *fp;
  int i, j, pid;
  char nom[64];

  pid = 0;

  for (j=0; j<pary; j++)
    for (i=0; i<parx; i++, pid++){
      sprintf(nom, "parallel%d.cfg", pid);
      fp = fopen(nom, "w");
      if ( ! fp ){
        ErrPrintf("erreur ouverture fichier parallel.cfg\n");
        exit(-4);
      }
      LogPrintf("creation fichier de configuration %s\n", nom);
      fprintf(fp, "PID %d\n", pid);
      if (!pid) fprintf(fp, "NN %d\n", parx*pary);
      //if (parx > 1){
        if ((i<parx-1) || (boundary == BC_PERIODIC)) fprintf(fp, "EST %d\n", (i<parx-1)? pid+1 : pid+1-parx);
        if ((i>0) || (boundary == BC_PERIODIC)) fprintf(fp, "OUEST %d\n", (i>0)? pid-1 : pid-1+parx);
      //}
      //if (pary > 1){
        if ((j<pary-1) || (boundary == BC_PERIODIC)) fprintf(fp, "SUD %d\n", (j<pary-1)? pid+parx : i);
        if ((j>0) || (boundary == BC_PERIODIC)) fprintf(fp, "NORD %d\n", (j>0)? pid-parx : i+parx*(pary-1));
      //}
      fclose(fp);
    }

}


void params_genesis()
{
  parametre("Csp_template", "template for the generation of the initial cellular space", &csp_template_str, PARAM_STRING);
#ifdef CYCLAGE_HOR
  parametre("Boundary", "boundary conditions (PERIODIC|OPEN|REINJECTION), PERIODIC by default", &boundary_str, PARAM_STRING);
#else
  parametre("Boundary", "boundary conditions (OPEN|REINJECTION), REINJECTION by default", &boundary_str, PARAM_STRING);
#endif
}

void usage()
{
  printf("genesis - %s model\n",MOD_NAME);

  printf("usage: \n genesis H L D [OPTIONS]\n");
  printf(" genesis -f PARAMETERS_FILE [OPTIONS]\n\n");

  param_usage();

  printf("\nCSP TEMPLATES\n");
  for(int i=0; i<nb_templates; i++) printf("  %s\n", t_templates[i].desc);
  
  printf("\nOPTIONS\n");
  printf("  -bin   \traw binary output\n");
  printf("  -g <n> \trandom seed\n");
#ifdef PARALLEL
  printf("  -par <n>x<m> \tparallel mode, grid <n>x<m>\n");
//  printf("  -node <n> \tnode number <n> in parallel mode\n");
#endif
  
  exit(-1);
}



int main(int argc, char **argv)
{
  int graine = 68374;
  int n=1, parx = 1, pary = 1;
  unsigned char opt_par = 0;

#ifdef LOG_FILE
  log_file = fopen("GENESIS.log","w");
#endif

  init_liste_params();
  init_template_list();
  params_genesis();
    
  if (argc < 3) usage();
    
  LogPrintf("genesis\n");
  
  if (!strcmp(argv[n],"-f")){
    char *ext;
    // read parameters
    lecture_parametres(argv[++n]);
    //assert(filename && strlen(bin_filename)>=4);
    genesis_parse();
  }
  else{
    if (argc < 4) usage();
    H = atoi(argv[n]);
    L = atoi(argv[++n]);
    P = atoi(argv[++n]);
  }
  
  for(n++; n<argc; n++){
    if (!strcmp(argv[n],"-g")){
      graine = atoi(argv[++n]);
    }
    else if (!strcmp(argv[n],"-bin")){
      opt_bin = 1;
    }
#ifdef PARALLEL
    else if (!strcmp(argv[n],"-par")){
      sscanf(argv[++n], "%dx%d", &parx, &pary);
      LogPrintf("geom = %s parx = %d \t pary = %d\n", argv[n], parx, pary);
      opt_par = 1;
      L *= parx;
      P *= pary;
    }
/*
    else if (!strcmp(argv[n],"-node")){
      node = atoi(argv[++n]);
    }
*/
#endif
    else{
      ErrPrintf("ERROR: unknown option %s\n", argv[n]);
    }
  }

  srand48(graine);

  if (opt_par) dump_par_conf(parx, pary);

  genesis();

  dump_terre(parx, pary);

  LogPrintf("genesis : done\n");

  return 1;
}



