/* ReSCAL - 3D space processes
 *
 * Copyright (C) 2011
 *
 * 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.
 */


#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <math.h>
#include <assert.h>
#include <pthread.h>
#include "defs.h"
#include "macros.h"
#include "format.h"
#include "space.h"
#include "doublets.h"
#include "cells.h"
#include "surface.h"
#ifdef LGCA
#include "lgca.h"
#endif
#ifdef PARALLEL
#include "synchro.h"
#endif
#include "transitions.h"
#include "view.h"
#include "simul.h"
#include "trace.h"

extern double csp_time;                  // temps reel simule
extern int pbc_mode;  //periodic boundary conditions

#ifdef PARALLEL
extern int mode_par;    //mode parallele
extern int proc_id;     //numero de process en mode parallele (0 = serveur)
extern int Tunnel[];    //bords partages avec d'autres process
#endif

extern int ava_norm;

#ifdef SURRECTIONS
extern int sur_mode;                // type de surrection (uniforme ou locale)
#endif

int H=0, L=0, P=0, HL=0, HLP=0;       // les dimensions de la terre
int LN, LS, LEO, LNS, HLN;    //couloir est-ouest (limite nord, limite sud, largeur nord-sud, ...)
int Hd2, Ld2, Pd2;
Cell  *TE=NULL;	           // la 'terre'
char *bin_filename=NULL; //nom du fichier BIN
char *csp_filename=NULL; //nom du fichier CSP

#ifdef REFDB_PTR
RefDoublets *RefDB=NULL;       // references des cellules de la terre vers les doublets actifs
#else
RefDoublets_Type *RefDB_Type=NULL;       // references des cellules de la terre vers les doublets actifs
RefDoublets_Ind *RefDB_Ind=NULL;       // references des cellules de la terre vers les doublets actifs
#endif
int graine = 273676;           // graine pour la generation des nombres aleatoires
int h_ceil=0;                 //epaisseur du plafond
int h_floor=0;                 //epaisseur du sol
int nb_pv=0;                  // nombre de plans verticaux (non-DUM)
int *psol;                    // indicateur des plans solides verticaux est-ouest
float csp_angle = 0.0;        //angle resultant de toutes les rotations

pthread_mutex_t mutex_csp;
//pthread_cond_t cond_csp;

void init_terre()
{
  /* Initialize mutex and condition variable objects */
  pthread_mutex_init(&mutex_csp, NULL);
  //pthread_cond_init (&cond_csp, NULL);
}

void lock_csp(int log_flag)
{
#ifdef CSP_MUTEX
  if (log_flag) LogPrintf("lock csp\n");
  pthread_mutex_lock(&mutex_csp);
  if (log_flag) LogPrintf("csp locked\n");
#endif
}

void unlock_csp(int log_flag)
{
#ifdef CSP_MUTEX
  if (log_flag) LogPrintf("unlock csp\n");
  pthread_mutex_unlock(&mutex_csp);
  if (log_flag) LogPrintf("csp unlocked\n");
#endif
}

void cree_terre()
{
  int i, j, k;
  Cell *aux;
  char *filename;

  if (pbc_mode) LogPrintf("cyclage horizontal\n");

#ifdef CELL_TIME
  LogPrintf("datation des cellules (mode 'sedimento')\n");
#endif

  lock_csp(0);

#ifdef PARALLEL
  //bin file
  csp_filename = NULL;
  AllocMemory(bin_filename, char*, 100);
  sprintf(bin_filename,"%s_%d_%d_%d_%d.bin", MOD_NAME, H, L, P, proc_id);
#endif

  filename = csp_filename ? csp_filename : bin_filename;
  if (csp_filename) read_csp_header(filename);

  LNS = P;
  LEO = L;
  H += 2;
  L += 2;  // pour les bords !
  if (!P) P = L; else P += 2;
  HL = H*L;
  HLP = HL*P;
  csp_set_bounds(1);

  AllocMemoryPrint("TE", TE, Cell, HLP);
  ResetMemory(TE, Cell, HLP);

#ifdef REFDB_PTR
  AllocMemoryPrint("RefDB", RefDB, RefDoublets, HLP);
  //ResetMemory(RefDB, RefDoublets, HLP);
#else
  AllocMemoryPrint("RefDB_Type", RefDB_Type, RefDoublets_Type, HLP);
  //ResetMemory(RefDB_Type, RefDoublets_Type, HLP);
  AllocMemoryPrint("RefDB_Ind", RefDB_Ind, RefDoublets_Ind, HLP);
  //ResetMemory(RefDB_Ind, RefDoublets_Ind, HLP);
#endif
  PrintTotalMemory();

#ifdef PARALLEL
  if (mode_par) graine += proc_id;
#endif
  srand48(graine);
  LogPrintf ("graine = %d\n", graine);

  // le centre de la terre
  Hd2 = (int) H/2;
  Ld2 = (int) L/2;
  Pd2 = (int) P/2;

  aux = TE;

#if defined(MODELE_CMB)
  for(k=0; k < P; k++){ // profondeur
    for(j=0; j < H; j++){ // hauteur
      for(i=0; i < L; i++, aux++){ //largeur
        if ((i == 0) || (j == 0) || (k == 0) || (i == L-1) || (j == H-1) || (k == P-1)){
          aux->celltype = BORD;
        }
        else 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);
        }
      }
    }
  }
