/* ReSCAL - Cells tracing
 *
 * 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.
 */


//mode TRACE_SRC : generation du bassin versant d'un point de sortie
//mode TRACE_AIRE : calcul de l'aire drainee (~debit) en chaque point (2D)
//mode TRACE_PAR : calcul du debit et de la pente en chaque point (2D), et sauvegarde du parcours effectif de chaque cellule
//mode TRACE_TRANS : comptage des transtions en chaque point (2D)
//mode TRACE3D_CEL : comptage des cellules d'un type donne qui transitent en chaque point (3D)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

#include "defs.h"
#include "macros.h"
#include "trace.h"
#include "cells.h"
#include "doublets.h"
#include "lgca.h"
#include "view.h"

extern int       H, L, P, HL, HLP;       // les dimensions de la terre
extern Cell  *TE;	                  // la 'terre'
extern double csp_time;                  // temps reel simule

int source_x, source_y, source_z, source_ix;   //point source
int src_flag = 0; //indicateur

void trace_src_init();
void trace_aire_init();
void trace_par_init();
void trace_trans_init();
void trace3d_cel_init();
void trace3d_time_init();

void trace_src_init_loop();
void trace_aire_init_loop();
void trace_par_init_loop();

void trace_src_end_loop();
void trace_aire_end_loop();
void trace_par_end_loop();

void trace_src_test(int ix);
void trace_aire_test(int ix);

int trace_src_point(int x, int z);
int trace_aire_point(int x, int z);

void trace_src_dump();
void trace_aire_dump();
void trace_par_dump();

void trace_dump();

#define NB_CEL_SRC 5 //nombre de cellules sources sur chaque colonne (injectees une par une)

void trace_init()
{
#ifdef TRACE_SRC
  trace_src_init();
#endif

#ifdef TRACE_AIRE
  trace_aire_init();
#endif

#ifdef TRACE_PAR
  trace_par_init();
#endif

#ifdef TRACE_PAR_COL
  trace_par_col_init();
#endif

#ifdef TRACE_TRANS
  trace_trans_init();
#endif

#ifdef TRACE3D_CEL
  trace3d_cel_init();
#endif
}

int trace_init_loop()
{
  static int x = 0, y = 1, z = 1, cpt_src = 0;
  int ix = 0;
  int res;

#ifdef TRACE_PAR
  // un seul point source
  if ((!x) || (cpt_src)) {
    x = L/2;
    z = L/2;
    LogPrintf("debut parcours x = %d   z = %d\n", x, z);
  }
  else{
    LogPrintf("fin parcours\n");
    return 1;
  }
#else
  // survol complet
  if (!cpt_src){
    if (!x) LogPrintf("NB_CEL_SRC = %d\n", NB_CEL_SRC);
    //on deplace le point source
    if (x<L-2){
      x++;
    }
    else if (z<P-2){
      x = 1;
      z++;
      LogPrintf("z = %d\n", z);
    }
    else{
      LogPrintf("fin trace\n");
      //view_dump_inter(0,IMG_PNG);
      dump_image(0,"png");
      trace_dump();
      return 1;
    }
  }
#endif

  cpt_src++;
  if (cpt_src >= NB_CEL_SRC) cpt_src = 0;

#ifdef MODELE_RIV
  //on pose une cellule source a la surface
  y = 1;
  ix = x + y*L + z*HL;
  while (TE[ix+L].celltype == EAU) {ix +=L; y++;}
  if (TE[ix].celltype == EAU){
    cellule_terre(BOUE, ix);
    //LogPrintf("point de depart : %d %d %d\n", x, y, z);
#if defined(TRACE_SRC)
    trace_src_test(ix);
#endif
  }
  else{
    ErrPrintf("ERROR: impossible to put source cell in (%d, %d)\n", x, z);
  }
#endif //MODELE_RIV

  source_x = x;
  source_y = y;
  source_z = z;
  source_ix = ix;

#ifdef TRACE_SRC
  trace_src_init_loop();
#endif

#ifdef TRACE_AIRE
  trace_aire_init_loop();
#endif

#ifdef TRACE_PAR
  trace_par_init_loop();
#endif

  return 0;
}

void trace_end_loop()
{
#ifdef TRACE_SRC
  trace_src_end_loop();
#endif

#ifdef TRACE_AIRE
  trace_aire_end_loop();
#endif

#ifdef TRACE_PAR
  trace_par_end_loop();
#endif
}

