/* ReSCAL - Lattice gas
 *
 * 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 "defs.h"
#include "macros.h"
#include "param.h"
#include "cells.h"
#include "surface.h"
#include "transitions.h"
#include "lgca.h"
#include "trace.h"
#include "simul.h"

//#undef CYCLAGE_HOR	//pas de cyclage sur le flux (optionnel) -> sinon essoufflement du flux (frottements)

#ifdef LGCA

//#ifdef MVT_VFOR
#if defined(CYCLAGE_HOR) || defined(MVT_VFOR) || defined(MVT_GFOR)
//#define CYCLAGE_HOR //mode MVT_VFOR => cyclage horizontal sur le flux
#define CYCLAGE_MVT
#endif

extern int pbc_mode;  //periodic boundary conditions

#ifdef MVT_VFOR
float intensite_vfor = 0.5; //instensite du forcage vertical
#endif

#ifdef MVT_GFOR
float intensite_gfor = 0.001; //instensite du forcage global
#endif


extern unsigned char opt_info;
extern int       H, L, HL, HLP;       // les dimensions de la terre
extern int LN, LS, LNS, HLN;    //couloir est-ouest (limite nord, limite sud, largeur nord-sud ...)
extern Cell  *TE;	           // la 'terre'
extern int Ncel[];                       // nombre de cellules par type
extern unsigned char opt_vel;       //affichage des vitesses moyennes
extern Vec2 *norm2d;              //normale a la surface
//extern float *norm2dx;          // normale a la surface  (composante X)
//extern float *norm2dy;          // normale a la surface  (composante Y)
extern int h_ceil;  //epaisseur du plafond

MvtField *CelMvt;	//cellules en mouvements
MvtField *CelMvt2;	//cellules en mouvements
unsigned char Phase[MAX_CELL];	//phase (fluide ou solide) des types de cellules
MvtField Collisions[SIZE_MVT_FIELD];  //tableau des collisions fuilde-fluide
#ifdef SOLSLOW
MvtField Collisions_SolSlow[16][SIZE_MVT_FIELD]; //tableau des collisions fluide-solide (cas d'une unique particule lente)
#endif
#ifdef SOLFAST
MvtField Collisions_SolFast[16][SIZE_MVT_FIELD]; //tableau des collisions fluide-solide (cas d'une unique particule rapide)
#endif
//char propag_flag=0;  //semaphore (propagation en cours)
//char CelMvt_lock=0;	//mutex
//char dump_flag=0;  //semaphore (dump suspendu)
//char dump_inter=0;  //intervale entre dumps
int *Velx_time=NULL, *Vely_time=NULL;	//vitesses locales instantanees
int *Velx_sum=NULL, *Vely_sum=NULL;	//vitesses locales additives
int *Velx=NULL, *Vely=NULL;	//vitesses locales moyennees (espace et temps)
int *Velx_interp=NULL, *Vely_interp=NULL;	//vitesses locales interpolees (espace) et moyennees (temps)
//int CH, CLEO, CLNS, CHL, CHLP, CLN, CLS;     // dimensions du reseau de collisions (CLEO/CLNS = largeur du reseau de collisions en est-ouest/nord-sud)
int CH, CLEO, CLNS, CHL, CHLP;     // dimensions du reseau de collisions
int VH, VL, VHL, VHLP;	// dimensions du tableau des vitesses
int HLN;    //offset dans le reseau classique
long nb_cel_fluide=0;	//nombre de cellules fluides
long *NbMvt=NULL;		//nombre de particules en mouvement dans chaque couche du reseau
long nb_mvt_in=0;		//nombre de particules entrant dans le reseau
long nb_mvt_out=0;		//nombre de particules sortant du reseau
#ifdef MVT_REGUL
long *NbMvt0=NULL;		//nombre de particules en mouvement au depart dans chaque couche du reseau
#endif
unsigned char NbMvtCell[SIZE_MVT_FIELD];  //nombre de particules en mouvement dans chaque etat d'une cellule
float densite=0;	//densite moyenne de particules en mouvement
float maxvel = 0.0;  //vitesse interpolee maximale apres stabilisation
float meanvel = 0.0;  //vitesse interpolee moyenne apres stabilisation
int col_iter=0;   //nombre de cycles de collisions
char *nom_fic_mvt = NULL; //nom du fichier VELOC (initial)
int lgca_speedup = 0; //initial speedup of the stabilization of lattice gas
int lgca_mixing = 1; //mixing of the flux (on the left boundary)

void lecture_mvt();
void do_mvt_in(int ix);
void do_regul_mvt();
void do_force_mvt();
void do_mixing();
void compute_vel();
void compute_vel_interp_ns();
//void do_dump_mvt(int inter);
void dump_densite();
void dump_mvt_in_out();
void dump_collisions();


void params_collisions()
{
  //fichier MVT
  parametre("Hpp_file", "initial lattice gas - optional", &nom_fic_mvt, PARAM_STRING);
  parametre_cache("Mvt_file", "Mvt_file parameter is obsolete - Please use Hpp_file parameter instead !", &nom_fic_mvt, PARAM_STRING);

#ifdef MVT_GFOR
  parametre("Lgca_gfor", "global forcing of the flow - optional", &intensite_gfor, PARAM_FLOAT);
  parametre_cache("Mvt_gfor", "Mvt_gfor parameter is obsolete - Please use Lgca_gfor parameter instead !", &intensite_gfor, PARAM_FLOAT);
#endif

  parametre("Lgca_speedup", "initial speedup of the stabilization of lattice gas - optional", &lgca_speedup, PARAM_INT);
}

void init_collisions()
{
  int i;

  //initialisation tableau Phase[]
  memset(Phase, SOLIDE, MAX_CELL);
#ifdef MODELE_CMB
  Phase[MOINS] = Phase[ZERO] = FLUIDE;
  //Phase[ZERO] = FLUIDE;
  //Phase[MOINS] = -1;
#elif defined(MODELE_DUN)
  Phase[EAUC] = Phase[EAUT] = FLUIDE;
#endif

//#ifdef CYCLAGE_MVT
//  Phase[BORD] = -1;
//#else
  Phase[BORD] = SOLIDE; //-1;//SOLIDE;
//#endif

  //taille du reseau de collisions
  CH=H*NB_MVT_VER;
  CLEO=L*NB_MVT_EO;
  CLNS=(LNS+DIST_MVT_NS-1)/DIST_MVT_NS;
  CHL=CH*CLEO;

  //largeur du couloir
  //for(i=0; (i<L) && plan_solide(i); i++);
  //CLN=i;
  //for(; (i<L) && !plan_solide(i); i++);
  //CLS=i-1;
  //CLNS=CLS-CLN+1;
  CHLP=CH*CLEO*CLNS;

  LogPrintf("init_collisions\n");
  LogPrintf("NB_MVT_VER = %d\n", NB_MVT_VER);
  LogPrintf("NB_MVT_EO = %d\n", NB_MVT_EO);
  LogPrintf("DIST_MVT_NS = %d\n", DIST_MVT_NS);
  LogPrintf("CH = %d\n", CH);
  LogPrintf("CLEO = %d\n", CLEO);
  LogPrintf("CLNS = %d\n", CLNS);
  LogPrintf("VSTEP_H = %d\n", VSTEP_H);
  LogPrintf("VSTEP_L = %d\n", VSTEP_L);
  LogPrintf("VSTEP_TIME = %d\n", VSTEP_TIME);
  VH = (CH + VSTEP_H - 1)/VSTEP_H;
  VL = (CLEO + VSTEP_L - 1)/VSTEP_L;
  LogPrintf("VH = %d\n", VH);
  LogPrintf("VL = %d\n", VL);
  VHL = VH*VL;
  VHLP = VHL*CLNS;
#if (VSTEP == 1)
  LogPrintf("VMEAN_REF = %d\n", VMEAN_REF);
#endif
#ifdef MVT_REGUL
  LogPrintf("regulation flux\n");
#endif
#ifdef MVT_HFOR
  LogPrintf("forcage flux par tapis roulant\n");
#endif
#ifdef MVT_VFOR
  LogPrintf("forcage flux sur le bord gauche, intensite = %f\n", intensite_vfor);
#endif
#ifdef MVT_GFOR
  LogPrintf("forcage global du flux, intensite = %f\n", intensite_gfor);
#endif
  lgca_mixing = !pbc_mode; //mixing when cellular space is not periodic
  if (lgca_mixing) LogPrintf("brouillage flux sur le bord gauche\n");
#ifdef SOLSLOW
  LogPrintf("rebonds realistes SOLSLOW\n");
#endif
#ifdef SOLFAST
  LogPrintf("rebonds realistes SOLFAST\n");
#endif
#ifdef SF_REAL
  LogPrintf("rebonds realistes SF_REAL\n");
//  LogPrintf("epaisseur du plafond = %d\n", h_ceil);
#endif
#ifdef PLAF_REAL
  LogPrintf("rebonds realistes PLAF_REAL\n");
//  LogPrintf("epaisseur du plafond = %d\n", h_ceil);
#endif

  //allocation tableau CelMvt[]
  AllocMemoryPrint("CelMvt", CelMvt, MvtField, CHLP);
  ResetMemory(CelMvt, MvtField, CHLP);

  //allocation tableau CelMvt2[]
  AllocMemoryPrint("CelMvt2", CelMvt2, MvtField, CHLP);
  ResetMemory(CelMvt2, MvtField, CHLP);

  //allocation tableau NbMvt[]
  AllocMemory(NbMvt, long, CLNS);
  ResetMemory(NbMvt, long, CLNS);

#ifdef MVT_REGUL
  //allocation tableau NbMvt0[]
  AllocMemory(NbMvt0, long, CLNS);
  ResetMemory(NbMvt0, long, CLNS);
#endif
  PrintTotalMemory();

  if (nom_fic_mvt){
    lecture_mvt();
  }
  else{
    init_mvt();
  }

  // initialisation tableau NbMvtCell[]
  for (i=0; i<SIZE_MVT_FIELD; i++){
    NbMvtCell[i] = 0;
    if (i & MVT_E) NbMvtCell[i]++;
    if (i & MVT_O) NbMvtCell[i]++;
    if (i & MVT_B) NbMvtCell[i]++;
    if (i & MVT_H) NbMvtCell[i]++;
    if (i & MVT_EB) NbMvtCell[i]++;
    if (i & MVT_EH) NbMvtCell[i]++;
    if (i & MVT_OB) NbMvtCell[i]++;
    if (i & MVT_OH) NbMvtCell[i]++;
  }

// nombre de cellules mobiles
#ifdef MVT_REGUL
  for(k=0; k<CLNS; k++){
    if (NbMvt[k]) //TODO
      NbMvt0[k] = NbMvt[k];
    else //approximation
#ifdef MODELE_DUN
      NbMvt0[k] = Ncel[EAUC]*DENSITE/CLNS;
#else
      NbMvt0[k] = Ncel[MOINS]*DENSITE/CLNS;
#endif
    LogPrintf("NbMvt0[%d] = %ld\n", k, NbMvt0[k]);
  }
#endif

  LogPrintf("nb_mvt_in = %ld\n", nb_mvt_in);

  // initialisation tableau Collisions[]
  for (i=0; i<SIZE_MVT_FIELD; i++){
    if (i & MVT_OUT){
      Collisions[i] = MVT_OUT;
    }
    else if (i & MVT_SOLID){
      Collisions[i] = MVT_SOLID;
      if (i & MVT_E) Collisions[i] |= MVT_O;
      //if (i & MVT_E) Collisions[i] |= MVT_H; //test
      if (i & MVT_O) Collisions[i] |= MVT_E;
      if (i & MVT_B) Collisions[i] |= MVT_H;
      if (i & MVT_H) Collisions[i] |= MVT_B;

      if (i & MVT_EB) Collisions[i] |= MVT_OH; //MVT_OH;
      if (i & MVT_EH) Collisions[i] |= MVT_OB; //MVT_OB;
      if (i & MVT_OB) Collisions[i] |= MVT_EH; //MVT_EH;
      if (i & MVT_OH) Collisions[i] |= MVT_EB; //MVT_EB;

      /*if (i & MVT_EB) Collisions[i] |= MVT_EH; //MVT_OH;
      if (i & MVT_EH) Collisions[i] |= MVT_EB; //MVT_OB;
      if (i & MVT_OB) Collisions[i] |= MVT_OH; //MVT_EH;
      if (i & MVT_OH) Collisions[i] |= MVT_OB; //MVT_EB;*/

      /*if (i & MVT_EB) Collisions[i] |= MVT_OB; //MVT_OH;
      if (i & MVT_EH) Collisions[i] |= MVT_EB; //MVT_OB;
      if (i & MVT_OB) Collisions[i] |= MVT_OH; //MVT_EH;
      if (i & MVT_OH) Collisions[i] |= MVT_EH; //MVT_EB;*/
    }
    else{
      Collisions[i] = 0;

      if ((i & MVT_SLOW) == MVT_E){
        if ((i & MVT_FAST) == MVT_OB) Collisions[i] |= (MVT_O | MVT_EB);
	else if ((i & MVT_FAST) == MVT_OH) Collisions[i] |= (MVT_O | MVT_EH);
	else Collisions[i] |= MVT_E;
      }
      else if ((i & MVT_SLOW) == MVT_O){
        if ((i & MVT_FAST) == MVT_EB) Collisions[i] |= (MVT_E | MVT_OB);
	else if ((i & MVT_FAST) == MVT_EH) Collisions[i] |= (MVT_E | MVT_OH);
	else Collisions[i] |= MVT_O;
      }
      else if ((i & MVT_SLOW) == MVT_B){
        if ((i & MVT_FAST) == MVT_EH) Collisions[i] |= (MVT_H | MVT_EB);
	else if ((i & MVT_FAST) == MVT_OH) Collisions[i] |= (MVT_H | MVT_OB);
	else Collisions[i] |= MVT_B;
      }
      else if ((i & MVT_SLOW) == MVT_H){
        if ((i & MVT_FAST) == MVT_EB) Collisions[i] |= (MVT_B | MVT_EH);
	else if ((i & MVT_FAST) == MVT_OB) Collisions[i] |= (MVT_B | MVT_OH);
	else Collisions[i] |= MVT_H;
      }
      else if ((i & MVT_SLOW) == (MVT_E | MVT_O)) Collisions[i] |= (MVT_H | MVT_B);
      else if ((i & MVT_SLOW) == (MVT_H | MVT_B)) Collisions[i] |= (MVT_E | MVT_O);
      else Collisions[i] |= (i & MVT_SLOW);

      if ((i & MVT_FAST) == (MVT_EB | MVT_OH)) Collisions[i] |= (MVT_EH | MVT_OB);
      else if ((i & MVT_FAST) == (MVT_EH | MVT_OB)) Collisions[i] |= (MVT_EB | MVT_OH);
      else if (!(Collisions[i] & MVT_FAST)) Collisions[i] |= (i & MVT_FAST);
    }
  }