#elif defined(MODELE_ICB)
  for(k=0; k < P; k++){ // profondeur
    for(j=0; j < H; j++){ // hauteur
      for(i=0; i < L; i++, aux++){ //largeur
        if ((i == 0) || (j == 0) || (k == 0) || (i == L-1) || (j == H-1) || (k == P-1))
          aux->celltype = BORD;
        else if (j < 0.8*H)
          aux->celltype   = ZERO;
        else
          aux->celltype   = PLUS;
         //aux->celltype   = (int) floor(drand48()*3);
      }
    }
  }
#elif defined(MODELE_DUN)
  for(k=0; k < P; k++){ // profondeur
    for(j=0; j < H; j++){ // hauteur
      for(i=0; i < L; i++, aux++){ //largeur
        if ((i == 0) || (j == 0) || (k == 0) || (i == L-1) || (j == H-1) || (k == P-1))
          aux->celltype = BORD;
        else if ((j == 1) || (j == H-2))
          aux->celltype   = DUM;
        else if (j > 0.9*H)
          aux->celltype   = GR;
        else
          aux->celltype   = EAUC;
         //aux->celltype   = (int) floor(drand48()*3);
      }
    }
  }
#else
  for(k=0; k < P; k++){ // profondeur
    for(j=0; j < H; j++){ // hauteur
      for(i=0; i < L; i++, aux++){ //largeur
        if ((i == 0) || (j == 0) || (k == 0) || (i == L-1) || (j == H-1) || (k == P-1))
          aux->celltype = BORD;
        else
          aux->celltype   = (int) floor(drand48()*2);
      }
    }
  }
#endif

#ifdef REFDB_PTR
  memset(RefDB, 0, sizeof(RefDoublets) * HLP);
#else
  memset(RefDB_Type, -1, sizeof(RefDoublets_Type) * HLP);
  memset(RefDB_Ind, -1, sizeof(RefDB_Ind) * HLP);
#endif

  if (filename) read_csp(filename);

  //signal_csp_ready();
  unlock_csp(0);

#ifdef LISSAGE
  init_lissage();
#endif
#ifdef CELL_COLOR
  //init_color();
#endif
}

#ifdef PARALLEL
void cree_tunnels()
{
  int i, j, k;
  Cell *aux;

  aux = TE;

  for(k=0; k < P; k++) // profondeur
    for(j=0; j < H; j++) // hauteur
      for(i=0; i < L; i++, aux++) //largeur
        if (((Tunnel[EST] >= 0) && (i == L-1)) || ((Tunnel[OUEST] >= 0) && (i == 0))
        || ((Tunnel[BAS] >= 0) && (j == H-1)) || ((Tunnel[HAUT] >= 0) && (j == 0))
        || ((Tunnel[SUD] >= 0) && (k == P-1)) || ((Tunnel[NORD] >= 0) && (k == 0)))
          aux->celltype = TUNNEL;

}
#endif

void apocalypse()
{
  LogPrintf("destruction terre !\n");
#ifdef REFDB_PTR
  if (RefDB) free(RefDB);
#else
  if (RefDB_Type) free(RefDB_Type);
  if (RefDB_Ind) free(RefDB_Ind);
#endif
  if (TE) free(TE);
}


void cellule_terre(unsigned char type, int ix)
{
  //changer l'etat d'une cellule de la terre
  elimine_doublet_est(ix);
  elimine_doublet_ouest(ix);
  elimine_doublet_bas(ix);
  elimine_doublet_haut(ix);
  elimine_doublet_sud(ix);
  elimine_doublet_nord(ix);

  modifie_cellule(type, ix);

  ajoute_doublet_est(ix);
  ajoute_doublet_ouest(ix);
  ajoute_doublet_bas(ix);
  ajoute_doublet_haut(ix);
  ajoute_doublet_sud(ix);
  ajoute_doublet_nord(ix);
}

int hcycle(int ix)
{
  int x, y, z;
  Calcule_xyz(ix, x, y, z);
  if (x == 0) x = L-2; else if (x == L-1) x = 1;
  if (z == 0) z = P-2; else if (z == P-1) z = 1;
  ix = x + y*L + z*HL;
  return ix;
}

#ifdef CELL_COLOR
float color_coef=0.1; //0.1;         // coefficient associe a la callback check_color

void init_color()
{
  int i, j, k;
  int nbcelcol=0;

  //LogPrintf("initialisation des couleurs %d %d\n", sizeof(Cell), sizeof(Cell2));
  //initialise la couleur des cellules
  for(k=1; k < P-1; k++) // profondeur
    for(j=1; j < H-1; j++) // hauteur
      for(i=1; i < L-1; i++){ //largeur
#ifdef MODELE_DUN
    TE[k*HL+j*L+i].color = 0;
    //TE[k*HL+j*L+i].color = (((TE[k*HL+j*L+i].celltype == GR) || (TE[k*HL+j*L+i].celltype == GRJ)) && (k<L/2) /*(j<H-5)*/ /*(i<L/2)*/) ? 1 : 0;
    //TE[k*HL+j*L+i].color = (((TE[k*HL+j*L+i].celltype == GR) || (TE[k*HL+j*L+i].celltype == GRJ)) && (drand48()<0.5)) ? 1 : 0;
    //TE[k*HL+j*L+i].color = ((TE[k*HL+j*L+i].celltype == GR) || (TE[k*HL+j*L+i].celltype == GRJ)) ? 1 : 0;
    //TE[k*HL+j*L+i].color = (TE[k*HL+j*L+i].celltype == IN) ? 1 : 0;
#endif
#ifdef MODELE_AVA
    TE[k*HL+j*L+i].color = (TE[k*HL+j*L+i].celltype == IN) ? 1 : 0;
#endif
    if (TE[k*HL+j*L+i].color) nbcelcol++;
  }
  LogPrintf("initialisation des couleurs : %d cellule(s) coloree(s)\n", nbcelcol);
}

