/* ReSCAL - Render and display module
 *
 * 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 <string.h>
#include <assert.h>
#include <math.h>

#include "defs.h"
#include "macros.h"
#include "view.h"
#include "cells.h"
#include "surface.h"
#include "trace.h"
#include "transitions.h"
#include "lgca.h"

extern int prog;    //main prog
extern int H, L, P, HL, HLP;       // les dimensions de la terre
extern int LEO, LN, LS, LNS, HLN;    //largeur est-ouest, limite nord, limite sud, largeur nord-sud
extern Cell  *TE;	           // la 'terre'
#ifdef ALTI
extern short *alti; // elevations locales du terrain
extern Vec2 *norm2d;
extern Vec3 *norm3d;          // normale 3d a la surface definie par alti
#endif
extern int *psol;                    // indicateur des plans solides verticaux est-ouest
extern int h_floor;
#ifdef LGCA
extern int VH, VL, VHL, CLNS ;	// dimensions du tableau des vitesses
extern int *Velx, *Vely;	//vitesses locales moyennees
extern int *Velx_sum, *Vely_sum;	//vitesses locales additives
extern int *Velx_interp, *Vely_interp;	//vitesses locales interpolees
extern int lgca_ready;
#endif

unsigned char opt_h=0;
unsigned char opt_cv=0;
unsigned char opt_lc=0;
unsigned char opt_tr=0;
unsigned char opt_al=0;
unsigned char opt_lal=0;
unsigned char opt_vel=0;
unsigned char opt_vss=0;
unsigned char opt_ls=0;
unsigned char opt_col=0;
unsigned char hh, typtr, typal;
int abs_cv=0, prof_cv=0;
unsigned char reorient_flag=0;
unsigned char horizontal_display=0;

unsigned char *color_alti = NULL; //affectation des couleurs selon l'altitude
unsigned char tr[MAX_CELL];// = {0 ,0 ,0 ,0 ,0 ,0,};  //transparence des cellules par type
unsigned char *image=NULL;
unsigned char *image_end=NULL;
int img_w, img_h;

//direction of the incident light source
float light_elevation = 45; // in degrees, from the horizon
float light_azimuth = -90; //in degrees, clockwise from north
Vec3 light;

//display mode of the current orientation
int vdir_mode = VDIR_NONE;

//intervalle pour les pointilles
int ilc = 2;

void set_light_xyz();

void bad_options(char *opt)
{
  ErrPrintf("erreur lecture option : %s\n", opt);

  exit(-4);
}

void show_view_options()
{
  printf("\nGRAPHICAL OPTIONS\n");
  printf("  -h <val> \t height of the horizontal cross section (H/2 by default)\n");
  printf("  -cv0 \t\t centered vertical cross sections\n");
  printf("  -cv <a>x<o> \t vertical cross sections at abscissa <a> and ordinate <o>\n");
  //printf("  -diag \t diagonal layers of cells \n");
  printf("  -lc \t\t show cut lines\n");
  printf("  -al <cel> \t height map of the cells of type <cel>\n");
  printf("  -lal <cel> \t low height map of the cells of type <cel>\n");
  printf("  -tr <cel> \t transparent cells of type <cel>\n");
//#if defined(MODELE_DUN) || defined(MODELE_RIV) || defined(MODELE_AVA)
#ifdef ALTI
  printf("  -ls <angle> \t surface rendering with ligth source elevation = <angle> degrees\n");
#endif
#ifdef LGCA
  if (prog == PROG_RESCAL)
    printf("  -vel \t\t view flow mean velocity on a vertical plan\n");
    printf("  -vel2 \t view flow mean velocity and vector field on a vertical plan\n");
#if defined(CGV)
  if (prog == PROG_RESCAL)
    printf("  -vss \t\t view shear stress\n");
#endif
#endif
#ifdef CELL_COLOR
  printf("  -vcol \t view colored cells\n");
#endif
  fflush(stdout);
}

int view_init(int argc, char *argv[], int *pL, int *pH)
{
  int i;

  hh = H>>1;
  memset(tr, 0, MAX_CELL);

  for(i=1; i<argc; i++){
    if (!strcmp(argv[i],"-cv0")){
      opt_cv = 1;
      abs_cv = L>>1;
      prof_cv = P>>1;
    }
    else if (!strcmp(argv[i],"-cv")){
      char c;
      int n=0;
      opt_cv = 1;
      if (i+1<argc) n = sscanf(argv[++i], "%d%c%d", &abs_cv, &c, &prof_cv);
      //LogPrintf("n=%d\n",n);
      if (n<3) bad_options("-cv <abs>x<prof>");
      abs_cv++;
      prof_cv++;
      //LogPrintf("abs_cv=%d   prof_cv=%d\n", abs_cv, prof_cv);
    }
    else if (!strcmp(argv[i],"-lc"))
      opt_lc = 1;
    else if (!strcmp(argv[i],"-h")){
      opt_h = 1;
      hh = atoi(argv[++i]);
      if (hh>H-2) hh=H-2;
    }
    else if (!strcmp(argv[i],"-tr")){
      opt_tr = 1;
      typtr = atoi(argv[++i]);
      tr[typtr] = 1;
    }
    else if (!strcmp(argv[i],"-al")){
      opt_al = 1;
      typal = atoi(argv[++i]);
    }
    else if (!strcmp(argv[i],"-lal")){
      opt_lal = 1;
      typal = atoi(argv[++i]);
    }
#ifdef ALTI
    else if (!strcmp(argv[i],"-ls")){
      opt_ls = 1;
      light_elevation = atof(argv[++i]);
      set_light_xyz();
    }
#endif
#ifdef LGCA
    else if ((prog == PROG_RESCAL) && (!strcmp(argv[i],"-vel")))
      opt_vel = 1;
    //LogPrintf("opt_vel = %d\n", opt_vel);
    else if ((prog == PROG_RESCAL) && (!strcmp(argv[i],"-vel2")))
      opt_vel = 2;
    //LogPrintf("opt_vel = %d\n", opt_vel);
#if defined(CGV)
    else if ((prog == PROG_RESCAL) && (!strcmp(argv[i],"-vss")))
      opt_vss = 1;
#endif
#endif
#ifdef CELL_COLOR
    else if (!strcmp(argv[i],"-vcol"))
      opt_col = 1;
#endif

    if (!opt_al && !opt_lal && !opt_tr && !opt_ls) opt_h = 1;
  }

  if (LNS == 1) opt_lc = opt_vss = 0;

  if (opt_vel && !opt_cv) prof_cv = P>>1;

  if (opt_cv || opt_vel){
    LogPrintf("abs_cv = %d\n", abs_cv);
    LogPrintf("prof_cv = %d\n", prof_cv);
  }

  img_w = (L<=IMG_W_MAX)? L : IMG_W_MAX;
  if (L>IMG_W_MAX) {LogPrintf("largeur de l'image fixee a IMG_W_MAX = %d\n", IMG_W_MAX); }
  img_h = P;//LNS+2;//L;
  if (opt_cv){
    if (LNS>1) img_w += H;
    img_h += H; //(LNS>1) ? 2*H : H;
  }
  if (opt_vel && (!opt_vss || !opt_cv || (L>P))) img_h += H;
#if defined(CGV)
  if (opt_vss){
    horizontal_display = (L<=P);
    if (horizontal_display)
      img_w += L;
    else
      img_h += LNS+2;
  }
#endif
  *pL = img_w;
  *pH = img_h;

  if (!image){
    /*image = (unsigned char *) malloc(img_w*img_h);
  	if (!image){
          ErrPrintf("erreur malloc image (img_w = %d   img_h = %d)\n", img_w, img_h);
		  exit(-4);
	  }*/
    AllocMemoryPrint("image",image, unsigned char, img_w*img_h);
    PrintTotalMemory();
    memset(image, BORD, img_w*img_h);
    image_end = image + img_w*img_h;
  }

  return 0;
}