#ifdef SOLSLOW
  // initialisation (d'une partie) du tableau Collisions_SolSlow[]
  int mvt_sol_E, mvt_sol_O, mvt_sol_H, mvt_sol_B, mvt_sol;
  LogPrintf("collisions realistes : solslow\n");
  memset(Collisions_SolSlow,0,sizeof(Collisions_SolSlow));
  for (mvt_sol=0; mvt_sol<16; mvt_sol++){
    mvt_sol_E = mvt_sol & (MVT_E >> 2);
    mvt_sol_O = mvt_sol & (MVT_O >> 2);
    mvt_sol_B = mvt_sol & (MVT_B >> 2);
    mvt_sol_H = mvt_sol & (MVT_H >> 2);

    Collisions_SolSlow[mvt_sol][MVT_E] = (!mvt_sol_H) ? MVT_H : (!mvt_sol_B) ? MVT_B : MVT_O;
    Collisions_SolSlow[mvt_sol][MVT_O] = (!mvt_sol_H) ? MVT_H : (!mvt_sol_B) ? MVT_B : MVT_E;
    Collisions_SolSlow[mvt_sol][MVT_B] = (!mvt_sol_E) ? MVT_E : (!mvt_sol_O) ? MVT_O : MVT_H;
    //Collisions_SolSlow[mvt_sol][MVT_B] = MVT_H;
    Collisions_SolSlow[mvt_sol][MVT_H] = (!mvt_sol_E) ? MVT_E : (!mvt_sol_O) ? MVT_O : MVT_B;
  }
#endif

#ifdef SOLFAST
  // initialisation (d'une partie) du tableau Collisions_SolFast[]
  int mvt_sol_EB, mvt_sol_OB, mvt_sol_EH, mvt_sol_OH;//, mvt_sol;
  LogPrintf("collisions realistes : solfast\n");
  memset(Collisions_SolFast,0,sizeof(Collisions_SolFast));
  for (mvt_sol=0; mvt_sol<16; mvt_sol++){
    mvt_sol_EB = mvt_sol & (MVT_EB >> 6);
    mvt_sol_EH = mvt_sol & (MVT_EH >> 6);
    mvt_sol_OB = mvt_sol & (MVT_OB >> 6);
    mvt_sol_OH = mvt_sol & (MVT_OH >> 6);

    Collisions_SolFast[mvt_sol][MVT_EB] = (!mvt_sol_EH) ? MVT_EH : (!mvt_sol_OB) ? MVT_OB : MVT_OH;
    Collisions_SolFast[mvt_sol][MVT_EH] = (!mvt_sol_OH) ? MVT_OH : (!mvt_sol_EB) ? MVT_EB : MVT_OB;
    Collisions_SolFast[mvt_sol][MVT_OB] = (!mvt_sol_OH) ? MVT_OH : (!mvt_sol_EB) ? MVT_EB : MVT_EH;
    Collisions_SolFast[mvt_sol][MVT_OH] = (!mvt_sol_EH) ? MVT_EH : (!mvt_sol_OB) ? MVT_OB : MVT_EB;
    //Collisions_SolFast[mvt_sol][MVT_EB] = MVT_EH;
    //Collisions_SolFast[mvt_sol][MVT_EH] = MVT_EB;
    //Collisions_SolFast[mvt_sol][MVT_OB] = MVT_OH;
    //Collisions_SolFast[mvt_sol][MVT_OH] = MVT_OB;
  }
#endif

  if (opt_info) dump_collisions();

  //dump_densite();

#ifdef DUMP_SIGNATURE
  //empreinte gazeuse
  if (opt_info) dump_signature_mvt();
#endif
}