// callback de controle en fonction de la couleur
int check_color(int ix, void *data)
{
  static int start = 1;
  //static float coef = 0.1;
  DataCheck *pdc = data;
  int chk_ok = 1;

  if (start){
    LogPrintf("controle en fonction de la couleur des grains - coef = %f - col = %d\n", pdc->coef, pdc->col);
    start = 0;
  }

  if (TE[ix].color == pdc->col){
    //calcul probabiliste
    //chk_ok = (drand48() < color_coef);
    chk_ok = (drand48() < pdc->coef);
  }
  //if (!chk_ok) {LogPrintf(" check_color : not ok (%d) - coef = %f\n", pdc->col, pdc->coef);}
  return chk_ok;
}
#endif

#if CELL_DATA
//char cdata_offset = sizeof(TE->celltype);
/*
//initialisation des donnees de la cellule ix
void init_cell_data(CellData data, int ix)
{
  memcpy((void *)(TE+ix) + CELL_TYPE_SIZE, (void *)(&data), CELL_DATA_SIZE);
}

//copier les donnees de la cellule ix dans data
void get_cell_data(CellData data, int ix)
{
  memcpy((void *)(&data), (void *)(TE+ix) + CELL_TYPE_SIZE, CELL_DATA_SIZE);
}
*/

//permuter les donnees de 2 cellules
void swap_cell_data(int ix, int ix2)
{
  Cell aux;
  char flag1 = 1;
  char flag2 = 1;

  //LogPrintf("offset = %d\n", offset);

  aux = TE[ix];
  //LogPrintf("%08x %08x\n", (&TE[ix]), (unsigned char*)(TE+ix)+offset)

  // les cellules BORD ne doivent pas etre modifiees
  if (TE[ix].celltype == BORD) flag1 = 0;
  if (TE[ix2].celltype == BORD) flag2 = 0;

#ifdef IN
  // les cellules IN ne doivent pas etre modifiees
  if (TE[ix].celltype == IN) flag1 = 0;
  if (TE[ix2].celltype == IN) flag2 = 0;
#endif

#ifdef OUT
  // les cellules OUT ne doivent pas etre modifiees
  if (TE[ix].celltype == OUT) flag1 = 0;
  if (TE[ix2].celltype == OUT) flag2 = 0;
#endif

  if (flag1) memcpy((void *)(TE+ix) + CELL_TYPE_SIZE, (void *)(TE+ix2) + CELL_TYPE_SIZE, CELL_DATA_SIZE); //TE[ix].celldata <- TE[ix2].celldata
  if (flag2) memcpy((void *)(TE+ix2) + CELL_TYPE_SIZE, (void *)(&aux) + CELL_TYPE_SIZE, CELL_DATA_SIZE); //TE[ix2].celldata <- aux.celldata
}

void update_inout_data(int ix, int ix2)
{
#ifdef CELL_COLOR
#ifdef IN
//#if 0
  if (TE[ix].celltype == IN) TE[ix].color = (drand48()<0.5) ? 1 : 0;
  if (TE[ix2].celltype == IN) TE[ix2].color = (drand48()<0.5) ? 1 : 0;
#endif
#endif
}

//copier les donnees de la cellule ix dans la cellule ix2
void copy_cell_data(int ix, int ix2)
{
  memcpy((void *)(TE+ix2) + CELL_TYPE_SIZE, (void *)(TE+ix) + CELL_TYPE_SIZE, CELL_DATA_SIZE); //TE[ix2].celldata <- aux.celldata
}

#endif  //CELL_DATA

// callback de controle en fonction de l'altitude
int check_alti(int ix, void *data)
{
  static int start = 1;
  if (start){
    LogPrintf("controle du flux entrant en fonction de l'altitude\n");
    start = 0;
  }
  //calcul de la position (x,y,z)
  int x, y, z;
  Calcule_xyz(ix, x, y, z);

  //calcul probabiliste suivant y
  float alea = drand48();
  float var = 1.0-(float)y/H; //probabilite plus forte si altitude elevee (y faible)
  //float var = 1.0-2.0*(float)y/H;
  int chk_ok = (alea < var);
  //if (!chk_ok) {LogPrintf(" check_alti : not ok\n");}
  return chk_ok;
}

#ifdef MODELE_DUN
// callback de controle en fonction de la cellule a l'ouest
int check_grain_seul(int index, void *data)
{
  return TE[index-1].celltype == EAUC;
}
#endif

#ifdef MODELE_LIFE
//threshold values for birth transitions
int B_min = 3;
int B_max = 3;
//threshold values for cells to stay alive
int S_min = 2;
int S_max = 3;