void view_palette(int *colors)
{
  int i;

#if defined(MODELE_DUN )

  colors[GR] = 0x00000000; //grain stable
  colors[GRJ] = 0x00ff0000; //grain instable jeune
  colors[GRV] = 0x00ff00ff; //grain instable vieux
  colors[EAUC] = 0x0000ffff; //eau
  colors[EAUT] = 0x000000ff; //turbulence
  colors[BORD] = 0x007f7f7f; //bord
  colors[DUM] = 0x00ffff00; //inerte
  colors[IN] = 0x007f7f00; //source de grains
  colors[OUT] = 0x009f0000; //sortie de grains
  colors[TUNNEL] = 0x00ffffff; //0x00009f00; //grain colore

#elif defined(MODELE_AVA)

  colors[0] = 0x00000000; //grain
  colors[1] = 0x0000ffff; //air
  colors[2] = 0x007f7f7f; //bord
  colors[3] = 0x00ffff00; //inerte
  colors[4] = 0x007f7f00; //source de grains
  colors[5] = 0x009f0000; //sortie de grains
  colors[9] = 0x00ffffff; //0x00009f00; //grain colore

#elif defined(MODELE_RIV)

  colors[0] = 0x00000000; //terre
  colors[1] = 0x00ffffff; //boue
  colors[2] = 0x00ff0000; //eau
  colors[3] = 0x0000ffff; //boue turbulente
  colors[4] = 0x0000c8c8; //pierre
  colors[5] = 0x007f7f7f; //bord
  colors[6] = 0x000000ff; //inerte
  colors[7] = 0x007f7f00; //source de grains
  colors[8] = 0x009f0000; //sortie de grains

#elif defined(MODELE_LIFE)

  colors[DEAD] = 0x00ffffff;
  colors[ALIVE] = 0x00000000;
  colors[BORD] = 0x007f7f7f;
  colors[DUM] = 0x000000ff;

#else
  colors[0] = 0x00000000; //PLUS
  colors[1] = 0x00ffffff; //ZERO
  colors[2] = 0x00ff0000; //MOINS
  colors[BORD] = 0x007f7f7f;
  colors[DUM] = 0x000000ff;
  colors[5] = 0x007f7f00; //IN
#ifdef MODELE_CHE
  colors[6] = 0x00007fff; //AIR
#else
  colors[6] = 0x00afafaf; //BT
#endif
  colors[7] = 0x0000c8c8; //PIERRE

#endif

#ifdef PARALLEL
  colors[TUNNEL] = 0x007f007f;
#endif

#ifdef MODELE_CMB
  colors[COUPE] = 0x00007f00;
#else
  colors[COUPE] = 0x00000000;
#endif

  //assert(H <= 256);
  /*for(i=0; i<H-2; i++) {
    red = 255*i/(H-3);
    green = 195*i/(H-3);
    blue = 127*i/(H-3);
    colors[255+i-(H-3)] = (red << 16) | (green << 8) | blue;
  }*/
  //if (opt_al || opt_lal){
  if (1){
    unsigned char red, green, blue;
    float coef;
    int nb_color_alti = 255 - MAX_CELL;
    for(i=0; i<nb_color_alti; i++) {
      coef = (float)(nb_color_alti-1-i)/(nb_color_alti-1);
      red = 255*coef + 64*(1-coef);
      green = 195*coef + 0*(1-coef)*195/255;
      blue = 127*coef + 0*(1-coef)*127/255;
      colors[255-i] = (red << 16) | (green << 8) | blue;
    }
    if (!color_alti){
      AllocMemory(color_alti, unsigned char, H);
      for(i=0; i<H; i++){
        color_alti[i] = 255 - i*(nb_color_alti-1)/(H-1);
        //LogPrintf("color_alti[%d] = %d\n", i, color_alti[i]);
      }
    }
  }
}