void trace_test(int ix)
{
#ifdef TRACE_SRC
  trace_src_test(ix);
#endif

#ifdef TRACE_AIRE
  trace_aire_test(ix);
#endif
}

int trace_point(int x, int z)
{
#ifdef TRACE_SRC
  if (trace_src_point(x, z)) return 1;
#endif

#ifdef TRACE_AIRE
  if (trace_aire_point(x, z)) return 2;
#endif

  return 0;
}

void trace_dump()
{
#ifdef TRACE_SRC
  trace_src_dump();
#endif

#ifdef TRACE_AIRE
  trace_aire_dump();
#endif

#ifdef TRACE_PAR
  trace_par_dump();
#endif
}

///////////////////////////////////////////////////////////////////////////////

#ifdef TRACE_SRC

#define DIST_SRC 100 //16 //5 //carre du rayon autour de la cellule test en mode TRACE_SRC

unsigned char *src_map;   //carte des points source
int test_x, test_z;


void trace_src_init()
{
  LogPrintf("init trace_src\n");
  LogPrintf("DIST_SRC = %d\n", DIST_SRC);

  AllocMemory(src_map, unsigned char, L*P);
  ResetMemory(src_map, unsigned char, L*P);
  test_x = 109; //149; //48; //109;//54; //161; //67;//422;//40;//36;//73;//1 + drand48()*(L-2);
  test_z = 190; //144; //190; //276;//121;//44;//63;//145;//1 + drand48()*(L-2);
  src_map[test_x+test_z*L] = NB_CEL_SRC;
  LogPrintf("test_x = %d\n", test_x);
  LogPrintf("test_z = %d\n", test_z);

  view_dump_init();
}

void trace_src_init_loop()
{
  src_flag = 0;
}

void trace_src_end_loop()
{
  if (src_flag) src_map[source_x+source_z*L]++;
}

void trace_src_test(int ix)
{
  if (TE[ix].celltype == BOUE){
    int x, y, z, dd;
    Calcule_xyz(ix, x, y, z);
    dd = (x - test_x)*(x - test_x) + (z - test_z)*(z - test_z);
    if (dd <= DIST_SRC){
    //if ((x == test_x) && (z == test_z)){
    //if (z == 1/*L-2*/){
      //src_map[source_x+source_z*L]++;
      src_flag = 1;
    }
    //trace[x+z*L] = 1;
  }
}

int trace_src_point(int x, int z)
{
  //return (src_map[x+z*L] == NB_CEL_SRC ) ?  1 : 0;
  return (src_map[x+z*L] > (NB_CEL_SRC >> 1)) ?  1 : 0;
  //return (src_map[x+z*L] > (NB_CEL_SRC >> 1)) || ((x == test_x) && (z == test_z)) ?  1 : 0;
}

void trace_src_dump()
{
  FILE *fp;
  int x,z;
  char nom[32];

  sprintf(nom, "BASSIN_%d.data", NB_CEL_SRC);

  fp = fopen(nom,"w");
  for (z=1; z<P-1; z++){
    for (x=1; x<L-1; x++){
      fprintf(fp, "%d\n", src_map[x+z*L]);
    }
  }
  fclose(fp);
}

#endif  //TRACE_SRC

///////////////////////////////////////////////////////////////////////////////

#ifdef TRACE_AIRE

unsigned char *trace_map;   //carte des points parcourus par une cellule source
int *aire_map;   //carte des aires drainees

void trace_aire_init()
{
  LogPrintf("init trace_aire\n");

  AllocMemory(trace_map, unsigned char, L*P);
  AllocMemory(aire_map, int, L*P);
  ResetMemory(aire_map, int, L*P);
  view_dump_init();
}

void trace_aire_init_loop()
{
  memset(trace_map, 0, L*P);
  trace_test(source_ix);
}

void trace_aire_end_loop()
{
  int i;
  for(i=0; i<L*P; i++) aire_map[i] += trace_map[i];

/*  int nn;
  for(i=0, nn=0; i<L*L; i++) nn += trace_map[i];
  LogPrintf("nn = %d\n", nn); */
}


void trace_aire_test(int ix)
{
  if (TE[ix].celltype == BOUE){
    int x, y, z;
    Calcule_xyz(ix, x, y, z);
    trace_map[x+z*L]=1;
    //LogPrintf("trace_test : %d %d %d\n", x, y, z);
  }
}