/// counts living neighbours
int count_neighbours(int ix)
{
  int nb = 0;

  if (TE[ix+1].celltype == ALIVE) nb++;
  if (TE[ix-1].celltype == ALIVE) nb++;
  if (TE[ix+HL].celltype == ALIVE) nb++;
  if (TE[ix-HL].celltype == ALIVE) nb++;
  if (TE[ix+HL+1].celltype == ALIVE) nb++;
  if (TE[ix+HL-1].celltype == ALIVE) nb++;
  if (TE[ix-HL+1].celltype == ALIVE) nb++;
  if (TE[ix-HL-1].celltype == ALIVE) nb++;

#ifdef CYCLAGE_HOR
  if (pbc_mode){
    if (TE[ix].celltype == BORD) {ErrPrintf("ERROR in count_neighbours: BORD cell\n"); exit(-1);}

    int ix2;
    ix2 = ix+1;
    if (TE[ix2].celltype == BORD){
      ix2 = hcycle(ix2);
      if (TE[ix2].celltype == ALIVE) nb++;
    }
    ix2 = ix-1;
    if (TE[ix2].celltype == BORD){
      ix2 = hcycle(ix2);
      if (TE[ix2].celltype == ALIVE) nb++;
    }
    ix2 = ix+HL;
    if (TE[ix2].celltype == BORD){
      ix2 = hcycle(ix2);
      if (TE[ix2].celltype == ALIVE) nb++;
    }
    ix2 = ix-HL;
    if (TE[ix2].celltype == BORD){
      ix2 = hcycle(ix2);
      if (TE[ix2].celltype == ALIVE) nb++;
    }
    ix2 = ix+HL+1;
    if (TE[ix2].celltype == BORD){
      ix2 = hcycle(ix2);
      if (TE[ix2].celltype == ALIVE) nb++;
    }
    ix2 = ix+HL-1;
    if (TE[ix2].celltype == BORD){
      ix2 = hcycle(ix2);
      if (TE[ix2].celltype == ALIVE) nb++;
    }
    ix2 = ix-HL+1;
    if (TE[ix2].celltype == BORD){
      ix2 = hcycle(ix2);
      if (TE[ix2].celltype == ALIVE) nb++;
    }
    ix2 = ix-HL-1;
    if (TE[ix2].celltype == BORD){
      ix2 = hcycle(ix2);
      if (TE[ix2].celltype == ALIVE) nb++;
    }
  }
#endif

  return nb;
}

/// check callback for the birth transition of a cell
int check_birth(int ix, void *data)
{
  int nb = count_neighbours(ix);

  return (nb >= B_min) && (nb <= B_max);
}

/// check callback for the death transition of a cell
int check_death(int ix, void *data)
{
  int nb = count_neighbours(ix);

  return (nb < S_min) || (nb > S_max);
}
#endif

#ifdef SURRECTIONS

unsigned char *surrection_map = NULL;

void init_surrection()
{
  int i, k, x1, x2, z1, z2;
  float alea;
  int Lmap = L-2;
  int Pmap = P-2;
  unsigned char *aux;
  Cell *pt;

  if (sur_mode == SUR_MODE_UNIFORM) {LogPrintf("surrection uniforme\n");}  // \t iter = %lu%09lu \n", md_iter, iter);
  else if (sur_mode == SUR_MODE_CONE) {LogPrintf("surrection conique\n");}  //\t iter = %lu%09lu \t", md_iter, iter);
  else if (sur_mode == SUR_MODE_TECTO) {LogPrintf("surrection tectonique\n");} //\t iter = %lu%09lu \t", md_iter, iter);

  if ((sur_mode == SUR_MODE_UNIFORM) && (surrection_map != NULL)) return;  //surrection deja initialisee

  if (!surrection_map){
    AllocMemory(surrection_map, unsigned char, Lmap*Pmap);
  }

  ResetMemory(surrection_map, unsigned char, Lmap*Pmap);

  x1 = 1;
  z1 = 1;
  x2 = Lmap-1;
  z2 = Pmap-1;

  if (sur_mode == SUR_MODE_UNIFORM){
    for(k=z1; k < z2; k++){ // profondeur
      aux = surrection_map + x1 + k*Lmap;
      pt = TE+1+x1+(H-2)*L+k*HL;
      for(i=x1; i < x2; i++, aux++){ //largeur
        *aux = ((pt++)->celltype != DUM);
      }
    }
  }

  if (sur_mode == SUR_MODE_CONE){
    int di, dk, r2;
    alea = drand48();
    r2 = alea*alea*Lmap*Pmap/8;
    LogPrintf("rayon = %d\n", (int)sqrt(r2));
    for(k=z1; k < z2; k++){ // profondeur
      aux = surrection_map + x1 + k*Lmap;
      for(i=x1; i < x2; i++, aux++){ //largeur
        di = i - Lmap/2;
        dk = k - Pmap/2;
        *aux = (di*di + dk*dk < r2);
      }
    }
  }

  if (sur_mode == SUR_MODE_TECTO){
    int prof;
    //alea = drand48();
    //prof = alea*Pmap/2;
    prof = Pmap;
    LogPrintf("profondeur = %d\n", prof);
    for(k=z1; k < z2; k++){ // profondeur
      aux = surrection_map + x1 + k*Lmap;
      for(i=x1; i < x2; i++, aux++){ //largeur
        int dk = Pmap*0.45 - k;
        *aux = ((dk >= 0) && (dk < prof));
      }
    }
  }
}

int surrection()
{
  int i, j, k, n;
  Cell *aux;
  unsigned char *auxmap;

  //initialisation du tableau surrection_map
  init_surrection();

  auxmap = surrection_map;
  n = 0;
  for(k=1; k < P-1; k++){ // profondeur
    for(i=1; i < L-1; i++, auxmap++){ //largeur
      if (*auxmap){
        aux = TE + i + L + k*HL;
        for(j=1; j < H-2; j++, aux+=L){ // hauteur
          if (aux->celltype != IN) *aux   = *(aux+L);
        }
        aux->celltype = ALTI;
#ifdef CELL_TIME
        aux->celltime = 0;
#endif
        n++;
      }
    }
  }

  //comptage des cellules
  init_Ncel();

  //unlock csp condition
  //signal_csp_ready();

  //recalcul des tableaux de doublets actifs
  init_db_pos();

#ifdef ALTI
  calcule_alti(ALTI, ALTI_MODE_BAS);
  if (ava_norm) calcule_normales();
#endif

#ifdef LGCA
  //mise-a-jour de la matrice de gaz-sur-reseau
  collisions_mod_terre();
#endif

  LogPrintf("surrection: %d cellule(s)\n", n);

  return n;
}