// draw line
void view_line(unsigned char *pix0, int dx, int dy, unsigned char color)
{
  int sdx,sdy,px,py,dxabs,dyabs,i;
  float slope;
  unsigned char *pix;

  /*if (dx*dx+dy*dy > 0){
    LogPrintf("view_line() dx = %d, dy = %d\n", dx, dy);
  }*/

  //horizontal clipping
  int ipix = (int)(pix0 - image);
  int xpix = ipix % img_w;
  Check_min_max(dx, -xpix, img_w-1-xpix);

  dxabs = abs(dx);
  dyabs = abs(dy);
  sdx = sgn(dx);
  sdy = sgn(dy);
  i = 0;
  if (dxabs >= dyabs) /* the line is more horizontal than vertical */
  {
    slope = (float)dy / (float)dx;
    while(1)
    {
      py = slope*i;
      pix = pix0+i+py*img_w;
      if ((pix >= image) && (pix < image_end)) *pix = color;
      if (i==dx) break;
      i+=sdx;
    }
  }
  else /* the line is more vertical than horizontal */
  {
    slope = (float)dx / (float)dy;
    while(1)
    {
      px = slope*i;
      pix = pix0+px+i*img_w;
      if ((pix >= image) && (pix < image_end)) *pix = color;
      if (i==dy) break;
      i+=sdy;
    }
  }
}

// draw filled rectangle
void view_rect(unsigned char *pix0, int width, int height, unsigned char color)
{
  unsigned char *pix;
  int y;

  assert(pix0 + (height-1)*img_w + width <= image_end);

  pix = pix0;

  for (y=0; y < height; y++, pix+=img_w)
    memset(pix, color, width);
}

void set_light_xyz()
{
  float alpha = light_azimuth*PI/180;
  float phi = light_elevation*PI/180;
  light.x = cosf(phi)*sinf(alpha); //1.0;
  light.y = sinf(phi); //0.0;
  light.z = -cosf(phi)*cosf(alpha);
}

void rotate_light(float angle)
{
  light_azimuth += angle;
  //LogPrintf("rotate_light: angle = %f   azimuth = %f\n", angle, light_azimuth);
  set_light_xyz();
}