int trace_aire_point(int x, int z)
{
  return (trace_map) ? trace_map[x+z*L] : 0;
}

void trace_aire_dump()
{
  FILE *fp;
  int aire_max=0;
  int x, z, amx, amz;

  fp = fopen("AIRE.data","w");
  for(z=1; z<L-1; z++){
    for(x=1; x<L-1; x++){
      aire_map[x+z*L] /= NB_CEL_SRC;
      fprintf(fp, "%d\n", aire_map[x+z*L]);
      if (aire_map[x+z*L] > aire_max){
        aire_max = aire_map[x+z*L];
        amx = x;
        amz = z;
      }
    }
  }
  fclose(fp);

  LogPrintf("aire_max = %d en (%d,%d)\n", aire_max, amx, amz);

  /*
  fp = fopen("AIRE.car","wb");
  fwrite(aire_map, sizeof(int), L*L, fp);
  fclose(fp);
  */
}
#endif  //TRACE_AIRE

///////////////////////////////////////////////////////////////////////////////

#ifdef TRACE_TRANS

int *trace_diffusion; //carte des transitions de diffusion
int *trace_dissolution; //carte des transitions de dissolution
int *trace_cristallisation; //carte des transitions de cristallisation
int *trace_erosion; //carte des transitions d'erosion
int *trace_chk_erosion; //carte des transitions d'erosion bloquees
int *trace_transport_eo; //carte du transport est-ouest
int *trace_transport; //carte du transport
int *trace_gravi; //carte des transitions de gravite

void trace_trans_init()
{
  LogPrintf("init trace_trans\n");

#ifdef MODELE_DUN
/*
  //AllocMemory(trace_diffusion, int, L*P);
  AllocMemory(trace_erosion, int, L*P);
  AllocMemory(trace_chk_erosion, int, L*P);
  //ResetMemory(trace_diffusion, int, L*P);
  ResetMemory(trace_erosion, int, L*P);
  ResetMemory(trace_chk_erosion, int, L*P);
*/
  AllocMemory(trace_transport_eo, int, L*P);
  ResetMemory(trace_transport_eo, int, L*P);
#elif defined(MODELE_AVA)
  AllocMemory(trace_transport, int, L*P);
  AllocMemory(trace_gravi, int, L*P);
  ResetMemory(trace_transport, int, L*P);
  ResetMemory(trace_gravi, int, L*P);
#else
  AllocMemory(trace_diffusion, int, L*P);
  AllocMemory(trace_dissolution, int, L*P);
  AllocMemory(trace_cristallisation, int, L*P);
  ResetMemory(trace_diffusion, int, L*P);
  ResetMemory(trace_dissolution, int, L*P);
  ResetMemory(trace_cristallisation, int, L*P);
#endif
}

void trace_trans(int tr, int ix)
{
  int x, y, z, ind;
  Calcule_xyz(ix, x, y, z);
  ind = x + z*L;
#ifdef MODELE_DUN
  /*if (tr <= 3)
    trace_diffusion[ind]++;*/
  /*if (tr == 0)
    trace_erosion[ind]++;*/
  if (tr == 4)
    trace_transport_eo[ind]++;
#elif defined(MODELE_AVA)
  if (tr == 0)
    trace_gravi[ind]++;
  //else if (tr >= 3)
  else if ((tr >= 3) && (tr <=3))
    trace_transport[ind]++;
#else
  if (tr <= 3)
    trace_diffusion[ind]++;
  else if (tr <= 11)
    trace_dissolution[ind]++;
  else if (tr <= 19)
    trace_cristallisation[ind]++;
#endif
}

void trace_trans_blocked(int tr, int ix)
{
  int x, y, z, ind;
  Calcule_xyz(ix, x, y, z);
  ind = x + z*L;
#ifdef MODELE_DUN
  if (tr == 0)
    trace_chk_erosion[ind]++;
#endif
}