#endif //SURRECTIONS


void translation(int dx, int dz)
{
  static Cell *auxmap = NULL;
  Cell *aux, *aux2;
  int i, j, k, i0, k0;

  if (!auxmap) AllocMemoryPrint("auxmap (translation)", auxmap, Cell, HLP);

  LogPrintf("translation : ( %d , %d )\n", dx, dz);

  aux = TE;
  aux2 = auxmap;

  for(k=0; k < P; k++){ // profondeur
    for(j=0; j < H; j++){ // hauteur
      for(i=0; i < L; i++, aux++, aux2++){ //largeur
        if ((aux->celltype != BORD) && (aux->celltype != DUM) && (aux->celltype != IN) && (aux->celltype != OUT)){
          i0 = i-dx;
          k0 = k-dz;
          if (pbc_mode){
            while (i0<1) i0+=L-2;
            while (i0>=L-1) i0-=L-2;
            while (k0<1) k0+=P-2;
            while (k0>=L-1) k0-=P-2;
            *aux2 = *(TE+i0+j*L+k0*HL);
          }
          else{
            if ((i0>=1) && (i0<L-1) && (k0>=1) && (k0<P-1))
              *aux2 = *(TE+i0+j*L+k0*HL);
#ifdef MODELE_DUN
            else
              aux2->celltype = EAUC;
#endif
          }
        }
        else
          *aux2   = *aux;
      }
    }
  }
  //LogPrintf("fin translation\n");
  memcpy(TE, auxmap, HLP*sizeof(Cell));

  //recalcul du nombre de cellules de chaque etat
  if (!pbc_mode) init_Ncel();

  //unlock csp condition
  //signal_csp_ready();

  //recalcul des tableaux de doublets actifs
  init_db_pos();

#ifdef ALTI
  calcule_alti(ALTI, ALTI_MODE_BAS);
  if (ava_norm) calcule_normales();
#endif

#ifdef LGCA
  //mise-a-jour de la matrice de gaz-sur-reseau
  collisions_mod_terre();
#endif
}

void rotation(float angle, char flags)
{
  static Cell *csp_tmp = NULL;
  static Cell *TE0 = NULL;
  Cell *aux, *aux2;
  float alpha;
  float co;
  float si;
  float di0, dk0, R;
  float di, dk;
  int i0, k0;
  int i, j, k;

  if (!csp_tmp) AllocMemoryPrint("csp_tmp (rotation)", csp_tmp, Cell, HLP);
  //ResetMemory(auxmap, Cell, HLP);

  if (flags & ROT_REORIENT_TEMP){
    if (!csp_angle) return;
    //LogPrintf("reorientation terre\n");
    angle = -csp_angle;
  }
  else if (flags & ROT_REORIENT_UNDO){
    if (!TE0) return;
    //LogPrintf("restauration terre\n");
    angle = csp_angle;
  }
  else{
    if (!angle) return;
    csp_angle = fmodf(csp_angle + angle, 360.0);
    LogPrintf("rotation angle: %f - orientation: %f - time: %f\n", angle, csp_angle, csp_time);
    //mode0 = rot_mode;
  }

  //co = 0.0;
  //si = 1.0;
  //LogPrintf("co = %.f   si = %.f\n", co, si);

  aux = TE;
  aux2 = csp_tmp;

  di0 = Ld2 - 0.5;
  dk0 = Pd2 - 0.5;
  R = (P<L) ? Pd2-10 : Ld2-10;
  alpha = angle*PI/180.0;
  co = cos(alpha);
  si = sin(alpha);

  // rotation sur table tournante
  if (flags & ROT_DISK){
    //int r2 = (Ld2-2)*(Ld2-2);
    //float r2 = (Ld2-12)*(Ld2-12); //en presence d'un trou cylindrique, prevoir une table tournante de rayon inferieur d'une cellule
    float r2 = R*R; //en presence d'un trou cylindrique, prevoir une table tournante de rayon inferieur d'une cellule
    for(k=0; k < P; k++){ // profondeur
      for(j=0; j < H; j++){ // hauteur
        for(i=0; i < L; i++, aux++, aux2++){ //largeur
          di = i - di0;
          dk = k - dk0;
          //if ((di*di + dk*dk < r2) && (aux->celltype != DUM) && (aux->celltype != IN) && (aux->celltype != OUT)){
          if (di*di + dk*dk < r2){  //modif oli
            i0 = roundf(di0 + co*di + si*dk);
            k0 = roundf(dk0 - si*di + co*dk);
            //aux2->celltype = (TE+i0+j*L+k0*HL)->celltype;
            *aux2 = TE[i0+j*L+k0*HL];
#ifdef MODELE_DUN
              //if (aux2->celltype == DUM) aux2->celltype = EAUC; //modif oli
#endif
          }
          else
            //aux2->celltype = aux->celltype;
            *aux2 = *aux;
        }
      }
    }
  }

  // rotation avec cyclage horizontal
  if (flags & ROT_CYCL){
    for(k=0; k < P; k++){ // profondeur
      for(j=0; j < H; j++){ // hauteur
        for(i=0; i < L; i++, aux++, aux2++){ //largeur
          di = i - di0;
          dk = k - dk0;
          if ((aux->celltype != BORD) && (aux->celltype != DUM) && (aux->celltype != IN)){
            i0 = roundf(di0 + co*di + si*dk);
            k0 = roundf(dk0 - si*di + co*dk);
#ifdef TRACE_PAR_COL
            char no_trace = 0;
            if ((i0<1) || (i0>=L-1) || (k0<1) || (k0>=P-1)) no_trace=1;
#endif
            while (i0<1) i0+=L-2;
            while (i0>=L-1) i0-=L-2;
            while (k0<1) k0+=P-2;
            while (k0>=P-1) k0-=P-2;
            //aux2->celltype = (TE+i0+j*L+k0*HL)->celltype;
            *aux2 = TE[i0+j*L+k0*HL];
#ifdef TRACE_PAR_COL
            if (no_trace) aux2->color = 0;
#endif

#ifdef MODELE_DUN
            if (aux2->celltype == DUM) aux2->celltype = EAUC;
#endif
          }
          else
            //aux2->celltype = aux->celltype;
            *aux2 = *aux;
        }
      }
    }
  }

  if (flags & ROT_REORIENT_TEMP){
    //save address of current cellular space
    TE0 = TE;
    //switch cellular space
    TE = csp_tmp;

#ifdef TRACE_PAR_COL
    //correction of the colors
    trace_par_col_rotation(0.0);
#endif
  }
  else if (flags & ROT_REORIENT_UNDO){
    //restore cellular space
    TE = TE0;
    TE0 = NULL;
  }
  else{
    //LogPrintf("fin rotation\n");
    memcpy(TE, csp_tmp, HLP*sizeof(Cell));

    //recalcul du nombre de cellules de chaque etat
    init_Ncel();

    //unlock csp condition
    //signal_csp_ready();

    //recalcul des tableaux de doublets actifs
    init_db_pos();

#ifdef LGCA
    //mise-a-jour de la matrice de gaz-sur-reseau
    collisions_mod_terre();
#endif

#ifdef TRACE_PAR_COL
    //correction of the colors
    trace_par_col_rotation(csp_angle);
#endif
  }

#ifdef ALTI
  calcule_alti(ALTI, ALTI_MODE_BAS);
  if (ava_norm) calcule_normales();
#ifdef ROTATE_LIGHT
  if (flags == ROT_CYCL) rotate_light(angle);
#endif
#endif
}