unsigned char * view()
{
  Cell *pt;
  unsigned char *pix, *cur_image;
  int i,j,k,imax;

  pix = cur_image = image;

  if (opt_al){
    // vue horizontale en mode altimap
    assert (typal != BORD);
    for(j=1; j < P-1; j++){
      pt = TE + j*HL;// + HLN;
      imax = Min(L,img_w);
      pix = cur_image + j*img_w;
      for(i=0; i < imax; i++, pt++){
        int k=1; //1;
//#ifdef TRACE_SRC
#if defined(TRACE_SRC) || defined(TRACE_AIRE)
        if (trace_point(i,j)){
          *pix++ = COUPE;
        }
        else
#endif
        if (!opt_tr){
          //if ((pt+2*L)->celltype == DUM) //couloir
          //if ((pt+5*L)->celltype == DUM) //couloir
          if (psol[j])
              *pix++ = DUM;
          else{
            while ((k < H) && ((pt+k*L)->celltype != typal)) k++;
            //while ((k < H) && ((pt+k*L)->celltype != typal) && ((pt+k*L)->celltype != DUM)) k++;
#ifdef CELL_COLOR
            *pix++ = (k < H)? (opt_col && ((pt+k*L)->color == 1)) ? TUNNEL : color_alti[k] : BORD;
#else
            *pix++ = (k < H)? color_alti[k] : BORD;
#endif
            //*pix++ = ((pt+L)->celltype == DUM)? DUM : (k < H)? color_alti[k] : BORD;
          }
        }
        else{
          while ((k < H) && tr[(pt+k*L)->celltype]) k++;
          *pix++ = (k < H) ? ((pt+k*L)->celltype == typal) ? color_alti[k]/*256-k*/ : (pt+k*L)->celltype : BORD;
          //while ((k < H) && ((pt+k*L)->celltype != typal) && ((pt+k*L)->celltype != typf)) k++;
          //*pix++ = (k < H)? ((pt+k*L)->celltype != typf)? 255-k : typf : BORD;
        }
      }
    }
  }
  else if (opt_lal){
    // vue horizontale en mode altimap basse
    assert (typal != BORD);
    for(j=1; j < P-1; j++){
      pt = TE + j*HL + 1;
      imax = Min(L-2,img_w-2);
      pix = cur_image + j*img_w + 1;
      for(i=0; i < imax; i++, pt++){
        int k=H-1-h_floor;
        while ((k > 0) && ((pt+(k-1)*L)->celltype == typal)) k--;
        *pix++ = (k < H-1-h_floor)? color_alti[k] : BORD;
        if (opt_tr){
          while ((k > 0) && tr[(pt+k*L)->celltype]) k--;
          if (k > 0) *(pix-1) = (pt+k*L)->celltype;
        }
      }
    }
  }
#ifdef ALTI
  else if (opt_ls){
  //affichage de la surface avec eclairage
    if (norm3d){
      int vcol;
      int kmax;
      //float *ptnx, *ptny, *ptnz;
#ifdef NORM2D
      Vec2 *ptn;
#else
      Vec3 *ptn;
#endif
      float ldiff = 0.0;
      //pix = cur_image+img_w;
      for(j=1; j < P-1; j++){
#ifdef NORM2D
        ptn = norm2d + (j-LN)*LEO;

#else
        ptn = norm3d + (j-LN)*LEO;

#endif
        pix = cur_image + j*img_w + 1;
        //imax = (img_w == L)? L-2 : img_w-2;
        imax = Min(L-2,img_w-1);
        for(i=0; i < imax; i++, ptn++, pix++){
          if (psol[j])
              *pix = DUM;
          else{
#ifdef NORM2D
            ldiff = (ptn->x)*light.x + (ptn->y)*light.y; //lumiere diffuse (2d)
#else
            ldiff = (ptn->x)*light.x + (ptn->y)*light.y + (ptn->z)*light.z; //lumiere diffuse
#endif
            if (ldiff < 0.0) ldiff = 0.0;
            vcol = ldiff*(H-1);
            *pix = color_alti[H-1-vcol];
            if (opt_tr){
              k = 1;
              pt = TE+HLN+j*HL+i+1;
              kmax = H-alti[j*LEO+i];
              while ((k <= kmax) && tr[(pt+k*L)->celltype]) k++;
              if (k < kmax) *pix = (pt+k*L)->celltype;
            }
          }
        }
        //pix++;
      }
    }
  }
#endif //ALTI
  else if (opt_tr){
  // vue horizontale de toutes les cellules avec transparence
    for(j=1; j < P-1; j++){
      pt = TE + j*HL;
      pix = cur_image + j*img_w;
      imax = Min(L,img_w);
      for(i=0; i < imax; i++, pt++){
        int k=1;
        while ((k < H) && tr[(pt+k*L)->celltype]) k++;
        *pix++ = (k < H) ? (pt+k*L)->celltype : BORD;
      }
    }
  }
/*#ifdef CELL_COLOR
  else if (opt_col){
    if (alti){
      pt = TE+HLN;
      for(j=0; j < LNS; j++, pt += (H-1)*L){
        pix++; pt++;
        for(i=0; i < L-2; i++, pix++, pt++){
          k = H-alti[j*(L-2)+i];
          *pix = ((pt+k*L)->color == 1) ? TUNNEL : color_alti[k];
        }
        pix++; pt++;
      }
    }
  }
#endif*/
  else{
    //coupe horizontale (au milieu par defaut)
    //extern float *grdv;
    //int val, coef = VSTEP_H*VSTEP_L*VSTEP_TIME/50;
    //if (grdv)
    for(j=1; j < P-1; j++){
      pt = TE + j*HL + L*(H-1-hh);
      pix = cur_image + j*img_w;
      for(i=0; i < img_w; i++, pt++){
#ifdef CELL_COLOR
        *pix++ = (opt_col && (pt->color))? TUNNEL : pt->celltype;
#else
        *pix++ = pt->celltype;
        //affichage du gradient de vitesse au sol
        //val = (int)((grdv[i+j*(L-2)]/coef)+20);
        //if ((i<L-2) && (j<L-2) && (val>=0) && (val<=H-1)) *pix++ = color_alti[H-1-val]; else *pix++ = 1;
#endif
      }
    }
  }
  memset(pix,BORD,img_w);
  //pix += img_w;
  cur_image += P*img_w;

  if (opt_cv){
    pix = cur_image;
    //coupe verticale, orientation est-ouest
    for(j=0; j < H; j++){
      pt = TE + HL*(prof_cv) + j*L;
      imax = Min(L,img_w);
      pix = cur_image + j*img_w;
      for(i=0; i < imax; i++, pt++){
#ifdef CELL_COLOR
        *pix++ = (pt->color) ? TUNNEL : pt->celltype;
#else
        *pix++ = pt->celltype;
#endif
      }
    }
    if (LNS>1){
    //coupe verticale, orientation nord-sud
      imax = Min(P,img_w);
      for(i=0; i < imax; i++){
        pt = TE + i*HL + abs_cv;
        pix = image + i*img_w + L;
        for(j=0; j < H; j++, pt+=L) {
        //imax = (P<=img_w) ? P : img_w;
#ifdef CELL_COLOR
          *pix++ = (pt->color) ? TUNNEL : pt->celltype;
#else
          *pix++ = pt->celltype;
#endif
        }
      }
      // empty square (antibug)
      view_rect(cur_image + L, H, H, BORD);
    }


#ifdef ROTATIONS
    if ((vdir_mode != VDIR_NONE) && lgca_ready){
      //display a representation of current orientation
      extern float csp_angle;
      int Hd2 = H/2;
      int l1 = Hd2;
      int l2 = l1/2;
      float alpha = 0.0;
      if (vdir_mode == VDIR_NORTH){
        //TODO
      } else if (vdir_mode == VDIR_WIND){
        alpha = -csp_angle*PI/180.0;
      }
      float co = cosf(alpha);
      float si = sinf(alpha);
      pix = cur_image + L + Hd2 + Hd2*img_w;
      pix -= (int)(0.5*l1*co);
      pix -= (int)(0.5*l1*si)*img_w;
      view_line(pix, l1*co, l1*si, COUPE);
      pix += (int)(0.5*l2*si);
      pix -= (int)(0.5*l2*co)*img_w;
      view_line(pix, -l2*si, l2*co, COUPE);
    }
#endif

    if (opt_lc){
      //lignes de coupe sur la vue horizontale
      //memset(image+L*(1+prof_cv-LN), COUPE, L);
      pix = image + img_w*(prof_cv) + 1;
      for (i=1; i<img_w-1; i+=ilc, pix+=ilc)
        if (i!=L-1) *pix = COUPE;

      pix = image + img_w + abs_cv;
      for (j=1; j<P-2+H; j+=ilc, pix+=ilc*img_w)
        if (j!=P-1) *pix = COUPE;

      //lignes de coupe horizontale sur une coupe verticale
      if (!opt_tr && !opt_al && !opt_lal && !opt_ls){
        //memset(cur_image+hh*L, COUPE, L);
        pix = cur_image + (H-1-hh)*img_w;
        imax = Min(L,img_w);
        for (i=0; i<imax; i+=ilc, pix+=ilc){
          *pix = COUPE;
        }
      }
    }
    cur_image += H*img_w;
  }
  else if (opt_vel && opt_lc){
    //ligne de coupe sur la vue horizontale
    pix = image + img_w*(prof_cv) + 1;
    for (i=1; i<img_w-1; i+=ilc, pix+=ilc)
      if (i!=L-1) *pix = COUPE;
  }

#ifdef LGCA
//#if 0
  unsigned char *img_col = cur_image;//image + img_w*(img_h-H);
  if (opt_vel) cur_image += H*img_w;
  if (horizontal_display)
    img_col = image + L + P*img_w + (opt_cv? H : 0);
  else if (opt_vss)
    view_rect(img_col + L, img_w - L, H, BORD);

#if (VSTEP_H > 1)
#if 0
  //affichage du controle lie au transport
  if (opt_vel){
  //if (0){
    unsigned char *pix0;
    int index, vcol;
    pix0 = img_col;
    index = HL*prof_cv;//(CLNS-1);//(CLNS>>1);
    for(i=0; (i<HL) ; i++, index++, pix0++){
      //*pix0 = check_mvt(index, NULL) ? COUPE : check_mvt_solid(index, NULL) ? DUM : BORD;
      //*pix0 = check_mvt_ouest(index, NULL) ? GRJ : check_mvt_est(index, NULL) ? EAUT : check_mvt_solid(index, NULL) ? DUM : BORD;
      //*pix0 = check_mvt_haut(index, NULL) ? GRJ : check_mvt_bas(index, NULL) ? EAUT : check_mvt_solid(index, NULL) ? DUM : BORD;
      *pix0 = check_mvt_haut(index, NULL) || check_mvt_bas(index, NULL) ? GRJ : check_mvt_ouest(index, NULL) || check_mvt_est(index, NULL) ? EAUT : check_mvt_solid(index, NULL) ? DUM : BORD;
    }
  }
#elif 1
  //affichage de la vitesse moyenne interpolee en chaque point
  if (opt_vel && Velx_interp && Vely_interp && color_alti){
    if ((prof_cv>=LN) && (prof_cv<=LS)){
      unsigned char *pix0;
      int index, vcol;
      int coef = VSTEP_H*VSTEP_L*VSTEP_TIME/50;
      pix0 = img_col;
      index = HL*prof_cv;//(CLNS-1);//(CLNS>>1);
      for(j=0; j<H; j++){
        pix0 = img_col + j*img_w;
        for(i=0; (i<L) ; i++, index++, pix0++){
          //vcol = abs(Velx_interp[index]/coef);
          vcol = (int)(sqrtf(Velx_interp[index]*Velx_interp[index]+Vely_interp[index]*Vely_interp[index])/coef);
          //vcol = (H/2)+4*((int)(sqrtf(Velx_interp[index]*Velx_interp[index]+Vely_interp[index]*Vely_interp[index])/coef)-75); //view_pression
          if (vcol >= H) vcol = H-1;
          //if (vcol < COUPE) LogPrintf("vcol=%d\n", vcol);
          //*pix0 = check_mvt_solid(index, NULL) ? DUM : color_alti[H-1-vcol];
          *pix0 = check_mvt_solid(index, NULL) ? BORD : (Velx_interp[index]<0) ? EAUT : color_alti[H-1-vcol];
          //*pix0 = check_mvt_solid(index, NULL) ? TE[index].celltype : (Velx_interp[index]<0) ? EAUT : color_alti[H-1-vcol];
          //*pix0 = check_mvt_solid(index, NULL) ? DUM : (vcol<0) ? EAUT : color_alti[H-1-vcol];
          //*pix0 = check_mvt_solid(index, NULL) ? DUM : (vcol<75) ? EAUT : (vcol>80) ? GRJ : color_alti[H-1-vcol];
        }
      }
    }
    else{
      for(j=0; j<H; j++) memset(img_col + j*img_w, DUM, L);
    }
  }

  if ((opt_vel >= 2) && Velx && Vely){
    //affichage du champ des vecteurs vitesses sur une coupe verticale est-ouest
    int iv, jv;
    int *pvx0, *pvy0, *pvx, *pvy;
    unsigned char *pix0;
    int vort;
    //int fv=VSTEP*VSTEP_TIME*NB_MVT_EO*NB_MVT_VER/4;
    int fv=VSTEP_H*VSTEP_L*VSTEP_TIME/20;
    int planz = (int)(prof_cv-LN+1)/DIST_MVT_NS; //(CLNS>>1); //LNS-1;//(CLNS>>1);
    Check_min_max(planz, 0, CLNS-1);
    //int fy=2*VSTEP_TIME*NB_MVT_EO*NB_MVT_VER;
    //memset(img_col, BORD, HL);
    pvx0 = Velx+planz*VHL;
    pvy0 = Vely+planz*VHL;
    pix0 = img_col; // + (VSTEP>>1)*img_w + (VSTEP>>1);
    //view_rect(img_col, L, H, BORD);
    for(jv=1; jv<VH-1; jv++){ //hauteur
      pix = pix0 + jv*VSTEP_L*img_w + VSTEP_H;
      pvx = pvx0 + jv*VL + 1;
      pvy = pvy0 + jv*VL + 1;
      for(iv=1; iv<VL-3; iv+=2, pvx+=2, pvy+=2, pix += 2*(VSTEP_L/NB_MVT_EO)){ //largeur
        if ((jv<VH-1) && (iv<VL-1)){
          //view_line(pix, (*pvx)>>1, (*pvy)>>1, 0);
          //view_line(pix, ((*pvx)/VSTEP_TIME)>>1, ((*pvy)/VSTEP_TIME)>>1, 0);
          view_line(pix, (*pvx)/fv, (*pvy)/fv, TUNNEL);
	  /*if ((iv>0) && (jv>0)){
	    //calcul vorticite
          vort = (*pvy) - (*(pvy-1)) - ((*pvx) - (*(pvx-VL)));
          if (vort > 10) *pix = *(pix+1) = *(pix+L) = *(pix+L+1) = MOINS;
          else if (vort < -10) *pix = *(pix+1) = *(pix+L) = *(pix+L+1) = DUM;
          }*/
        }
      }
    }
  }
#else
  //affichage des particules en mouvement
    unsigned char *pix0;
    int ii0, ii, x, y;
    pix0 = img_col;
    ii0 = (HL*NB_MVT_EO)*(int)((prof_cv-LN)/DIST_MVT_NS);//+(L*NB_MVT_EO*0.15);
    //ii0 = (HL*NB_MVT_EO)*(L>>1)+(L*(NB_MVT_EO-1));
    for(y=0; y<H ; y++, ii0+=L*NB_MVT_EO){
      for(x=0, ii=ii0; x<L; x++, pix0++, ii+=NB_MVT_EO){
        *pix0 = mvt_ouest(ii) ? EAUT : mvt_est(ii) ? GRJ : mvt(ii) ? EAUC : mvt_solid(ii) ? DUM : GR;
        //*pix0 = mvt_haut(ii) ? GRJ : mvt_est(ii) ? EAUT : mvt_solid(ii) ? GR : EAUC;
        //*pix0 = mvt_haut(ii) ? GRJ : mvt_bas(ii) ? EAUT : mvt_solid(ii) ? DUM : BORD;
        //*pix0 = mvt_haut(ii) || mvt_bas(ii) ? GRJ : mvt_ouest(ii) || mvt_est(ii) ? EAUT : mvt_solid(ii) ? DUM : BORD;
      }
    }
#endif

  //affichage des normales
  if (opt_vel && opt_vss && alti && norm3d && (prof_cv>=LN) && (prof_cv<=LS) && !reorient_flag){
  //if (0){
    float ln=10;
#ifdef CGV
    extern float dist_grdv;
    ln = dist_grdv;
#endif
    unsigned char *pix0;
#ifdef NORM2D
    Vec2 *pn;
#else
    Vec3 *pn;
#endif
    short *alt, x, y, z, stepx;
    z = prof_cv-LN;
    stepx = 10;//10; //VSTEP*2; //*2;
    alt = alti+z*LEO+(stepx>>1);
#ifdef NORM2D
    pn = &norm2d[z*LEO+(stepx>>1)];
#else
    pn = &norm3d[z*LEO+(stepx>>1)];
#endif
    for(x=(stepx>>1); x<LEO; x+=stepx, alt+=stepx, pn+=stepx){
    //for(x=1; x<L-2; x+=2, alt+=2, pnx+=2, pny+=2){
      y = (H-1)-(*alt);
      pix = img_col + y*img_w + x + 1;
      view_line(pix, (int)(ln*(pn->x))/*/NB_MVT_EO*/, (int)(-ln*(pn->y)), GRJ);
      //if (*pnx > 0.5) LogPrintf("nx = %f   ny = %f\n", *pnx, *pny);
    }
  }
//#endif
//#endif
#else // !(VSTEP_H > 1)
  //affichage de la vitesse moyenne par cellule
  if (opt_vel && Velx_sum && Vely_sum){
    int *pvx, *pvy;
    unsigned char *pix0;
    int index, vv;
    pvx = Velx_sum;
    pvy = Vely_sum;
    pix0 = img_col;
    index = HL*prof_cv;
    for(i=0; i<HL; i++, index++, pvx++, pvy++){
      *pix0++ = check_mvt_mean(index, NULL) ? EAUT : check_mvt_solid(index, NULL) ? DUM : 0;
    }
  }
#endif //LGCA


#ifdef CGV
  //affichage du gradient de vitesse en surface
  extern float *grdv;
  extern float grdvc;

  if (horizontal_display){
    // horizontal disposition
    cur_image = image + L + LN*img_w;
    if (opt_cv) cur_image += H;
  }

  if (opt_vss){
    if (grdv && lgca_ready){
      memset(cur_image, BORD, L);
      //pix = cur_image+img_w;
      float *ptgv;
      int vcol, jj, jc;
      for(j=0; j < LNS; j++){
        if (!psol[LN+j]){ //test si plan solide
          ptgv = grdv + j*LEO;
          pix = cur_image+(j+1)*img_w;
          *pix++ = BORD;
          for(i=0; i < LEO; i++, ptgv++){
            vcol = (*ptgv)*(H/4)/grdvc;
            *pix++ = (vcol < 0) ? EAUT : (vcol < H-1) ? color_alti[H-1-vcol] : GRJ;
          }
          *pix++ = BORD;
        }
        //else pix += L;
      }
      memset(cur_image + LNS*img_w, BORD, L);
      if (opt_lc && (prof_cv>=LN) && (prof_cv<=LS)){
        //ligne de coupe
        //extern int CHL;
        //int cx, cy, cz;
        //Calcule_cxyz(HL*prof_cv+1, cx, cy, cz);
        //int cz = (int)(prof_cv-LN)/DIST_MVT_NS;
        //memset(cur_image+img_w+img_w*cz*DIST_MVT_NS, COUPE, L);
        //memset(cur_image+img_w*(1+prof_cv-LN), COUPE, L);
        pix = cur_image+img_w*(1+prof_cv-LN);
        for(i=0; i<L; i+=ilc, pix +=ilc) *pix = COUPE;
      }
    }
    else{
      view_rect(cur_image, LEO, LNS, BORD);
    }
  }
#ifndef NORM2D
  extern float dist_grdv;
  //affichage des normales vues d'en haut
  if (opt_vss && alti && norm3d && !reorient_flag){
  //if (0){
    float ln=10;
    ln = dist_grdv;
    unsigned char *pix0;
    Vec3 *pn;
    short x, z, stepx, stepz;
    stepx = 10;
    stepz = 10;
    for(z=(stepz>>1); z<LNS; z+=stepz){
      pn = &norm3d[z*LEO+(stepx>>1)];
      for(x=(stepx>>1); x<LEO; x+=stepx, pn+=stepx){
        pix = cur_image + (z+1)*img_w + x + 1;
        view_line(pix, (int)(ln*(pn->x))/*/NB_MVT_EO*/, (int)(ln*(pn->z)), GRJ);
      }
    }
  }
#endif //NORM2D
#endif //CGV

#endif //LGCA

  return image;
}