void trace_trans_dump()
{
  FILE *fp;
static int cpt = 0;
  char nom[128];

#ifdef MODELE_DUN
  /*sprintf(nom, "DIFFUS%03d.car", cpt++);
  fp = fopen(nom,"wb");
  fwrite(trace_diffusion, sizeof(int), L*P, fp);
  fclose(fp);*/
/*
  sprintf(nom, "EROSION%03d.car", cpt);
  fp = fopen(nom,"wb");
  fwrite(trace_erosion, sizeof(int), L*P, fp);
  fclose(fp);

  sprintf(nom, "EROSION_CHK%03d.car", cpt);
  fp = fopen(nom,"wb");
  fwrite(trace_chk_erosion, sizeof(int), L*P, fp);
  fclose(fp);
*/
  sprintf(nom, "TRANSPORT%03d.car", cpt);
  fp = fopen(nom,"wb");
  fwrite(trace_transport_eo, sizeof(int), L*P, fp);
  fclose(fp);

  /*
  //memset(trace_diffusion, 0, L*P*sizeof(int));
  memset(trace_erosion, 0, L*P*sizeof(int));
  memset(trace_chk_erosion, 0, L*P*sizeof(int));
  */
  memset(trace_transport_eo, 0, L*P*sizeof(int));

  cpt++;
#elif defined(MODELE_AVA)
  sprintf(nom, "TRANSPORT%03d.car", cpt);
  fp = fopen(nom,"wb");
  fwrite(trace_transport, sizeof(int), L*P, fp);
  fclose(fp);

  sprintf(nom, "GRAVI%03d.car", cpt);
  fp = fopen(nom,"wb");
  fwrite(trace_gravi, sizeof(int), L*P, fp);
  fclose(fp);

  //remise a zero
  memset(trace_transport, 0, L*P*sizeof(int));
  memset(trace_gravi, 0, L*P*sizeof(int));

  cpt++;
#else
  fp = fopen("DIFFUS.car","wb");
  fwrite(trace_diffusion, sizeof(int), L*P, fp);
  fclose(fp);

  fp = fopen("DISSOL.car","wb");
  fwrite(trace_dissolution, sizeof(int), L*P, fp);
  fclose(fp);

  fp = fopen("CRISTA.car","wb");
  fwrite(trace_cristallisation, sizeof(int), L*P, fp);
  fclose(fp);

  //remise a zero
  memset(trace_diffusion, 0, L*P*sizeof(int));
  memset(trace_dissolution, 0, L*P*sizeof(int));
  memset(trace_cristallisation, 0, L*P*sizeof(int));
#endif
}

void trace_trans_quit()
{
#ifdef MODELE_DUN
  //free(trace_diffusion);
  /*
  free(trace_erosion);
  free(trace_chk_erosion);
  */
  free(trace_transport_eo);
#elif defined(MODELE_AVA)
  free(trace_transport);
  free(trace_gravi);
#else
  free(trace_diffusion);
  free(trace_dissolution);
  free(trace_cristallisation);
#endif
}
#endif  //TRACE_TRANS

///////////////////////////////////////////////////////////////////////////////

#ifdef TRACE_PAR

//option
#define OPTI_PAR        //optimisation du parcours (un seul passage par chaque colonne)

//constantes
#define LG_PAR_MAX 10000  //longueur maximale du parcours sauvegarde
#define BUF_SIZE LG_PAR_MAX     //taille du tableau
#ifdef OPTI_PAR
# define INT_PENTE 10      //intervalle entre 2 positions pour le calcul de la pente
#else
# define INT_PENTE 100      //intervalle entre 2 positions pour le calcul de la pente
#endif
#define OFFSET_PENTE (INT_PENTE >> 1)

//structures
typedef struct un_element_de_parcours{
  int x,y,z;  //position
  float t;    //temps
  //float p;    //pente locale
} par_elt;

par_elt parcours[BUF_SIZE]; //portion de parcours d'une cellule test
#ifdef OPTI_PAR
int *par_map;               //indices du premier passage en chaque colonne
#endif
double *pente_map;          //estimation de la pente en chaque point
int *nb_pas;                //nombre de passages en chaque point
int cur_elt = 0;            //indice courant
int lg_par = 0;             //longueur du parcours


void elimine_cellules_actives()
{
  int i, j, k, n;
  Cell *aux = TE;
  // on enleve les cellules actives avant de calculer les differents parcours
  n=0;
  for(k=0; k < P; k++) // profondeur
    for(j=0; j < H; j++) // hauteur
      for(i=0; i < L; i++, aux++) //largeur
        if (aux->celltype == BOUE) {aux->celltype = EAU; n++;}

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

  LogPrintf ("cellules actives eliminees : %d\n", n);
}