void rotation90(int n)
{
  static Cell *auxmap = NULL;
  Cell *aux, *aux2;
  int i, j, k;
  int i0, k0;
  float a = n*PI/2;
  float co = cos(a);
  float si = sin(a);
  float di0 = Ld2 - 0.5;
  float dk0 = Ld2 - 0.5;
  float di, dk;

  if (!auxmap) AllocMemoryPrint("auxmap (rotation90)", auxmap, Cell, HLP);
  ResetMemory(auxmap, Cell, HLP);

  csp_angle = fmodf(csp_angle + n*90, 360.0);
  LogPrintf("rotation90 angle : %d - orientation : %f\n", n*90, csp_angle);

  aux = TE;
  aux2 = auxmap;
  for(k=0; k < P; k++){ // profondeur
    for(j=0; j < H; j++){ // hauteur
      for(i=0; i < L; i++, aux++, aux2++){ //largeur
        if ((aux->celltype != BORD) && (aux->celltype != DUM) && (aux->celltype != IN) && (aux->celltype != OUT)){
          di = i - di0;
          dk = k - dk0;
          i0 = roundf(di0 + co*di + si*dk);
          k0 = roundf(dk0 - si*di + co*dk);
          //i0 = di0 + co*(float)di - si*(float)dk;
          //k0 = dk0 + si*(float)di + co*(float)dk;
          while (i0<1) i0+=L-2;
          while (i0>=L-1) i0-=L-2;
          while (k0<1) k0+=P-2;
          while (k0>=P-1) k0-=P-2;
          *aux2 = TE[i0+j*L+k0*HL];
        }
        else
          *aux2 = *aux;
      }
    }
  }
  //LogPrintf("fin rotation\n");
  memcpy(TE, auxmap, HLP*sizeof(Cell));

  //recalcul des tableaux de doublets actifs
  init_db_pos();

  //recalcul du nombre de cellules de chaque etat
  init_Ncel();

#ifdef ALTI
  calcule_alti(ALTI, ALTI_MODE_BAS);
  if (ava_norm) calcule_normales();
#endif

#ifdef LGCA
  //mise-a-jour de la matrice de gaz-sur-reseau
  collisions_mod_terre();
#endif
}


void conditions_bords()
{
  Cell *aux;
  int j, k, k1, k2, ll;

  //ouverture aleatoire d'une partie du bord gauche (sans DUM)
  ll = P/4;
  k1 = drand48()*(P-ll);
  k2 = k1 + ll;
  LogPrintf ("bord ouvert entre %d et %d\n", k1, k2);

  for(k=1; k < P-1; k++) // profondeur
    for(j=1; j < H-1; j++){ // hauteur
      aux = TE+1+j*L+k*HL;
      if ((k >= k1) && (k <= k2) && (j < H-2))
        aux->celltype = (aux+1)->celltype; //ouverture du bord gauche
      else
        aux->celltype = DUM;
    }

  //recalcul des tableaux de doublets actifs
  init_db_pos();
}