// initialisation du tableau MvtCel des particules de gaz sur reseau
void init_mvt()
{
  int i, j, k, ic, alea;
  Cell *aux;
  MvtField *aux_mvt;

  LogPrintf("initialisation du gaz sur reseau\n");

  // initialisation tableau CelMvt[]
  MvtField mvt_slow[] = {MVT_E, MVT_O, MVT_B, MVT_H};
  MvtField mvt_fast[] = {MVT_EB, MVT_EH, MVT_OB, MVT_OH};
  aux = TE+HLN;
  aux_mvt = CelMvt;
  nb_cel_fluide = 0;
  for(k=0; k<CLNS; k++, aux+=HL*(DIST_MVT_NS-1)){ //profondeur
    for(j=0; j<H; j++){ //hauteur
      for(i=0; i<L; i++, aux++){ //largeur
        for(ic=0; ic<NB_MVT_EO; ic++, aux_mvt++){
          if (Phase[aux->celltype] == FLUIDE){
            *aux_mvt=0;
            nb_cel_fluide++;
            //if (drand48() < 0.5)
            /*if (j>0.6*H)*/{
              //do_mvt_in(k*HL+j*L+i);
              /* *mvt = MVT_E;
              nb_mvt++;
              float aleat = drand48();
              if (aleat < 0.25) {*mvt |= MVT_EB; nb_mvt++;}
              else if (aleat < 0.5) {*mvt |= MVT_EH; nb_mvt++;}*/

              // 1 particule
              /*int alea8=drand48()*8;
              *aux_mvt = (alea8<1) ? MVT_E : 0;
              // *aux_mvt = (alea8<4) ? mvt_slow[alea8] : (alea8<8) ? mvt_fast[alea8-4] : 0;
                      if (*aux_mvt) nb_mvt_in++;*/

              // 1 particule lente
              alea=drand48()*4;
              if (alea<4){
                *aux_mvt = mvt_slow[alea];
                nb_mvt_in++;
              }
              //*aux_mvt = MVT_E;
              //*aux_mvt = (alea4) ? MVT_E : 0;

              // 1 particule rapide
              //alea4=drand48()*4;
              //*aux_mvt |= mvt_fast[alea4];
              //nb_mvt_in +=2;

              // 1/2 particule rapide
              alea=drand48()*8;
              if (alea<4){ // 1/2 particule rapide
                *aux_mvt |= mvt_fast[alea];
                nb_mvt_in++;
              }

              // 2 particules lentes
              /*int alea2=drand48()*2;
              *aux_mvt = mvt_slow[alea2]; // 1 particule lente horizontale
              alea2=drand48()*2;
              *aux_mvt |= mvt_slow[2+alea2]; // 1 particule lente verticale
              nb_mvt_in +=2;*/

              // 2 particules rapides
              /*alea2=drand48()*2;
              *aux_mvt |= mvt_fast[alea2]; // 1 particule rapide vers le bas
              alea2=drand48()*2;
              *aux_mvt |= mvt_fast[2+alea2]; // 1 particule rapide vers le haut
              nb_mvt_in +=2;*/

              // 6 particules
              /*int alea4=drand48()*4;
              *aux_mvt = MVT_SLOW ^ mvt_slow[alea4]; // 3 particules lentes
              alea4=drand48()*4;
              *aux_mvt |= MVT_FAST ^ mvt_fast[alea4]; // 3 particules rapides
              nb_mvt_in +=6;*/
            }
          }
          else if (Phase[aux->celltype] == SOLIDE)
            *aux_mvt = MVT_SOLID;

          //if ((j==0) || (j==H-1)) //bords horizontaux
          //if ((j==1) || (j==H-2)) //bords horizontaux
          //  *aux_mvt = MVT_SOLID;
//#if defined(MVT_REGUL) && !defined(CYCLAGE_HOR)
#if !defined(CYCLAGE_MVT)
#ifdef MVT_REGUL
          else if ((i==0) || (i==L-1)) //bords est-ouest
            *aux_mvt = MVT_OUT;
#else
          else if (i==L-1) // bord est
            *aux_mvt = MVT_OUT;
#endif
#endif
        }
      }
    }
  }

  //speed up the stabilization of the flow
  if (lgca_speedup){
    //int nb_cyc_for=VSTEP_TIME*10;
    LogPrintf("speed up of the lattice gas: %d\n", lgca_speedup);
    for(i=0; i<lgca_speedup; i++) do_force_mvt();
  }
}

// modification globale de la terre
void collisions_mod_terre()
{
  int i, j, k, ic;
  Cell *aux;
  MvtField *aux_mvt;

  nb_cel_fluide = 0;
  aux = TE+HLN;
  aux_mvt = CelMvt;
  for(k=0; k<CLNS; k++, aux+=HL*(DIST_MVT_NS-1)){ //profondeur
    for(j=0; j<H; j++){ //hauteur
      for(i=0; i<L; i++, aux++){ //largeur
        for(ic=0; ic<NB_MVT_EO; ic++, aux_mvt++){
          if (Phase[aux->celltype] == SOLIDE){
            SetMask(*aux_mvt, MVT_SOLID);
          }
          else{
            UnsetMask(*aux_mvt, MVT_SOLID);
            nb_cel_fluide++;
          }
        }
      }
    }
  }
}

void lecture_mvt()
{
  FILE *fp;
  int i, ix;

  LogPrintf("lecture fichier HPP : %s\n",nom_fic_mvt);

  fp = fopen(nom_fic_mvt,"r");
  if ( ! fp ){
    ErrPrintf("ERROR: cannot open HPP file %s\n", nom_fic_mvt);
    exit(-4);
  }

  fread(CelMvt, sizeof(MvtField), CHLP, fp);

  fclose(fp);

  // initialisation nb_cel_fluide
  for (ix=0; ix<CHLP; ix++){
    if (!(CelMvt[ix] & MVT_SOLID)) nb_cel_fluide++;
  }
  /*nb_cel_fluide = 0;*/

  // moyennage de la vitesse initiale
  /*for (i=0; i<VSTEP_TIME; i++)*/ compute_vel();
}


int mvt(int index)
{
  //return (CelMvt[index] & (MVT_E | MVT_O | MVT_B | MVT_H));
  return (CelMvt[index] & MVT_ALLDIR);
  //return (NbMvtCell[CelMvt[index]] > densite);
}

int mvt_mean(int index)
{
//#if (VSTEP == 1)
#if 0
  int vx = Velx_sum[index];
  int vy = Vely_sum[index];
  return (vx*vx+vy*vy) > VMEAN_REF; /* /64 */
#else
  return mvt(index);
#endif
}

int mvt_solid(int ii)
{
  return (CelMvt[ii] & MVT_SOLID);
}

int mvt_bas(int ii)
{
  return (CelMvt[ii] & (MVT_B | MVT_EB | MVT_OB));
}

int mvt_haut(int ii)
{
  return (CelMvt[ii] & (MVT_H | MVT_EH | MVT_OH));
}

int mvt_est(int ii)
{
  //return (CelMvt[ii] & (MVT_EB | MVT_EH));
  return (CelMvt[ii] & (MVT_E | MVT_EB | MVT_EH));
}

int mvt_ouest(int ii)
{
  //return (CelMvt[ii] & MVT_O);
  //return (CelMvt[ii] & (MVT_OB | MVT_OH));
  return (CelMvt[ii] & (MVT_O | MVT_OB | MVT_OH));
}

int check_mvt(int index, void *data)
{
  int ii;
  Calcule_cix(index, ii);
  return (CelMvt[ii] & MVT_ALLDIR);
}

int check_mvt_solid(int index, void *data)
{
  int ii;
  Calcule_cix(index, ii);
  return (CelMvt[ii] & MVT_SOLID);
}

int check_mvt_bas(int index, void *data)
{
  int ii;
  Calcule_cix(index, ii);
  /*
  int val = 0;
  if (CelMvt[index] & MVT_B) val++;
  if (CelMvt[index] & MVT_EB) val++;
  if (CelMvt[index] & MVT_OB) val++;
  if (CelMvt[index] & MVT_H) val--;
  if (CelMvt[index] & MVT_EH) val--;
  if (CelMvt[index] & MVT_OH) val--;

  return (val > 0);
  */
  return (CelMvt[ii] & (MVT_B | MVT_EB | MVT_OB));
}

int check_mvt_haut(int index, void *data)
{
  int ii;
  Calcule_cix(index, ii);
  /*
  int val = 0;
  if (CelMvt[ii] & MVT_H) val++;
  if (CelMvt[ii] & MVT_EH) val++;
  if (CelMvt[ii] & MVT_OH) val++;
  if (CelMvt[ii] & MVT_B) val--;
  if (CelMvt[ii] & MVT_EB) val--;
  if (CelMvt[ii] & MVT_OB) val--;

  return (val > 0);
  */
  return (CelMvt[ii] & (MVT_H | MVT_EH | MVT_OH));
}

int check_mvt_est(int index, void *data)
{
  int ii;
  Calcule_cix(index, ii);
  /*
  int val = 0;
  if (CelMvt[index] & MVT_E) val++;
  if (CelMvt[index] & MVT_EB) val++;
  if (CelMvt[index] & MVT_EH) val++;
  if (CelMvt[index] & MVT_O) val--;
  if (CelMvt[index] & MVT_OB) val--;
  if (CelMvt[index] & MVT_OH) val--;

  return (val > 0);
  */
  return (CelMvt[ii] & (MVT_E | MVT_EB | MVT_EH));
  //return (CelMvt[ii] & MVT_E);
}

int check_mvt_ouest(int index, void *data)
{
  int ii;
  Calcule_cix(index, ii);
  /*int val = 0;
  if (CelMvt[ii] & MVT_O) val++;
  if (CelMvt[ii] & MVT_OB) val++;
  if (CelMvt[ii] & MVT_OH) val++;
  if (CelMvt[ii] & MVT_E) val--;
  if (CelMvt[ii] & MVT_EB) val--;
  if (CelMvt[ii] & MVT_EH) val--;

  return (val > 0);*/
  return (CelMvt[ii] & (MVT_O | MVT_OB | MVT_OH));
}

int check_mvt_est_et_bas(int index, void *data)
{
  return check_mvt_est(index, NULL) && check_mvt_bas(index, NULL);
}

int check_mvt_EB(int index, void *data)
{
  int ii;
  Calcule_cix(index, ii);
  //return (CelMvt[ii] & MVT_B | MVT_E | MVT_EB); //bug
  //return (CelMvt[ii] & (MVT_B | MVT_E | MVT_EB));
  return (CelMvt[ii-10*CLEO] & (MVT_B | MVT_E | MVT_EB)); //magouille
}

int check_no_mvt_est_et_bas(int index, void *data)
{
  return !check_mvt_est_et_bas(index, NULL);
}

int check_mvt_mean_bas(int index, void *data)
{
#if (VSTEP > 1)
  int cx, cy, cz;
  Calcule_cxyz(index, cx, cy, cz);
  return (Vely[(cx/VSTEP_L)+(cy/VSTEP_H)*VL+cz*VHL] > 0);
#else
  return (Vely[index] >0);
#endif
}

int check_mvt_mean_est(int index, void *data)
{
#if (VSTEP > 1)
  int cx, cy, cz;
  Calcule_cxyz(index, cx, cy, cz);
  return (Velx[(cx/VSTEP_L)+(cy/VSTEP_H)*VL+cz*VHL] > 0);
#else
  return (Velx[index] >0);
#endif
}

int check_no_mvt_mean_bas(int index, void *data)
{
  return !check_mvt_mean_bas(index, NULL);
}