void trace_par_init()
{
  LogPrintf("init trace_par\n");
  elimine_cellules_actives();
#ifdef OPTI_PAR
  AllocMemory(par_map, int, L*P);
  LogPrintf("optimisation du parcours : OPTI_PAR\n");
#endif
  AllocMemory(pente_map, double, L*P);
  AllocMemory(nb_pas, int, L*P);
  ResetMemory(pente_map, double, L*P);
  ResetMemory(nb_pas, int, L*P);
}

void trace_par_init_loop()
{
  //initialisation du parcours de la cellule test
  par_elt *pe = parcours;
  pe->x = source_x;
  pe->y = source_y;
  pe->z = source_z;
  pe->t = csp_time;
  //pe->p = 0.0;
#ifdef OPTI_PAR
  memset(par_map, -1, L*P*sizeof(int));
  par_map[source_x+source_z*L] = 0;
#endif
  lg_par = cur_elt = 1;
}

double calcul_pente(int i) //(par_elt *pe1, par_elt *pe2)
{
  double pente;

  /*if ((i-OFFSET_PENTE>=0) && (i+OFFSET_PENTE<cur_elt))*/{
    int i1, i2, dx, dy, dz, dh, dd;
    i1 = i-OFFSET_PENTE; if (i1<0) i1 = 0;
    i2 = i+OFFSET_PENTE; if (i2>=cur_elt) i2 = cur_elt-1;
    dx = parcours[i2].x - parcours[i1].x;
    dy = parcours[i2].y - parcours[i1].y;
    dz = parcours[i2].z - parcours[i1].z;

    //dh = dx*dx + dz*dz;
    dd = dx*dx + dy*dy + dz*dz;

    /*if (dh)
      pente = abs(dy) / sqrt(dh);
    else
      pente = (dy) ? 100.0 : 0.0;*/
    if (dd)
      pente = asin( /*abs*/(dy) / sqrt(dd) ) * (180.0 / 3.14159);
    else
      pente = 0.0;
  }
  //else pente = 0.0;

  return pente;
}

void dump_parcours()
{
  //sauvegarder une portion du parcours
  FILE *fp;
  par_elt *pe;
  float pente;
  static int num = 1;

  int i;
  fp = fopen("PARCOURS.data","a");
  fprintf(fp, "### PARCOURS %d ###\n", num++);
  for (i=0, pe = parcours; i<cur_elt; i++, pe++){
    pente = calcul_pente(i);  //((i-OFFSET_PENTE>=0) && (i+OFFSET_PENTE<cur_elt)) ? calcul_pente(pe-OFFSET_PENTE, pe+OFFSET_PENTE) : 0.0;
    fprintf(fp, "%03d %03d %03d %f %f\n", pe->x-1, pe->z-1, H-1-pe->y, pe->t, pente);
  }
  fclose(fp);
}

void trace_par(int ix)
{
  //memoriser la position de la cellule test
  if (TE[ix].celltype == BOUE){
    if (lg_par<LG_PAR_MAX){
      int x, y, z;
      par_elt *pe = parcours + cur_elt;
      Calcule_xyz(ix, x, y, z);
  #ifdef OPTI_PAR
      if (par_map[x+z*L] >= 0){
        //on est deja passe par la ...
        //donc on efface la fin du parcours ...
        int i, premier;
        premier = par_map[x+z*L];
        assert (cur_elt >= premier);
        for (i = cur_elt-1, pe = parcours+i; i>premier; i--, pe--)
          par_map[pe->x + pe->z*L] = -1;
        //et on recommence ! (a la nouvelle altitude)
        //lg_par -= (cur_elt - premier);
        cur_elt = premier;//+1;
        assert (pe == parcours+premier);
        assert ((pe->x == x) && (pe->z == z));
        //LogPrintf("premier=%d x=%d y=%d temps=%f\n",premier,x-1,z-1,temps);
      }
  #endif
      pe->x = x;
      pe->y = y;
      pe->z = z;
      pe->t = (float)csp_time;
      //pe->p = (cur_elt>=INT_PENTE) ? calcul_pente(pe, pe-INT_PENTE) : 0.0;
  #ifdef OPTI_PAR
      par_map[x+z*L] = cur_elt;
  #endif
      cur_elt++;
      lg_par++;
      if (cur_elt >= BUF_SIZE){
        //sauvegarder une portion du parcours
        ErrPrintf("WARNING: out of buffer parcours (size=%d) !\n", BUF_SIZE);
        dump_parcours();
        //trace_par_fin();
        cur_elt = 0;
      }
    }
    else{
      // parcours trop long, on remplace la cellule de boue par une cellule d'eau
      ErrPrintf("ERROR: cell course is too long (length=%d) !\n", lg_par);
      cellule_terre(EAU, ix);
      trace_par_end_loop();
     }
  }
}