#ifdef LISSAGE
void init_lissage()
{
  int i, j, k;
  Cell *aux;

  LogPrintf("init_lissage\n");
#ifdef MODELE_RIV
  // on depose une couche de 5 cellules de boue
  for(k=1; k < P-1; k++){ // profondeur
    for(i=1; i < L-1; i++, aux++){ //largeur
      aux = TE + i + k*HL;
      while (((aux+L)->celltype != TERRE) && ((aux+L)->celltype != BORD)){
        aux += L;
        if (aux->celltype == BOUE) aux->celltype = EAU;
      }
      for (j=0; (j<5) && (aux->celltype == EAU); j++, aux-=L){
        aux->celltype = BOUE;
      }
    }
  }
#endif
}

void lissage_terre()
{
  int i, n;
  Cell *aux;

#ifdef MODELE_RIV
  // on change la boue piegee en terre
  aux = TE;
  n = 0;
  for (i=0; i<HLP; i++, aux++){
    if (aux->celltype == BOUE){
      aux->celltype = TERRE;
      n++;
    }
  }
#endif
  LogPrintf("lissage_terre : %d cellules\n", n);
}
#endif  //LISSAGE

int dummy_h_plan(int y)
{
  int i,k,ix;
  int res=1;

  //test si le plan horizontal y est entierement constitue de cellules DUM|BORD|IN|OUT
  for(k=0; (k<P) && (res); k++){ // profondeur
    ix = y*L + k*HL;
    for(i=0; (i<L) && (res); i++, ix++){ //largeur
      res = (TE[ix].celltype == DUM) || (TE[ix].celltype == BORD) || (TE[ix].celltype == IN) || (TE[ix].celltype == OUT);
    }
  }
  //LogPrintf("plan_h_solide %d : res = %d   ix = %d\n", y, res, ix);

  return res;
}

int dummy_v_plan(int z)
{
  int ix;
  int res=1;

  //test si le plan vertical z est entierement constitue de cellules DUM|BORD
  ix=z*HL;
  while ((ix<(z+1)*HL) && res){
    res = (TE[ix].celltype == DUM) || (TE[ix].celltype == BORD);
    ix++;
  }
  //LogPrintf("plan_v_solide %d : res = %d   ix = %d\n", z, res, ix);

  return res;
}


void calcule_couloir()
{
  int i;

  AllocMemory(psol,int,P);
  //ResetMemory(psol,int,L);

  //initialisation du tableau psol
  //LogPrintf("calcule_couloir: P=%d\n",P);
  for (i=0; i<P; i++){
    psol[i] = dummy_v_plan(i);
    if (!psol[i]) nb_pv++;
    /*if (psol[i]) LogPrintf("plan solide : z = %d\n", i); */
  }

  //limite nord du couloir
  for(i=0; (i<P) && psol[i]; i++);
  LN=i;

  //limite sud du couloir
  //for(; (i<L) && !plan_solide(i); i++);
  for(i=P-1; (i>0) && psol[i]; i--);
  LS=i;

  //largeur du couloir
  LNS=LS-LN+1;
  /*if (LNS==L-2){
    LN=0;
    LS=L-1;
    LNS=L;
  }*/
  HLN=LN*HL;
  //if (LNS>0){
  if (1){
    LogPrintf("couloir : LN = %d   LS = %d   LNS = %d   nb_pv = %d\n", LN, LS, LNS, nb_pv);
    if (LNS<=0){ ErrPrintf("ERROR: corridor not found\n");  exit(-1);}

  }
}

//calcul epaisseur du plafond en DUM
void compute_h_ceil()
{
  int j;

  for(j=1; (j<H) && dummy_h_plan(j); j++); // hauteur

  h_ceil = j-1;

  LogPrintf("h_ceil = %d\n", h_ceil);

  //return j;
}

//calcul epaisseur du sol en DUM
void compute_h_floor()
{
  int j;

  for(j=H-2; (j>0) && dummy_h_plan(j); j--); // hauteur

  h_floor = H-2-j;

  LogPrintf("h_floor = %d\n", h_floor);

  //return j;
}

#ifdef DUMP_RUGOSI