void collisions_modcell(int type, int index)
{
  //int ii=(index-HLN)*NB_MVT_EO;
  int x, y, z, ii, ic, cz;
  Calcule_xyz(index, x, y, z);
  cz = (z-LN)/DIST_MVT_NS;
  if (z != cz*DIST_MVT_NS + LN){
    //ErrPrintf("collisions_modcell : mauvais plan ! (%d)\n", z);
    return;
  }
  ii = x*NB_MVT_EO + y*CLEO + cz*CHL;
  //Calcul_cix(index, ii);

  for(ic=0; ic<NB_MVT_EO; ic++, ii++){
    if (Phase[type] == SOLIDE){
      if (!(CelMvt[ii] & MVT_SOLID)) nb_cel_fluide--;
      SetMask(CelMvt[ii], MVT_SOLID);
      //if (mvt(index)) {LogPrintf("vert piege !!\n");}
    }
    else{
      if (CelMvt[ii] & MVT_SOLID) nb_cel_fluide++;
      UnsetMask(CelMvt[ii], MVT_SOLID);
    }
  }
}


void do_collisions()
{
  int i, j, k, ii, ii2;

#ifdef SF_REAL
  //const float PI = 3.14159;
  // angle limite pour les rebonds = PI/8
  // si l'angle avec la normale < angle limite, alors direction inversee, sinon changement de direction a angle droit
  //const float TP8 = tan(PI/8.0); //~0.41
  const float CP8 = cos(PI/8.0);
  const float SP8 = sin(PI/8.0);
  char sfmod = 0;
#endif


  //LogPrintf("do_collisions\n");
  for(ii=k=0; k<CLNS; k++){ // profondeur
    //LogPrintf("k=%d ii=%d\n", k, ii);
    ii+=CLEO;
    for(j=1; j<CH-1; j++){ // hauteur
      ii += NB_MVT_EO;
      for(i=NB_MVT_EO; i<CLEO-NB_MVT_EO; i++, ii++){ //largeur
        CelMvt2[ii] = Collisions[CelMvt[ii]];
        //CelMvt2[ii] = CelMvt[ii];

#ifdef PLAF_REAL
        if ((j <= h_ceil) && (NbMvtCell[CelMvt[ii] & MVT_FAST]==1)){
          //rebond realiste au plafond
          MvtField mvt2 = 0;
          switch (CelMvt[ii] & MVT_FAST)
          {
            case MVT_EH:
              mvt2 = MVT_EB;
              break;
            case MVT_OH:
              mvt2 = MVT_OB;
              break;
          }
          if (mvt2){
            UnsetMask(CelMvt2[ii], MVT_FAST);
            SetMask(CelMvt2[ii], mvt2);
          }
        }
#endif

  //int mvt_sol_slow = 0;//debug
  //int mvt_sol_fast = 0;//debug
  //int mvt_sfr_slow = 0;//debug
  //int mvt_sfr_fast = 0;//debug

#ifdef SOLSLOW
        if (mvt_solid(ii) && (NbMvtCell[CelMvt[ii] & MVT_SLOW]==1)){ //Collisions ameliorees (fuilde-solide)
          //EST
          ii2 = ii+1;
#ifdef CYCLAGE_HOR
          if (pbc_mode && (i==CLEO-NB_MVT_EO-1)) ii2 -= CLEO-2*NB_MVT_EO;
#endif
          //OUEST
          int mvt_sol_E = mvt_solid(ii2) ? (MVT_E >> 2) : 0;
          ii2 = ii-1;
#ifdef CYCLAGE_HOR
          if (pbc_mode && (i==NB_MVT_EO)) ii2 += CLEO-2*NB_MVT_EO;
#endif
          int mvt_sol_O = mvt_solid(ii-1) ? (MVT_O >> 2) : 0;
          //BAS
          int mvt_sol_B = mvt_solid(ii+CLEO) ? (MVT_B >> 2) : 0;
          //HAUT
          int mvt_sol_H = mvt_solid(ii-CLEO) ? (MVT_H >> 2) : 0;
          int mvt_sol = mvt_sol_E | mvt_sol_O | mvt_sol_B | mvt_sol_H;
          UnsetMask(CelMvt2[ii], MVT_SLOW);
          SetMask(CelMvt2[ii], Collisions_SolSlow[mvt_sol][CelMvt[ii] & MVT_SLOW]);
          //mvt_sol_slow = Collisions_SolSlow[mvt_sol][CelMvt[ii] & MVT_SLOW];  //debug
        }
#endif

#ifdef SOLFAST
        //if ((CelMvt[ii] & MVT_SOLID) && (NbMvtCell[CelMvt[ii] & MVT_FAST]==1)){ //Collisions ameliorees (fuilde-solide)
        if (mvt_solid(ii) && (NbMvtCell[CelMvt[ii] & MVT_FAST]==1)){ //Collisions ameliorees (fuilde-solide)
        //if (mvt_solid(ii) && (NbMvtCell[CelMvt[ii] & MVT_FAST]>=1)){ //Collisions ameliorees (fuilde-solide)
          /*int mvt_sol_E = CelMvt[ii+1] & MVT_SOLID;
          int mvt_sol_O = CelMvt[ii-1] & MVT_SOLID;
          int mvt_sol_B = CelMvt[ii+CLEO] & MVT_SOLID;
          int mvt_sol_H = CelMvt[ii-CLEO] & MVT_SOLID;*/

          //EST-BAS
          ii2 = ii+CLEO+1;
#ifdef CYCLAGE_HOR
          if (pbc_mode && (i==CLEO-NB_MVT_EO-1)) ii2 -= CLEO-2*NB_MVT_EO;
#endif
          int mvt_sol_EB = mvt_solid(ii2) ? (MVT_EB >> 6) : 0;  //CelMvt[ii+CLEO+1] & MVT_SOLID;

          //EST-HAUT
          ii2 = ii-CLEO+1;
#ifdef CYCLAGE_HOR
          if (pbc_mode && (i==CLEO-NB_MVT_EO-1)) ii2 -= CLEO-2*NB_MVT_EO;
#endif
          int mvt_sol_EH = mvt_solid(ii2) ? (MVT_EH >> 6) : 0;  //CelMvt[ii-CLEO+1] & MVT_SOLID;

          //OUEST-BAS
          ii2 = ii+CLEO-1;
#ifdef CYCLAGE_HOR
          if (pbc_mode && (i==NB_MVT_EO)) ii2 += CLEO-2*NB_MVT_EO;
#endif
          int mvt_sol_OB = mvt_solid(ii2) ? (MVT_OB >> 6) : 0;  //CelMvt[ii+CLEO-1] & MVT_SOLID;

          //OUEST-HAUT
          ii2 = ii-CLEO-1;
#ifdef CYCLAGE_HOR
          if (pbc_mode && (i==NB_MVT_EO)) ii2 += CLEO-2*NB_MVT_EO;
#endif
          int mvt_sol_OH = mvt_solid(ii2) ? (MVT_OH >> 6) : 0;  //CelMvt[ii-CLEO-1] & MVT_SOLID;
          //int mvt_sol = (mvt_sol_E << 3) | (mvt_sol_O << 2) | (mvt_sol_B << 1) | mvt_sol_H;
          //int mvt_sol = (mvt_sol_EB << 3) | (mvt_sol_OB << 2) | (mvt_sol_EH << 1) | mvt_sol_OH;
          int mvt_sol = mvt_sol_EB | mvt_sol_EH | mvt_sol_OB | mvt_sol_OH;
          UnsetMask(CelMvt2[ii], MVT_FAST);
          SetMask(CelMvt2[ii], Collisions_SolFast[mvt_sol][CelMvt[ii] & MVT_FAST]);
          //mvt_sol_fast = Collisions_SolFast[mvt_sol][CelMvt[ii] & MVT_FAST];  //debug
        }
#endif

#ifdef SF_REAL  //Collisions realistes (solide-fluide)
        const float slow_val = SP8; //0.01; //SP8; //~0.38
        const float fast_val = CP8-SP8; //0.01; //CP8-SP8; //~0.54
        // normale precalculee
        int ni = (i/NB_MVT_EO) - 1 + (k*DIST_MVT_NS)*(L-2);
        //LogPrintf("ni = %d \n", ni);
        float nx = norm2d[ni].x;
        float ny = norm2d[ni].y;
        //cas d'une particule lente
        //if (0){
        if (mvt_solid(ii) && (NbMvtCell[CelMvt[ii] & MVT_SLOW]==1)){
          MvtField mvt2 = 0;
          sfmod = 0;
          // calcul de la direction apres collision
          switch (CelMvt[ii] & MVT_SLOW)
          {
            case MVT_E:
              //LogPrintf("i=%d  k=%d   nx = %f   ny = %f\n", i, k, nx, ny);
              //if (fabs(ny) > slow_coef*fabs(nx)){
              if (fabs(ny) > slow_val){
              //if (ny){
                mvt2 = (ny > 0) ? MVT_H : MVT_B;
                sfmod = 1;
              }
              break;

            case MVT_O:
              if (fabs(ny) > slow_val){
              //if (ny){
                mvt2 = (ny > 0) ? MVT_H : MVT_B;
                sfmod = 1;
              }
              break;

            case MVT_B:
              /*if (fabs(nx) > slow_val){
              //if (nx){
                mvt2 = (nx > 0) ? MVT_E : MVT_O;
                sfmod = 1;
          }*/
              break;

            case MVT_H:
              //on conserve le rebond classique (MVT_B), on suppose qu'il s'agit du plafond ...
              break;
          }
          if (sfmod){
            // verification de l'absence de particule solide dans la direction calculee
            switch (mvt2)
            {
              case MVT_E:
                ii2 = ii+1;
#ifdef CYCLAGE_HOR
                if (pbc_mode && (i==CLEO-NB_MVT_EO-1)) ii2 -= CLEO-2*NB_MVT_EO;
#endif
                break;
              case MVT_O:
                ii2 = ii-1;
#ifdef CYCLAGE_HOR
                if (pbc_mode && (i==NB_MVT_EO)) ii2 += CLEO-2*NB_MVT_EO;
#endif
                break;
              case MVT_B:
                ii2 = ii+CLEO;
                break;
              case MVT_H:
                ii2 = ii-CLEO;
                break;
            }
            if (mvt_solid(ii2)) sfmod = 0;
          }
          if (sfmod){
            // modification des collisions par defaut
            //LogPrintf("collision modifiee : %d\n", mvt2);
            UnsetMask(CelMvt2[ii], MVT_SLOW);
            SetMask(CelMvt2[ii], mvt2);
            //mvt_sfr_slow = mvt2;  //debug
          }
        }

        //cas d'une particule rapide
        if (mvt_solid(ii) && (NbMvtCell[CelMvt[ii] & MVT_FAST]==1)){
          MvtField mvt2 = 0;
          sfmod = 0;
          // calcul de la direction apres collision
          switch (CelMvt[ii] & MVT_FAST)
          {
            case MVT_EB:
              if (fabs(nx+ny) > fast_val){
              //if (nx+ny){
              //if (1){
                mvt2 = (nx+ny > 0) ? MVT_EH : MVT_OB;
                sfmod = 1;
              }
              break;

            case MVT_EH:
              //if (!(ny-nx))
              //  continue; //on conserve le rebond classique (MVT_OB)
              //else
                //mvt2 = (ny-nx > 0) ? MVT_EH : MVT_OB;
              mvt2 = (j <= h_ceil) ? MVT_EB : MVT_OH; //selon l'altitude on considere que l'on rebondit sur le plafond ou sur une surface quasi-verticale
              sfmod = 1;
              break;

            case MVT_OB:
              if (fabs(ny-nx) > fast_val){
              //if (ny-nx){
                mvt2 = (ny-nx > 0) ? MVT_OH : MVT_EB;
                sfmod = 1;
              }
              break;

            case MVT_OH:
              mvt2 = (j <= h_ceil) ? MVT_OB : MVT_EH; //selon l'altitude on considere que l'on rebondit sur le plafond ou sur une surface quasi-verticale
              sfmod = 1;
          }
          if (sfmod){
            // verification de l'absence de particule solide dans la direction calculee
            switch (mvt2)
            {
              case MVT_EB:
                ii2=ii+1+CLEO;
#ifdef CYCLAGE_HOR
                if (pbc_mode && (i==CLEO-NB_MVT_EO-1)) ii2 -= CLEO-2*NB_MVT_EO;
#endif
                break;
              case MVT_EH:
                ii2=ii+1-CLEO;
#ifdef CYCLAGE_HOR
                if (pbc_mode && (i==CLEO-NB_MVT_EO-1)) ii2 -= CLEO-2*NB_MVT_EO;
#endif
                break;
              case MVT_OB:
                ii2=ii-1+CLEO;
#ifdef CYCLAGE_HOR
                if (pbc_mode && (i==NB_MVT_EO)) ii2 += CLEO-2*NB_MVT_EO;
#endif
                break;
              case MVT_OH:
                ii2=ii-1-CLEO;
#ifdef CYCLAGE_HOR
                if (pbc_mode && (i==NB_MVT_EO)) ii2 += CLEO-2*NB_MVT_EO;
#endif
                break;
            }
            if (mvt_solid(ii2)) sfmod = 0;
          }
          if (sfmod){
            // modification des collisions par defaut
            //LogPrintf("collision modifiee : %d (k=%d   j=%d   i=%d)\n", mvt2, k, j, i);
            UnsetMask(CelMvt2[ii], MVT_FAST);
            SetMask(CelMvt2[ii], mvt2);
            //mvt_sfr_fast = mvt2;  //debug
          }
        }
#endif

        //comparaison SOL vs SF_REAL
        //if (mvt_sol_slow != mvt_sfr_slow) {LogPrintf("col_iter = %d   ii = %d   CelMvt[ii] = %d   CelMvt2[ii] = %d   mvt_sol_slow = %d   mvt_sfr_slow = %d   nx = %f   ny = %f   ni = %d\n", col_iter, ii, CelMvt[ii], CelMvt2[ii], mvt_sol_slow, mvt_sfr_slow, nx, ny, ni);}
        //if (mvt_sol_fast != mvt_sfr_fast) {LogPrintf("col_iter = %d   ii = %d   CelMvt[ii] = %d   CelMvt2[ii] = %d   mvt_sol_fast = %d   mvt_sfr_fast = %d   nx = %f   ny = %f   ni = %d\n", col_iter, ii, CelMvt[ii], CelMvt2[ii], mvt_sol_fast, mvt_sfr_fast, nx, ny, ni);}

#if !defined(CYCLAGE_MVT)
        //if ((CelMvt[ix] & MVT_OUT) && mvt(ix)) nb_mvt--;
        if ((CelMvt[ii] & MVT_OUT) && mvt(ii)){
          NbMvt[k] -= NbMvtCell[CelMvt[ii]];
          nb_mvt_out += NbMvtCell[CelMvt[ii]];
        }
    //if ((CelMvt[ix] & MVT_OUT) && (CelMvt[ix] & (MVT_E | MVT_O | MVT_B | MVT_H))) nb_mvt--;
#endif
      }
      ii += NB_MVT_EO;
    }
    ii+=CLEO;
  }
  col_iter++;
}