void trace_par_end_loop()
{
  if (cur_elt){
    int i, n, ixz;
    par_elt *pe;
    double pente;
/*#ifdef OPTI_PAR
    LogPrintf("lg_par_opti = %d \t lg_par = %d\n", cur_elt, lg_par);
#else
    LogPrintf("lg_par = %d\n", lg_par);
#endif*/

    //sauvegarder le parcours complet
    dump_parcours();

    //calculer la pente moyenne a chaque pas
    for(i=0, pe=parcours; i<cur_elt; i++, pe++){
      ixz = pe->x + L*pe->z;
      n = nb_pas[ixz];
      pente = calcul_pente(i);
      pente_map[ixz] = (pente_map[ixz]*n + pente)/(n+1);
      nb_pas[ixz]++;
    }
    cur_elt = 0;
  }
}

void trace_par_dump()
{
  FILE *fp;
  int x,z;
  char nom[32];

#ifdef OPTI_PAR
  fp = fopen("DEBIT_OP.data","w");
#else
  fp = fopen("DEBIT.data","w");
#endif
  for (z=1; z<P-1; z++){
    for (x=1; x<L-1; x++){
      fprintf(fp, "%d\n", nb_pas[x+z*L]);
    }
    //fprintf(fp, "\n");
  }
  fclose(fp);

#ifdef OPTI_PAR
  sprintf(nom, "PENTE_OP_%d.data", INT_PENTE);
#else
  sprintf(nom, "PENTE_%d.data", INT_PENTE);
#endif
  fp = fopen(nom, "w");
  for (z=1; z<P-1; z++){
    for (x=1; x<L-1; x++){
      fprintf(fp, "%f\n", pente_map[x+z*L]);
    }
    //fprintf(fp, "\n");
  }
  fclose(fp);
}

#endif  //TRACE_PAR

///////////////////////////////////////////////////////////////////////////////

#ifdef TRACE_PAR_COL

//constants
#define NB_CEL_PAR_COL_MAX 100 //max number of colors
#define LG_PAR_COL_MAX 10000  //max length of each course
#define BUF_SIZE_COL 10 //LG_PAR_COL_MAX  //size of buffer for each course

extern float csp_angle;

int nb_par_col = 0; //number of colors

//structures
typedef struct un_element_de_parcours_color{
  int x,y,z;  //position
  int tr; //transition
  float t;    //time
} par_elt_col;

par_elt_col parcours[NB_CEL_PAR_COL_MAX][BUF_SIZE_COL]; //courses of all colored cells

int cur_elt[NB_CEL_PAR_COL_MAX];            //current indices for each course
int lg_par[NB_CEL_PAR_COL_MAX];             //length for each course


void trace_par_col_init()
{
  int i,j,ix,cm1;
  par_elt_col *pe;

  LogPrintf("init trace_par_col\n");

  for (cm1=0; cm1<NB_CEL_PAR_COL_MAX; cm1++){
    cur_elt[cm1] = -1; //default value
  }

  //initialization of all courses
  for (ix=0; ix<HLP; ix++){
    if (TE[ix].color){
      cm1 = TE[ix].color - 1;
      if (cm1 >= NB_CEL_PAR_COL_MAX){
        ErrPrintf("ERROR: trace_par_col_init - color=%d not allowed\n", TE[ix].color);
        exit(-1);
      }
      if (cm1 >= nb_par_col) nb_par_col = TE[ix].color;
      cur_elt[cm1]=0;
      //par_elt_col *pe = parcours[cm1];
      if (lg_par[cm1]==0)
        trace_par_col(ix,-1);
      else{
        ErrPrintf("ERROR: trace_par_col_init - two cells with the same color(%d), this is not allowed in TRACE_PAR_COL mode !\n", TE[ix].color);
        exit(-1);
      }
    }
  }

  LogPrintf("%d color(s) detected\n", nb_par_col)
}