// Update the position of vertical layers
int update_cv(int x, int y)
{
  int ret=0;
  //if ((x>=0) && (x<L) && (y>=0) && (y < P+2*H)){
  if ((y < P) || ((x < L) && (y < P+H))){
    if ((x < L) && (y < P)){
      abs_cv = x;
      prof_cv = y;
    }
    else if (y > P){
      abs_cv = x;
      if (opt_h) hh = P+H-y;
    }
    else{
      prof_cv = y;
    }
    if (abs_cv < 1) abs_cv = 1;
    if (abs_cv > L-2) abs_cv = L-2;
    if (prof_cv < 1) prof_cv = 1;
    if (prof_cv > P-2) prof_cv = P-2;
    ret = 1;
  }
  return ret;
}

#ifdef USE_LIBPNG

#include <png.h>

png_color* png_palette;

void view_dump_init()
{
  int colors[256];
  int i;

  LogPrintf("view_dump_init\n");

  AllocMemory(png_palette, png_color, 256);
  ResetMemory(png_palette, png_color, 256);

  view_palette(colors);

  for(i=0; i<256; i++){
    png_palette[i].red   = (colors[i] & 0x00ff0000) >> 16;
    png_palette[i].green = (colors[i] & 0x0000ff00) >> 8;
    png_palette[i].blue  = (colors[i] & 0x000000ff);
  }
}