void dump_rugosi(/*double csp_time*/)
{
  static int cpt=0;
	double lmean, hmean, rugH, rugB, rugHmean, rugBmean;
	Cell *haut, *bas;
	int i,j,top,bot;
	static int *Vtop=NULL, *Vbot=NULL;
	char nom[256];
	FILE *fp;

	if ( ! Vtop ) Vtop = (int *) malloc(sizeof(int)*L*P);
	if ( ! Vbot ) Vbot = (int *) malloc(sizeof(int)*L*P);

	lmean = hmean = rugH = rugB = rugHmean = rugBmean = 0.0;

  for(j=1; j<P-1; j++)
    for(i=1; i<L-1; i++){
	  	haut = TE + j*HL + L + i;
		  bas  = TE + j*HL + L*(H-2) + i;

  		top = 0;
#ifdef MODELE_CMB
	  	while( haut->celltype == PLUS ) {
#else
	  	while( haut->celltype != PLUS ) {
#endif
		  	haut += L;
			  top++;
		  }
		  Vtop[i+j*L] = top;
		  rugHmean += top;

		  bot = H-2;
#ifdef MODELE_CMB
		  while( bas->celltype  == MOINS) {
#else
		  while( bas->celltype  == PLUS) {
#endif
			  bas  -= L;
			  bot--;
		  }
		  Vbot[i+j*L] = bot;
		  rugBmean += bot;

		  hmean += (bot - top + 1);
	  }

	rugHmean /= (L-2)*(P-2);
	rugBmean /= (L-2)*(P-2);
  hmean /= (L-2)*(P-2);
	lmean = (rugHmean + rugBmean)/2;

	for(j=1; j<P-1; j++)
	  for(i=1; i<L-1; i++){
		  rugH += (Vtop[i+j*L] - rugHmean)*(Vtop[i+j*L] - rugHmean);
		  rugB += (Vbot[i+j*L] - rugBmean)*(Vbot[i+j*L] - rugBmean);
	  }

	rugH = sqrt( rugH/((L-2)*(P-2)) );
	rugB = sqrt( rugB/((L-2)*(P-2)) );


	sprintf(nom,"RUGOS.log");
	fp = fopen(nom, "ab");
	if( ! fp ){ErrPrintf("ouverture fic. dump rugosite impossible\n");exit(-5);}
	fprintf(fp,"%e %15.6lf %15.6f %15.6lf %15.6lf\n",csp_time,lmean,hmean,rugH,rugB);
	fclose(fp);

  if ( ((cpt++) & 7) == 7 ){ /* 1 fois sur 8 */
    // fichier interface TOP
    //sprintf(nom,"TOP_%d_%d_%04lu%09lu_%e.log",H-2,L-2,md_iter,iter,csp_time);
    sprintf(nom,"TOP_%d_%d.log",L-2,P-2);
    fp = fopen(nom, "w");
  	if( ! fp ){ErrPrintf("ouverture fichier dump TOP impossible\n");exit(-5);}
  	for(j=1; j<P-1; j++){
  	  for(i=1; i<L-1; i++){
        fprintf(fp,"%d ",(H-2)-Vtop[i+j*L]);
      }
      fprintf(fp,"\n");
    }
	  fclose(fp);

    // fichier interface BOT
    //sprintf(nom,"BOT_%d_%d_%04lu%09lu_%e.log",H-2,L-2,md_iter,iter,csp_time);
    sprintf(nom,"BOT_%d_%d.log",L-2,P-2);
    fp = fopen(nom, "w");
  	if( ! fp ){ErrPrintf("ouverture fichier dump BOT impossible\n");exit(-5);}
	  for(j=1; j<P-1; j++){
	    for(i=1; i<L-1; i++){
        fprintf(fp,"%d ",(H-2)-Vbot[i+j*L]);
      }
      fprintf(fp,"\n");
    }
	  fclose(fp);
  }

}
#endif


void dump_terre(char dump_type, int cpt, int unit)
{
  char filename[100];
  char str[100];
  char *ext;

  ext = (dump_type == DUMP_CSP) ? "csp" : "bin";
#ifdef PARALLEL
  extern int proc_id;
  sprintf(str,"_%d",proc_id);
#else
  *str = 0;
#endif
  if (unit == UNIT_T0)
    sprintf(filename,"%s%05d_t0%s.%s",MOD_NAME,cpt,str,ext);
  else
    sprintf(filename,"%s%04d%s.%s",MOD_NAME,cpt,str,ext);

  //sprintf(filename,"TERRE_%d_%d_%04lu%09lu_%e.bin",H-2,L-2,md_iter,iter,csp_time);

  LogPrintf("write CSP: %s, csp_time = %f (t0)\n", filename, csp_time);

  write_csp(dump_type, filename);

#ifdef APPEL_GMT
  if (dump_type == DUMP_BIN){
    char command[100];
    sprintf(command, "./bin2map %s %d %d >> GMT.log", filename, H-2, L-2);
    system(command);
    //sprintf(command, "mv ero.jpeg GMT%04d.jpeg", cpt);
    //system(command);
  }
#endif

#ifdef APPEL_MCC
  if (dump_type == DUMP_BIN){
    char command[100];
    fp = fopen("PARAMS.txt","w");
    fprintf(fp,"datafile='%s'\n",filename);
    fprintf(fp,"h=%d;\n",H-2);
    fprintf(fp,"lx=%d;\n",L-2);
    fprintf(fp,"ly=%d;\n",P-2);
    fprintf(fp,"lim_alt=1;\n");
    fprintf(fp,"seuil_dune=10;\n");
    fprintf(fp,"seuil_dens=0.2;\n");
    fprintf(fp,"ezon=8;\n");
    fprintf(fp,"flag_isole_dune=1;\n");
    fprintf(fp,"show_figs=1;\n");
    fprintf(fp,"save_figs=1;\n");
    fprintf(fp,"save_stats=1;\n");
    fclose(fp);
    //sprintf(command, "cat PARAMS.tail >> PARAMS.txt");
    //system(command);
    sprintf(command, "./contour_dunes >> MAT.log");
    system(command);
  }
#endif

  compress(filename,1);
}


#ifdef DUMP_SIGNATURE
void dump_signature(int ii)
{
  FILE *fp;
  int i;
  unsigned int sig, *aux;

  //calcul de la signature
  sig=0;
  aux = (unsigned int*)TE;
  for(i=0; i<HLP*sizeof(Cell)/sizeof(unsigned int); i++, aux++) sig += (*aux); //sig = sig ^ (*aux);

  //dump
  fp = fopen("SIGN_CSP.log","a");
  if ( ! fp ){
    ErrPrintf("erreur ouverture fichier dump signature CSP\n");
    exit(-4);
  }
  fprintf(fp, "%08x : %08x\n", ii, sig);
  fclose(fp);

}
#endif