void trace_par_col_rotation(float angle)
{
  int ix,col,cur;
  par_elt_col *pe;
  int x, y, z;

  //LogPrintf("trace_par_col_rotation - angle=%f\n",angle);

  //remove the colors
  for (ix=0; ix<HLP; ix++){
    TE[ix].color = 0;
  }

  //restore the colors
  for (col=1; col<=NB_CEL_PAR_COL_MAX; col++){
    cur = cur_elt[col-1];
    if (cur>=0){ //the color is present
      //LogPrintf("color %d is present\n", col);
      if (cur==0) cur = BUF_SIZE_COL;
      //read the last position
      pe = &parcours[col-1][cur-1];
      //rotation
      x = pe->x;
      y = pe->y;
      z = pe->z;
      //LogPrintf("x=%d y=%d z=%d\n",x,y,z);
      if (angle) Rotate_xz(x, z, angle);
      //LogPrintf("x=%d y=%d z=%d\n",x,y,z);
      //new position
      ix = z*HL + y*L + x; 
      //update the color
      if ((x>=1) && (x<L-1) && (y>=1) && (y<H-1) && (z>=1) && (z<P-1)) {
        if (TE[ix].color == 0) {
#ifdef MODELE_DUN
          if ((TE[ix].celltype == GR) || (TE[ix].celltype == GRJ)){ 
            TE[ix].color = col;
          }
          else if (TE[ix].celltype == EAUC ){
            LogPrintf("WARNING: trace_par_col_rotation - not a grain (%d), GR added for color %d\n", TE[ix].celltype, col); /*exit(-1);*/
            TE[ix].celltype = GR; //adding a grain !
            TE[ix].color = col;
          }
          else{
            ErrPrintf("WARNING: trace_par_col_rotation - not a grain (%d), trace is stopped for color %d\n", TE[ix].celltype, col); /*exit(-1);*/
            cur_elt[col-1]=-1; //trace is stopped for this grain
          }
#else
          TE[ix].color = col;
#endif
        }
        else{
          ErrPrintf("WARNING: trace_par_col_rotation - cell in %d is already colored (%d), cannot assign color %d\n", ix, TE[ix].color,col);
          cur_elt[col-1]=-1; //trace is stopped for this grain
        }
      }
      else{
        ErrPrintf("WARNING: trace_par_col_rotation - colored (%d) grain is going outside\n", col);
        cur_elt[col-1]=-1; //trace is stopped for this grain
      }
    }
  } 
}

void dump_par_col(int col)
{
  FILE *fp;
  par_elt_col *pe;
  //static int num = 1;
  char name[100];
  int j;

  //save course for one color from current buffer
  sprintf(name,"PAR_COL%03d.data", col);
  fp = fopen(name,"a");
  for (j=0, pe = parcours[col-1]; j<cur_elt[col-1]; j++, pe++){
    fprintf(fp, "%05d %05d %05d %03d %e\n", pe->x-1, pe->z-1, H-1-pe->y, pe->tr, pe->t);
  }
  fclose(fp);
}

void trace_par_col(int ix, int tr)
{
  if (TE[ix].color){
    int cm1 = TE[ix].color-1;
    //LogPrintf("trace_par_col: ix=%d, tr=%d, col=%d, csp_angle=%f\n", ix, tr, TE[ix].color, csp_angle);
    if (cur_elt[cm1]<0) return;
    if (lg_par[cm1]<LG_PAR_COL_MAX){
      int x, y, z;
      par_elt_col *pe = parcours[cm1] + cur_elt[cm1];
      Calcule_xyz(ix, x, y, z);
      //LogPrintf("trace_par_col: x=%d, y=%d, z=%d\n", x, y, z);
      if (csp_angle) Rotate_xz(x, z, -csp_angle);
      pe->x = x;
      pe->y = y;
      pe->z = z;
      pe->tr = tr;
      pe->t = (float)csp_time;
      cur_elt[cm1]++;
      lg_par[cm1]++;
      if (cur_elt[cm1] >= BUF_SIZE_COL){
        //save course
        //ErrPrintf("WARNING: out of buffer par_col (col=%d, size=%d) !\n", TE[ix].color, BUF_SIZE_COL);
        dump_par_col(TE[ix].color);
        cur_elt[cm1] = 0;
      }
    }
    else{
      //course too long
      ErrPrintf("ERROR: cell course is too long (col=%d, length=%d) !\n", TE[ix].color, lg_par[cm1]);
      //color is removed
      TE[ix].color = 0;
    }
  }
}