void view_dump_png(char *filename)
{
  unsigned char *pix;
  int j;

  LogPrintf("view_dump_png: writing file %s\n", filename);

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

  png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
  if (!png_ptr){
    ErrPrintf("ERROR: cannot create structure png_struct\n");
    exit(-4);
  }

  png_infop info_ptr = png_create_info_struct(png_ptr);
  if (!info_ptr)
  {
    png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
    ErrPrintf("ERROR: cannot create structure png_info\n");
    exit(-4);
  }

  if (setjmp(png_jmpbuf(png_ptr)))
  {
    png_destroy_write_struct(&png_ptr, &info_ptr);
    fclose(fp);
    ErrPrintf("ERROR: cannot set jump in view_dump_png()\n");
    exit(-4);
  }

  png_init_io(png_ptr, fp);

  png_set_IHDR(png_ptr, info_ptr, img_w, img_h,
   8, PNG_COLOR_TYPE_PALETTE, PNG_INTERLACE_NONE,
   PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);

  png_set_PLTE(png_ptr, info_ptr, png_palette, 256);

  png_write_info(png_ptr, info_ptr);
  //png_write_png(png_ptr, info_ptr, 0, NULL);

  pix = view();

  for(j=0; j<img_h; j++, pix+=img_w)
    png_write_row(png_ptr, (png_byte*)pix);

  //png_write_flush(png_ptr);

  png_write_end(png_ptr, NULL);

  png_destroy_write_struct(&png_ptr, &info_ptr);

  fclose(fp);
}