void do_propagations()
{
  int i, j, k, ii, ii0;
  int n=0;

  //while (CelMvt_lock) {LogPrintf("lock\n");} //c devenu inutile ...

  //CelMvt_lock = 1;

  //memset(CelMvt, 0, sizeof(unsigned char) * HL);

  ii=0;
  for(k=0; k<CLNS; k++) // profondeur
    for(j=0; j<CH; j++) // hauteur
      for(i=0; i<CLEO; i++, ii++){ //largeur
        if ((i>=NB_MVT_EO) && (i<CLEO-NB_MVT_EO) && (j>0) && (j<CH-1)){ //on ne propage pas sur les bords ...
          CelMvt[ii]=0;
          //SOLID
          if (CelMvt2[ii] & MVT_SOLID)
            CelMvt[ii] |= MVT_SOLID;
          //OUT
          if (CelMvt2[ii] & MVT_OUT)
            CelMvt[ii] |= MVT_OUT;
          //EST
          ii0 = ii-1;
#ifdef CYCLAGE_MVT
          if (i==NB_MVT_EO) ii0 += CLEO-2*NB_MVT_EO;
#endif
          if (CelMvt2[ii0] & MVT_E)
            CelMvt[ii] |= MVT_E;
          //else
            //CelMvt[ii] &= (-1) ^ MVT_E;

          //OUEST
          ii0 = ii+1;
#ifdef CYCLAGE_MVT
          if (i==CLEO-NB_MVT_EO-1) ii0 -= CLEO-2*NB_MVT_EO;
#endif
          if (CelMvt2[ii0] & MVT_O)
            CelMvt[ii] |= MVT_O;
          //else
            //CelMvt[ii] &= (-1) ^ MVT_O;

          //BAS
          ii0 = ii-CLEO;
          if (CelMvt2[ii0] & MVT_B)
            CelMvt[ii] |= MVT_B;

          //HAUT
          ii0 = ii+CLEO;
          if (CelMvt2[ii0] & MVT_H)
            CelMvt[ii] |= MVT_H;

          //EST+BAS
          ii0 = ii-CLEO-1;
#ifdef CYCLAGE_MVT
          if (i==NB_MVT_EO) ii0 += CLEO-2*NB_MVT_EO;
#endif
          if (CelMvt2[ii0] & MVT_EB)
            CelMvt[ii] |= MVT_EB;

          //EST+HAUT
          ii0 = ii+CLEO-1;
#ifdef CYCLAGE_MVT
          if (i==NB_MVT_EO) ii0 += CLEO-2*NB_MVT_EO;
#endif
          if (CelMvt2[ii0] & MVT_EH)
            CelMvt[ii] |= MVT_EH;

          //OUEST+BAS
          ii0 = ii-CLEO+1;
#ifdef CYCLAGE_MVT
          if (i==CLEO-NB_MVT_EO-1) ii0 -= CLEO-2*NB_MVT_EO;
#endif
          if (CelMvt2[ii0] & MVT_OB)
            CelMvt[ii] |= MVT_OB;

          //OUEST+HAUT
          ii0 = ii+CLEO+1;
#ifdef CYCLAGE_MVT
          if (i==CLEO-NB_MVT_EO-1) ii0 -= CLEO-2*NB_MVT_EO;
#endif
          if (CelMvt2[ii0] & MVT_OH)
            CelMvt[ii] |= MVT_OH;
        }
        //else
          //CelMvt[ii] = MVT_SOLID; //bords
        //if ((CelMvt[ii] & MVT_SOLID) && (mvt(ii))) {LogPrintf("vert piege !!\n"); }
      }

  //CelMvt_lock = 0;


  //Retroaction transport/matiere
  /*
  ix = 0;
  for(j=0; j<H; j++){ //hauteur
    for(i=0; i<L; i++, ix++){ //largeur
      if (CelMvt[ix] & ((-1) ^ MVT_SOLID)) { if (TE[ix].celltype == MOINS) {TE[ix].celltype = ZERO; n++;} }
      else if (TE[ix].celltype == ZERO) {TE[ix].celltype = MOINS; n++;}
    }
  }
  */

  //LogPrintf("propagation de %d cellules\n", n);

#ifdef MVT_REGUL
  do_regul_mvt();
#endif
#if defined(MVT_HFOR) || defined(MVT_VFOR) || defined(MVT_GFOR) //forcage flux
  do_force_mvt();
#endif

  if (lgca_mixing) do_mixing();

#ifdef CGV
  compute_vel();
#else
  /*if (opt_vel)*/ compute_vel();
#endif

#ifdef DUMP_SIGNATURE
  //empreinte gazeuse
  if (opt_info && !(col_iter & 0xf))dump_signature_mvt(); //1 cycle sur 16
#endif

  if (opt_info) dump_mvt_in_out();
}