#endif  //TRACE_PAR_COL

///////////////////////////////////////////////////////////////////////////////

#ifdef TRACE3D_CEL

//constantes
#ifdef MODELE_DUN
//#define CEL_TEST  EAUT
#define CEL_TEST  GRJ
#endif

#ifdef LGCA
//#define TRACE_MVT
#endif

int *cpt_cel; //compteur de cellules (en chaque point)
#ifdef TRACE_MVT
int *cpt_mvt_est; //compteur de cellules mvt_est (en chaque point)
int *cpt_mvt_ouest; //compteur de cellules mvt_ouest (en chaque point)
int *cpt_check_mvt_est; //compteur de tests reussis pour le controle du transport vers l'est (en chaque point)
int *cpt_check_mvt_bas; //compteur de tests reussis pour le controle du transport vers le bas (en chaque point)
#endif //TRACE_MVT
int nb_trace3d_cel=0; // nombre d'appels de la procedure trace3d_cel()

void trace3d_cel_init()
{
  LogPrintf("init trace3d_cel\n");

  AllocMemory(cpt_cel, int, HLP);
  ResetMemory(cpt_cel, int, HLP);

#ifdef TRACE_MVT
  AllocMemory(cpt_mvt_est, int, HLP);
  ResetMemory(cpt_mvt_est, int, HLP);

  AllocMemory(cpt_mvt_ouest, int, HLP);
  ResetMemory(cpt_mvt_ouest, int, HLP);

  AllocMemory(cpt_check_mvt_est, int, HLP);
  ResetMemory(cpt_check_mvt_est, int, HLP);

  AllocMemory(cpt_check_mvt_bas, int, HLP);
  ResetMemory(cpt_check_mvt_bas, int, HLP);
#endif //TRACE_MVT
}

void trace3d_cel()
{
  int i;
  for (i=0; i<HLP; i++){
    if (TE[i].celltype == CEL_TEST) cpt_cel[i]++;
#ifdef TRACE_MVT
    if (mvt_est(i)) cpt_mvt_est[i]++;
    if (mvt_ouest(i)) cpt_mvt_ouest[i]++;
    if (check_mvt_est(i)) cpt_check_mvt_est[i]++;
    if (check_mvt_bas(i)) cpt_check_mvt_bas[i]++;
#endif //TRACE_MVT
  }

  nb_trace3d_cel++;
}

void trace3d_cel_dump()
{
  FILE *fp;
  static int cpt = 0;
  //char nom[128];

  //sprintf(nom, "CEL%03d.bin", cpt++);
  fp = fopen("CEL.bin","wb");
  fwrite(cpt_cel, sizeof(int), HLP, fp);
  fclose(fp);

#ifdef TRACE_MVT
  fp = fopen("MVT_EST.bin","wb");
  fwrite(cpt_mvt_est, sizeof(int), HLP, fp);
  fclose(fp);

  fp = fopen("MVT_OUEST.bin","wb");
  fwrite(cpt_mvt_ouest, sizeof(int), HLP, fp);
  fclose(fp);

  fp = fopen("CHK_MVT_EST.bin","wb");
  fwrite(cpt_check_mvt_est, sizeof(int), HLP, fp);
  fclose(fp);

  fp = fopen("CHK_MVT_BAS.bin","wb");
  fwrite(cpt_check_mvt_bas, sizeof(int), HLP, fp);
  fclose(fp);
#endif //TRACE_MVT

  fp = fopen("TRACE3D_CEL.data","w");
  fprintf(fp, "%d\n", nb_trace3d_cel);
  fclose(fp);

  //memset(cpt_cel, 0, HLP*sizeof(int));
  cpt++;
}

void trace3d_cel_quit()
{
  free(cpt_cel);
#ifdef TRACE_MVT
  free(cpt_mvt_est);
  free(cpt_mvt_ouest);
  free(cpt_check_mvt_est);
  free(cpt_check_mvt_bas);
#endif //TRACE_MVT
}

#endif //TRACE3D_CEL