void dump_image(char *filename, char *format)
{
  if (!strcmp(format,"png"))
    view_dump_png(filename);
  else{
    ErrPrintf("ERROR: jpeg format not supported\n");
    exit(-1);
  }
}

#elif defined(USE_GD)

#include <gd.h>

gdImagePtr gd_img=NULL;
int gd_col[256];

void view_dump_init()
{
  int colors[256];
  int i, red, green, blue;

  gd_img = gdImageCreate(img_w, img_h);

  view_palette(colors);

  for(i=0; i<256; i++){
    red   = (colors[i] & 0x00ff0000) >> 16;
    green = (colors[i] & 0x0000ff00) >> 8;
    blue  = (colors[i] & 0x000000ff);
    gd_col[i] = gdImageColorAllocate(gd_img, red, green, blue);
  }
}

void view_gd()
{
  unsigned char *pix;
  int i, j;

  assert(gd_img);

  pix = view();

  for(j=0; j<img_h; j++)
    for(i=0; i<img_w; i++)
      gdImageSetPixel(gd_img, i, j, gd_col[*pix++]);
}

void view_dump_png(char *nom)
{
  FILE *fp;

  view_gd();
  fp = fopen(nom,"w");
  if ( ! fp ){
    ErrPrintf("erreur ouverture fichier png\n");
    exit(-4);
  }
  gdImagePng(gd_img, fp);
  fclose(fp);
}