void do_mvt_in(int ix)
{
  float aleat = drand48();
  int k = (int)(ix/CHL);
  //MvtField mvt0 = CelMvt[ix];
  //injection d'une particule lente
  if (aleat < 0.5*DENSITE)
  if (!(CelMvt[ix] & (MVT_SOLID | MVT_OUT | MVT_E | MVT_O))){
    CelMvt[ix] |= MVT_E;
    NbMvt[k]++;
    nb_mvt_in++;
  }

  //injection d'une particule rapide
  if (aleat < 0.25*DENSITE){
    if (!(CelMvt[ix] & (MVT_SOLID | MVT_OUT | MVT_EB | MVT_OH))){
      CelMvt[ix] |= MVT_EB;
      NbMvt[k]++;
      nb_mvt_in++;
    }
  }
  else if (aleat < 0.5*DENSITE){
    if (!(CelMvt[ix] & (MVT_SOLID | MVT_OUT | MVT_EH | MVT_OB))){
      CelMvt[ix] |= MVT_EH;
      NbMvt[k]++;
      nb_mvt_in++;
    }
  }
}

#ifdef MVT_REGUL
void do_regul_mvt()
{
  int k, ix, y, cpt;

  //LogPrintf("do_regul_mvt : %d -> %d\n", NbMvt, NbMvt0);

  for(k=0; k<CLNS; k++){
    //on injecte des particules sur le bord gauche
    cpt = 0;
    while ((NbMvt[k] < NbMvt0[k]) && (cpt<CH*2/*CH/2*/)){
      y = CH*drand48();
      ix = NB_MVT_EO+y*CLEO+k*CHL;
      do_mvt_in(ix);
      cpt++;
    }

    if (NbMvt[k] < NbMvt0[k]*0.9) LogPrintf("NbMvt[%d] = %ld\n", k, NbMvt[k])
  }

  //LogPrintf("fin regul\n");
}
#endif

void do_force_mvt()
{
  int x, y, z;
  int ix;
  // injection flux maximum sur le bord gauche
  /*for (y=0; y<CH; y++){
    ix = 1+y*CLEO;
    do_mvt_in(ix);
  }*/
#ifdef MVT_HFOR
  // forcage flux ouest-est sur la couche superieure
  //for(y=1; y<CH/2; y++)

  for(z=0; z<CLNS; z++){
    for(x=0; x<CLEO; x++){
      ix = x+CLEO*NB_MVT_VER*(h_ceil+1)+z*CHL;
      if ((CelMvt[ix] & (MVT_E | MVT_O)) == MVT_O){
        UnsetMask(CelMvt[ix], MVT_O);
        SetMask(CelMvt[ix], MVT_E);
      }
      if ((CelMvt[ix] & (MVT_EH | MVT_OH)) == MVT_OH){
        UnsetMask(CelMvt[ix], MVT_OH);
        SetMask(CelMvt[ix], MVT_EH);
      }
      if ((CelMvt[ix] & (MVT_EB | MVT_OB)) == MVT_OB){
        UnsetMask(CelMvt[ix], MVT_OB);
        SetMask(CelMvt[ix], MVT_EB);
      }
    }
  }


  // forcage flux est-ouest sur la couche superieure
  /*  for(x=0; x<L; x++){
      ix = x+L;
      if ((CelMvt[ix] & (MVT_E | MVT_O)) == MVT_E){
        CelMvt[ix] &= (-1) ^ MVT_E;
        CelMvt[ix] |= MVT_O;
      }
      if ((CelMvt[ix] & (MVT_EH | MVT_OH)) == MVT_EH){
        CelMvt[ix] &= (-1) ^ MVT_EH;
        CelMvt[ix] |= MVT_OH;
      }
      if ((CelMvt[ix] & (MVT_EB | MVT_OB)) == MVT_EB){
        CelMvt[ix] &= (-1) ^ MVT_EB;
        CelMvt[ix] |= MVT_OB;
      }
  }*/

  // forcage flux ouest-est sur la couche inferieure
  /*for(x=0; x<L; x++){
    ix = x+HL-2*L;
    if ((CelMvt[ix] & (MVT_E | MVT_O)) == MVT_O){
      CelMvt[ix] &= (-1) ^ MVT_O;
      CelMvt[ix] |= MVT_E;
    }
    if ((CelMvt[ix] & (MVT_EH | MVT_OH)) == MVT_OH){
      CelMvt[ix] &= (-1) ^ MVT_OH;
      CelMvt[ix] |= MVT_EH;
    }
    if ((CelMvt[ix] & (MVT_EB | MVT_OB)) == MVT_OB){
      CelMvt[ix] &= (-1) ^ MVT_OB;
      CelMvt[ix] |= MVT_EB;
    }
  }*/
#endif

#ifdef MVT_VFOR
  // forcage partiel du flux sur le bord gauche
  int cpt_vfor=0;

  for(z=0; z<CLNS; z++){
    for(x=0; x<1; x++){
      ix = z*CHL+5*CLEO+NB_MVT_EO+x;
      //for (y=5; y<CH-5; y++, ix+=CLEO){
      for (y=0; y<CH; y++, ix+=CLEO){
        //if (drand48()<1.0){
        //forcage
        if (drand48() < intensite_vfor){
        //if (1){
          if ((CelMvt[ix] & (MVT_E | MVT_O)) == MVT_O){
            UnsetMask(CelMvt[ix], MVT_O);
            SetMask(CelMvt[ix], MVT_E);
            cpt_vfor++;
          }

          if ((CelMvt[ix] & (MVT_EH | MVT_OH)) == MVT_OH){
            UnsetMask(CelMvt[ix], MVT_OH);
            SetMask(CelMvt[ix], MVT_EH);
            cpt_vfor++;
          }

          if ((CelMvt[ix] & (MVT_EB | MVT_OB)) == MVT_OB){
            UnsetMask(CelMvt[ix], MVT_OB);
            SetMask(CelMvt[ix], MVT_EB);
            cpt_vfor++;
          }
        }
        //if (!mvt_solid(ix)) CelMvt[ix] = MVT_E;
        //if (!mvt_solid(ix+1)) CelMvt[ix+1] = MVT_E;
        //if (!mvt_solid(ix+2)) CelMvt[ix+2] = MVT_E;
      }
    }
  }
#endif //MVT_VFOR

#ifdef MVT_GFOR
  // forcage global du flux
  int cpt_for = 0;
  int nb_mvt_for = nb_cel_fluide * intensite_gfor;

  //LogPrintf ("forcage sur %d cellules MVT\n", nb_mvt_for);

  while (nb_mvt_for > 0){
    ix = drand48()*CHLP;
    if ((CelMvt[ix] & (MVT_E | MVT_O)) == MVT_O){
      UnsetMask(CelMvt[ix], MVT_O);
      SetMask(CelMvt[ix], MVT_E);
      cpt_for++;
    }

    if ((CelMvt[ix] & (MVT_EH | MVT_OH)) == MVT_OH){
      UnsetMask(CelMvt[ix], MVT_OH);
      SetMask(CelMvt[ix], MVT_EH);
      cpt_for++;
    }

    if ((CelMvt[ix] & (MVT_EB | MVT_OB)) == MVT_OB){
      UnsetMask(CelMvt[ix], MVT_OB);
      SetMask(CelMvt[ix], MVT_EB);
      cpt_for++;
    }
    nb_mvt_for--;
  }
  //LogPrintf ("forcage de %d particules MVT\n", cpt_for);
#endif
}

//#ifdef BROUILLAGE_MVT
void do_mixing()
{
  int i, y, z, y1, y2;
  int cix, cix1, cix2, mvt;

  for(z=0; z<CLNS; z++){
    for (i=0; i<(CH/2); i++){

      //brouillage directionnel (bord gauche)
      y = 5+(CH-10)*drand48();
      cix = z*CHL+y*CLEO+NB_MVT_EO;
      mvt = CelMvt[cix] & (MVT_EB | MVT_EH);
      if (mvt == MVT_EB){
        UnsetMask(CelMvt[cix], MVT_EB);
        SetMask(CelMvt[cix], MVT_EH);
      }
      else if (mvt == MVT_EH){
        UnsetMask(CelMvt[cix], MVT_EH);
        SetMask(CelMvt[cix], MVT_EB);
      }

      //brouillage positionnel (bord gauche)
      y1 = 5+(CH-10)*drand48();
      y2 = 5+(CH-10)*drand48();
      cix1 = z*CHL+y1*CLEO+NB_MVT_EO;
      cix2 = z*CHL+y2*CLEO+NB_MVT_EO;
      if (!mvt_solid(cix1) && !mvt_solid(cix2)){
        //permutation
        mvt = CelMvt[cix1];
        CelMvt[cix1] = CelMvt[cix2];
        CelMvt[cix2] = mvt;
      }
    }
  }
}
//#endif