void view_dump_jpeg(char *nom)
{
  FILE *fp;

  view_gd();
  fp = fopen(nom,"w");
  if ( ! fp ){
    ErrPrintf("erreur ouverture fichier jpg\n");
    exit(-4);
  }
  gdImageJpeg(gd_img, fp, 50); /* 50 -> medium quality and compression */
  fclose(fp);
}

void dump_image(char *filename, char *format)
{
  if (!strcmp(format,"png"))
    view_dump_png(filename);
  else
    view_dump_jpeg(filename);
}

#else

#include <gtk/gtk.h>

extern gint area_w, area_h;
extern GtkWidget *drawingarea;

void dump_image(char *filename, char *format)
{
  //LogPrintf("dump_image: %s - %s\n",filename,format);
  GdkPixbuf *pixbuf = gdk_pixbuf_get_from_drawable(NULL,drawingarea->window,NULL,0,0,0,0,area_w,area_h);
  gdk_pixbuf_save(pixbuf, filename, format, NULL, NULL);
}

#endif //USE_GD

void dump_image_inter(int inter, char *format)
{
  static int cpt_inter=0;
  static int cpt_snap=0;
  static char nom[100];
  static char str[100];
  static char flag_sec=0;
  int nmin;
  //char *ext;

#ifdef PARALLEL
  extern int proc_id;
  sprintf(str,"_%d",proc_id);
#else
  *str = 0;
#endif

  //ext = (format==IMG_PNG)? "png" : "jpg";
  if (inter){
    if (!cpt_inter) flag_sec = inter % 60; //fmodf(inter,60);
    nmin = cpt_inter/60;
    if (!flag_sec)
      sprintf(nom,"%s%04d%s.%s", MOD_NAME, nmin, str, format);
    else
      sprintf(nom,"%s%04d-%02d%s.%s", MOD_NAME, nmin, cpt_inter - 60*nmin, str, format);
    cpt_inter += inter;
  }
  else{
    sprintf(nom,"SNAP%03d.%s", cpt_snap++, format);
  }

  dump_image(nom, format);
}

void view_quit()
{
  if (opt_al || opt_lal) free(color_alti);
  if (image) free(image);
#ifdef USE_GD
  if (gd_img) free(gd_img);
#endif
}