void compute_vel()
{
  static int cpt=0;
  //static int cpt_dump=0;
  int i, j, k, ix, iv, jv, imax, jmax, cpt1;
  int *pvx, *pvy;

  if (!Velx){
    AllocMemoryPrint("Velx", Velx, int, VHLP);
    ResetMemory(Velx, int, VHLP);
  }
  if (!Vely){
    AllocMemoryPrint("Vely", Vely, int, VHLP);
    ResetMemory(Vely, int, VHLP);
  }
  if (!Velx_time){
    if (VSTEP_TIME > 1){
      AllocMemoryPrint("Velx_time", Velx_time, int, VHLP*VSTEP_TIME);
      ResetMemory(Velx_time, int, VHLP*VSTEP_TIME);
    }
    else{
      Velx_time = Velx;
    }
  }
  if (!Vely_time){
    if (VSTEP_TIME > 1){
      AllocMemoryPrint("Vely_time", Vely_time, int, VHLP*VSTEP_TIME);
      ResetMemory(Vely_time, int, VHLP*VSTEP_TIME);
    }
    else{
      Vely_time = Vely;
    }
  }

  if (VSTEP_TIME > 1){
    if (!Velx_sum){
      AllocMemoryPrint("Velx_sum", Velx_sum, int, VHLP);
      ResetMemory(Velx_sum, int, VHLP);
    }
    if (!Vely_sum){
      AllocMemoryPrint("Vely_sum", Vely_sum, int, VHLP);
      ResetMemory(Vely_sum, int, VHLP);
      PrintTotalMemory();
    }
  }

#ifdef MVT_INTERP
  if (!Velx_interp){
    AllocMemoryPrint("Velx_interp", Velx_interp, int, HLP);
    ResetMemory(Velx_interp, int, HLP);
  }
  if (!Vely_interp){
    AllocMemoryPrint("Vely_interp", Vely_interp, int, HLP);
    ResetMemory(Vely_interp, int, HLP);
    PrintTotalMemory();
  }
#endif


  //while (CelMvt_lock) {LogPrintf("lock2\n"); };

  //CelMvt_lock = 1;
  //return;

  //moyennage dans l'espace
  //LogPrintf("compute_vel: moyennage dans l'espace\n");
  int nb_mvt_fluide;
  float fluide_ratio;
  pvx = Velx_time + VHLP*cpt;
  pvy = Vely_time + VHLP*cpt;
  for(k=0; k<CLNS; k++){ //profondeur
    for(jv=0; jv<VH; jv++){ //hauteur
      for(iv=0; iv<VL; iv++, pvx++, pvy++){ //largeur
        // moyenne des vitesses dans un rectangle de taille VSTEP_H x VSTEP_L
        ix = k*CHL + jv*CLEO*VSTEP_H + iv*VSTEP_L;
        if (jv < VH-1) jmax = VSTEP_H; else jmax = CH-jv*VSTEP_H;
        if (iv < VL-1) imax = VSTEP_L; else imax = CLEO-iv*VSTEP_L;
        *pvx = *pvy = 0;
        nb_mvt_fluide = 0;
        for(j=0; j<jmax; j++, ix+=CLEO-imax){ //hauteur
          for(i=0; i<imax; i++, ix++){ //largeur
            //if (CelMvt[ix] & (MVT_E | MVT_O | MVT_B | MVT_H)){
            if (CelMvt[ix] & MVT_SLOW){
              if (CelMvt[ix] & MVT_E) (*pvx)++;
              if (CelMvt[ix] & MVT_O) (*pvx)--;
              //if (CelMvt[ix] & MVT_O) (*pvx)++; //pression
              if (CelMvt[ix] & MVT_B) (*pvy)++;
              if (CelMvt[ix] & MVT_H) (*pvy)--;
              //if (CelMvt[ix] & MVT_H) (*pvy)++; //pression
            }
            if (CelMvt[ix] & MVT_FAST){
              if (CelMvt[ix] & MVT_EB) {(*pvx)++; (*pvy)++;}
              if (CelMvt[ix] & MVT_EH) {(*pvx)++; (*pvy)--;}
              //if (CelMvt[ix] & MVT_EH) {(*pvx)++; (*pvy)++;} //pression
              if (CelMvt[ix] & MVT_OB) {(*pvx)--; (*pvy)++;}
              //if (CelMvt[ix] & MVT_OB) {(*pvx)++; (*pvy)++;} //pression
              if (CelMvt[ix] & MVT_OH) {(*pvx)--; (*pvy)--;}
              //if (CelMvt[ix] & MVT_OH) {(*pvx)++; (*pvy)++;} //pression
            }
            if (!(CelMvt[ix] & MVT_SOLID)){
              nb_mvt_fluide++;
            }
          }
        }
        //moyennage sur les cellules fluides uniquement
        fluide_ratio = (float)nb_mvt_fluide/(imax*jmax);
        if (fluide_ratio < 0.2) fluide_ratio = 0.2;
        *pvx = (*pvx)/fluide_ratio;
        *pvy = (*pvy)/fluide_ratio;
      }
    }
  }

  //CelMvt_lock = 0;

  //moyennage dans le temps
  //LogPrintf("compute_vel: moyennage dans le temps\n");
  cpt1 = cpt+1;
  if (cpt1==VSTEP_TIME) cpt1 = 0;
  if (VSTEP_TIME > 1){
    for(i=0; i<VHLP; i++){
      Velx_sum[i] += Velx_time[cpt*VHLP+i] - Velx_time[cpt1*VHLP+i];
      Vely_sum[i] += Vely_time[cpt*VHLP+i] - Vely_time[cpt1*VHLP+i];
      Velx[i] = Velx_sum[i]; ///VSTEP_TIME;
      Vely[i] = Vely_sum[i]; ///VSTEP_TIME;
    }
    //if (Velx[0]) {LogPrintf("Velx[0]=%d\tcpt=%d\n", Velx[0], cpt); }
  }

  meanvel = 0;
  maxvel = 0;
  int vx, vy;
  float vel;
  for(i=0; i<VHLP; i++){
    vx = Velx[i];
    vy = Vely[i];
    vel = sqrt(vx*vx + vy*vy);
    meanvel += vel;
    if (maxvel<vel) maxvel = vel;
  }
  meanvel = meanvel/VHLP;

#ifdef MVT_INTERP
  //interpolation dans l'espace
  float fx0, fx1, fy0, fy1;
  float f00, f01, f10, f11;
  int Velx00, Velx01, Velx10, Velx11;
  int Vely00, Vely01, Vely10, Vely11;
  imax = VSTEP/NB_MVT_EO;
  //interpolation dans chaque plan de collisions
  //LogPrintf("interpolation dans chaque plan de collisions\n");
  for(k=0; k<CLNS; k++){ //profondeur
    for(jv=0; jv<VH-1; jv++){ //hauteur
      pvx = Velx+k*VHL+jv*VL;
      pvy = Vely+k*VHL+jv*VL;
      for(iv=0; iv<VL-2; iv++, pvx++, pvy++){ //largeur
        // interpolation des vitesses dans un rectangle de taille (imax)xVSTEP
        ix = (LN+k*DIST_MVT_NS)*HL + jv*L*VSTEP + iv*imax;// + (VSTEP>>1)*L + (VSTEP>>1);
        Velx00 = *pvx;
        Velx10 = *(pvx+VL);
        Vely00 = *pvy;
        Vely10 = *(pvy+VL);
        if (imax > 1){
          Velx01 = *(pvx+1);
          Velx11 = *(pvx+VL+1);
          Vely01 = *(pvy+1);
          Vely11 = *(pvy+VL+1);
        }
        for(j=0; j<VSTEP; j++, ix+=L-imax){ //hauteur
          fy1 = (float)j/VSTEP;
          fy0 = 1.0-fy1;
          if (imax > 1){ //interpolation entre les 4 sommets du rectangle englobant
            for(i=0; i<imax; i++, ix++){ //largeur
              fx1 = (float)i/imax;
              fx0 = 1.0-fx1;
              f00 = fy0*fx0;
              f01 = fy0*fx1;
              f10 = fy1*fx0;
              f11 = fy1*fx1;
              Velx_interp[ix] = Velx00*f00 + Velx01*f01 + Velx10*f10 + Velx11*f11;
              Vely_interp[ix] = Vely00*f00 + Vely01*f01 + Vely10*f10 + Vely11*f11;
            }
          }
          else{ //interpolation verticale entre 2 sommets
            Velx_interp[ix] = Velx00*fy0 + Velx10*fy1;
            Vely_interp[ix] = Vely00*fy0 + Vely10*fy1;
            ix+=imax;
          }
        }
      }
      // interpolation a cote du bord droit
      // on repete la derniere valeur interpolee
      // c'est plus facile et on reduit les effets de bord ...
      ix = (LN+k*DIST_MVT_NS)*HL + jv*L*VSTEP + iv*imax;// + (VSTEP>>1)*L + (VSTEP>>1);
      for(j=0; j<VSTEP; j++, ix+=L){ //hauteur
        for(i=0; i<L-iv*imax; i++){ //largeur
          Velx_interp[ix+i] = Velx_interp[ix-1];
          Vely_interp[ix+i] = Vely_interp[ix-1];
        }
      }
    }
  }
  
#if DIST_MVT_NS > 1
  compute_vel_interp_ns();
#endif
#endif


  cpt = cpt1;

  /*cpt_dump++;
  if (cpt_dump>=(H+L)/40) cpt_dump=0;
  if (!cpt_dump) dump_mvt(1);*/
  //LogPrintf("compute_vel: fin\n");
}

void compute_vel_interp_ns()
{
  int k, kv, i, ix;
  int ix0, ix1, kvmax;
  float fz0, fz1;
    
  //interpolation nord-sud dans les plans intermediaires
  //LogPrintf("interpolation nord-sud dans les plans intermediaires\n");

  for(k=0; k<CLNS; k++){ //profondeur
    kvmax = (k<CLNS-1) ? DIST_MVT_NS : LNS-k*DIST_MVT_NS;
    for(kv=1; kv<kvmax; kv++){ //profondeur intermediaire
      ix = (LN+k*DIST_MVT_NS+kv)*HL;
      ix0 = (LN+k*DIST_MVT_NS)*HL;
      ix1 = (k<CLNS-1) ? (LN+(k+1)*DIST_MVT_NS)*HL : ix0;
      fz1 = (float)kv/DIST_MVT_NS;
      fz0 = 1.0-fz1;
      for(i=0; i<HL; i++, ix++, ix0++, ix1++){
        Velx_interp[ix] = Velx_interp[ix0]*fz0 + Velx_interp[ix1]*fz1;
        Vely_interp[ix] = Vely_interp[ix0]*fz0 + Vely_interp[ix1]*fz1;
      }
    }
  }

}

void dump_mvt_in_out()
{
  static int cpt = 0;
  static int step = 0;
  FILE *fp;

  if (!cpt && !step){
    fp = fopen("MVT_IO.log","w");
    fprintf(fp,"      \t mvt_in    \t mvt_out   \t ratio\n");
    fclose(fp);
  }

  //fichier in-out
  if (step >= VSTEP_TIME){
    fp = fopen("MVT_IO.log","a");
    fprintf(fp,"%04d: \t%09ld \t%09ld \t%f\n", cpt, nb_mvt_in, nb_mvt_out, (nb_mvt_in) ? ((float)nb_mvt_out/nb_mvt_in) : 0);
    fclose(fp);
    step = nb_mvt_in = nb_mvt_out = 0;	//reset
    cpt++;
  }
  else
    step++;
}

void dump_densite()
{
  static int cpt = 0;
  FILE *fp;
  int ix, n;
  //int nb_cel_fluide_dump = 0;

  //recalcul nb_mvt, nb_mvt_sol et nb_cel_fluide
  int nb_mvt = 0; //nombre de particules mobiles
  int nb_mvt_sol = 0; //nombre de particules mobiles localisees dans une cellule solide
  for (ix=0; ix<CHLP; ix++){
    n = NbMvtCell[CelMvt[ix]];
    nb_mvt += n;
    if ((CelMvt[ix] & MVT_SOLID)) nb_mvt_sol += n;
    //if (!(CelMvt[ix] & MVT_SOLID)) nb_cel_fluide_dump++;
  }
  //densite = (float)nb_mvt / nb_cel_fluide;
  densite = (float)(nb_mvt - nb_mvt_sol) / nb_cel_fluide;

  //fichier densite
  fp = fopen("DENSI.log","a");
  if (!cpt) fprintf(fp,"\tmvt cells \tfluid nodes \tdensity \tsolid nodes\n");
  fprintf(fp,"%04d: \t%09ld \t%09ld \t%f \t%09ld\n", cpt++, (long)nb_mvt, nb_cel_fluide, densite, (long)nb_mvt_sol);
  //fprintf(fp,"%03d : \t%09ld \t%09ld (%09ld) \t%f\n", cpt++, (long)nb_mvt, nb_cel_fluide, nb_cel_fluide_dump, densite);
  fclose(fp);
}

void dump_vel()
{
  static int cpt = 0;
  static char start = 1;
  FILE *fp;

  if ((!Velx_interp) || (!Vely_interp)) {cpt++; return;}

#if 0
  //recalcul maxvel et meanvel
  int ix, vx, vy, vv, ncel;
  maxvel = 0.0;
  meanvel = 0.0;
  //float meanvel2 = 0.0;
  ncel = 0;
  for (ix=0; ix<HLP; ix++){
    vx = Velx_interp[ix];
    vy = Vely_interp[ix];
    vv = vx*vx+vy*vy;
    if (vv>maxvel) maxvel=vv;
    if (Phase[TE[ix].celltype] == FLUIDE) {meanvel += vv; /*meanvel2 += sqrt(vv);*/ ncel++;}
  }
  maxvel=sqrt(maxvel);
  meanvel=sqrt(meanvel/ncel); //approximation (en toute rigueur, il faut calculer la moyenne des racines carrees => plus de temps de calcul)
  //meanvel2/=ncel; //valeur exacte
#endif

  //fichier VEL.log
  fp = fopen("VEL.log","a");
  if (start){
    start = 0;
    fprintf(fp,"      \t maxvel \t meanvel\n");
  }
  fprintf(fp,"%04d: \t  %04d  \t  %04d\n", cpt++, (int)maxvel, (int)meanvel);
  //fprintf(fp,"%03d : \t  %04d  \t  %04d  \t  %04d\n", cpt++, (int)maxvel, (int)meanvel, (int)meanvel2);
  //fprintf(fp,"nb_mvt = %lu\n", nb_mvt);
  fclose(fp);
}

void dump_mvt(int cpt, int unit)
{
  static char filename[100];
  FILE *fp;

#if 1
  // binary  file
  //sprintf(filename,"HPP_%d_%d_%d.bin",CLNS,CH,CLEO);
  if (unit == UNIT_COMP)
    sprintf(filename,"HPP%04d.bin", cpt);
  else
    sprintf(filename,"HPP%05d_t0.bin", cpt);

  fp = fopen(filename,"w");
  if ( ! fp ){
    ErrPrintf("ERROR: cannot open binary file %s\n", filename);
    exit(-4);
  }


  //while (CelMvt_lock) {LogPrintf("lock3\n"); };

  //CelMvt_lock = 1;

  fwrite(CelMvt, sizeof(MvtField), CHLP, fp);

  //CelMvt_lock = 0;

  fclose(fp);
#endif

  //free(buf);
  /*if (inter){
    char command[100];
    sprintf(command, "gzip %s", nom);
    system(command);
  }*/

  // fichier vitesses
//#if (VSTEP > 1)
#if 1
  if (Velx && Vely){
    if (unit == UNIT_COMP)
      sprintf(filename,"HPP%04d_%d_%d_%d.vel", cpt, CLNS, VH, VL);
    else
      sprintf(filename,"HPP%05d_%d_%d_%d.vel", cpt, CLNS, VH, VL);

    fp = fopen(filename,"w");
    if ( ! fp ){
      ErrPrintf("ERROR: cannot open VEL file %s\n", filename);
      exit(-4);
    }

    //compute_vel();

    fwrite(Velx, sizeof(int), VHLP, fp);
    fwrite(Vely, sizeof(int), VHLP, fp);

    fclose(fp);
  }
#endif
}

void dump_collisions()
{
  FILE *fp;
  int mvt_in, mvt_out;

  fp = fopen("LGCA.log","w");

  //dump des regles de collision entre particules fluides
  fprintf(fp, "          IN            |          OUT\n");
  fprintf(fp, "E  W  B  T EB ET WB WT  |  E  W  B  T EB ET WB WT\n");
  fprintf(fp, "-------------------------------------------------\n");
  for(mvt_in=0; mvt_in<SIZE_MVT_FIELD; mvt_in+=4){
    mvt_out = Collisions[mvt_in];
    fprintf(fp, "%d  %d  %d  %d  %d  %d  %d  %d  |  %d  %d  %d  %d  %d  %d  %d  %d\n",
      (mvt_in & MVT_E)? 1 : 0, (mvt_in & MVT_O)? 1 : 0, (mvt_in & MVT_B)? 1 : 0, (mvt_in & MVT_H)? 1 : 0,
      (mvt_in & MVT_EB)? 1 : 0, (mvt_in & MVT_EH)? 1 : 0, (mvt_in & MVT_OB)? 1 : 0, (mvt_in & MVT_OH)? 1 : 0,
      (mvt_out & MVT_E)? 1 : 0, (mvt_out & MVT_O)? 1 : 0, (mvt_out & MVT_B)? 1 : 0, (mvt_out & MVT_H)? 1 : 0,
      (mvt_out & MVT_EB)? 1 : 0, (mvt_out & MVT_EH)? 1 : 0, (mvt_out & MVT_OB)? 1 : 0, (mvt_out & MVT_OH)? 1 : 0);
  }

#ifdef SOLSLOW
  //dump des regles de collision solide-fluide pour les particules lentes
  //selon la configuration du voisinage
  int slow[4] = {MVT_E, MVT_O, MVT_B, MVT_H};
  int mvt_sol_E, mvt_sol_O, mvt_sol_B, mvt_sol_H, mvt_sol, i;
  fprintf(fp, "\n   SOLIDE    |      IN      |      OUT\n");
  fprintf(fp, "E  O  B  H   | E  O  B  H   | E  O  B  H\n");
  for (mvt_sol=0; mvt_sol<16; mvt_sol++){
    mvt_sol_E = (mvt_sol & (MVT_E >> 2))? 1 : 0;
    mvt_sol_O = (mvt_sol & (MVT_O >> 2))? 1 : 0;
    mvt_sol_B = (mvt_sol & (MVT_B >> 2))? 1 : 0;
    mvt_sol_H = (mvt_sol & (MVT_H >> 2))? 1 : 0;
    for (i=0; i<4; i++){
      mvt_out = Collisions_SolSlow[mvt_sol][slow[i]];
      //mvt_in = MVT_SOLID | (1 << (2+i));
      //mvt_out = Collisions[mvt_in];
      fprintf(fp, "% d  %d  %d  %d  |  %d  %d  %d  %d  |  %d  %d  %d  %d\n",
              mvt_sol_E, mvt_sol_O, mvt_sol_B, mvt_sol_H,
              (i==0)? 1 : 0, (i==1)? 1 : 0, (i==2)? 1 : 0, (i==3)? 1 : 0,
              (mvt_out & MVT_E)? 1 : 0, (mvt_out & MVT_O)? 1 : 0, (mvt_out & MVT_B)? 1 : 0, (mvt_out & MVT_H)? 1 : 0);
    }
  }
#endif

#ifdef SOLFAST
  //dump des regles de collision solide-fluide pour les particules rapides (en diagonale)
  //selon la configuration du voisinage
  int fast[4] = {MVT_EB, MVT_EH, MVT_OB, MVT_OH};
  int mvt_sol_EB, mvt_sol_OB, mvt_sol_EH, mvt_sol_OH;//, mvt_sol, i;
  fprintf(fp, "\n   SOLIDE    |      IN      |      OUT\n");
  fprintf(fp, "EB EH OB OH  | EB EH OB OH  | EB EH OB OH\n");
  for (mvt_sol=0; mvt_sol<16; mvt_sol++){
    mvt_sol_EB = (mvt_sol & (MVT_EB >> 6))? 1 : 0;
    mvt_sol_EH = (mvt_sol & (MVT_EH >> 6))? 1 : 0;
    mvt_sol_OB = (mvt_sol & (MVT_OB >> 6))? 1 : 0;
    mvt_sol_OH = (mvt_sol & (MVT_OH >> 6))? 1 : 0;
    for (i=0; i<4; i++){
      mvt_out = Collisions_SolFast[mvt_sol][fast[i]];
      //mvt_in = MVT_SOLID | (1 << (6+i));
      //mvt_out = Collisions[mvt_in];
      fprintf(fp, "% d  %d  %d  %d  |  %d  %d  %d  %d  |  %d  %d  %d  %d\n",
              mvt_sol_EB, mvt_sol_EH, mvt_sol_OB, mvt_sol_OH,
              (i==0)? 1 : 0, (i==1)? 1 : 0, (i==2)? 1 : 0, (i==3)? 1 : 0,
              (mvt_out & MVT_EB)? 1 : 0, (mvt_out & MVT_EH)? 1 : 0, (mvt_out & MVT_OB)? 1 : 0, (mvt_out & MVT_OH)? 1 : 0);
    }
  }
#endif
  fclose(fp);
}



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

  //calcul de la signature
  sig=0;
  aux = (unsigned int*)CelMvt;
  for(i=0; i<CHLP*sizeof(MvtField)/sizeof(unsigned int); i++, aux++) sig += (*aux); //sig = sig ^ (*aux);

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

}
#endif

#endif //LGCA
