#include<stdio.h>
#include <stdlib.h> /* pour malloc() */
#include<math.h>
#include <sys/stat.h> /* */
#include <string.h> /* pour memcpy() */

#include"lire_fort_util.h"

/* compilation:
gcc fort2tec.c -o fort2tec -O3 -Wall
*/

/*--------------------------------------------------------------------*/

void aide(char *progname)
{
  fprintf(stderr,"Utilisation de %s :\n",progname);
  fprintf(stderr,"%s 'fichierfort' 'fichiertec' \n",progname);
  fprintf(stderr,"\t ou 'fichierfort' est le nom du fichier a lire\n");
  fprintf(stderr,"\t et 'fichiertec' est le nom du fichier Tecplot produit\n");
  fprintf(stderr,"\t (sans 'fichiertec', le fichier de sortie sera 'champs.tec')\n");
}

/*--------------------------------------------------------------------*/

void check_cmdline(int argc,char **argv){
/* verification de l'exactitude de l'appel du programme */

if(argc==1){ /* il manque le nom du fichier de donnees a lire */
  fprintf(stderr,"!!! Erreur (fatale)!!!\n");
  aide(argv[0]);
  exit(0); /* arret du programme*/
}

if(argc>=4){ /* il y a trop d'arguments */
  fprintf(stderr,"!!! Erreur (trop d arguments) !!!\n");
  aide(argv[0]);
  exit(0); /* arret du programme*/
}

}

/* ------------------------------------------------------------------- */

void build_pgl(double *pgl,int taille){
/********************************************************/
/* construction des points de Gauss-Lobatto sur [-1;+1] */
/********************************************************/
/* Variables en entree: 				*/
/*	taille	: nombre de points 			*/
/********************************************************/
/* Variables en sortie:					*/
/*	pgl[]	: abscisses (sur [-1;+1]) des points	*/
/********************************************************/
char func_name[]="build_pgl";
int i;
int demitaille;

if(taille<2){
  fprintf(stdout,"%s: taille:%d !!!\n",func_name,taille);
  exit(0);
}
/* points de Gauss Lobatto (sur [-1;+1]) */
/* construction brute:
for(i=0;i<taille;i++){
   pgl[i]=-cos(i*M_PI/(taille-1));
} */
pgl[0]=-1;
pgl[(taille-1)]=1;
if((taille%2)==0) { /* taille est pair */
  demitaille=taille/2;
  for(i=1;i<demitaille;i++){
    pgl[i]=-cos(i*M_PI/(taille-1));
    pgl[(taille-1)-i]=-pgl[i];
  }
}
else { /* taille est impair */
  demitaille=(taille-1)/2;
  for(i=1;i<demitaille;i++){
    pgl[i]=-cos(i*M_PI/(taille-1));
    pgl[(taille-1)-i]=-pgl[i];
  }
  pgl[demitaille]=0;
}


}

/*--------------------------------------------------------------------*/

int main(int argc,char **argv){
int Nx; /* degre polynomial suivant x */
int Ny; /* degre polynomial suivant y */
int Nz; /* degre polynomial suivant z */
int taille_x; /* taille_x=Nx+1 */
int taille_y; /* taille_y=Ny+1 */
int taille_z; /* taille_z=Nz+1 */
double Ax,Ay,Az;
double *pgl_x,*pgl_y,*pgl_z; /* points de Gauss lobatto dans [-1,1] */
double *pos_x,*pos_y,*pos_z; /* points dans [-Ax;Ax] et [-Ay;Ay] [-Az;Az] */

/* variables "fortran" lues (c.f. writedisk.f) */
int NDIM,NVAR;
double PR,RA;
double DT;

double **champ2D_U,**champ2D_W,**champ2D_T;
double ***champ3D_U,***champ3D_V,***champ3D_W,***champ3D_T;
char *infile;
char *outfile;
char *outfile_default="champs.tec";
FILE *fout;

int i,j,k;

/* 1. Initialisations */
/* 1.1 Verification de la ligne de commande*/
check_cmdline(argc,argv);

/* 1.2 Lecture des donnees (1ere partie, en-tete du fichier fort) */
infile=argv[1];
fprintf(stdout,"Fichier d'entrees %s\n",infile);
if (argc == 2) { /* user did not specify any output file name*/
  outfile=outfile_default;
}
else {
  outfile=argv[2];
}
fprintf(stdout,"Fichier de sorties %s\n",outfile);

lire_fort_header(infile,&NDIM,&NVAR,&Nx,&Ny,&Nz,
	&Ax,&Ay,&Az,&PR,&RA,&DT);

/* 1.3 Calcul de la taille des systemes */
taille_x=Nx+1;
taille_y=Ny+1;
if (NDIM == 3) {
  taille_z=Nz+1;
}

/* 1.4 Allocations dynamiques */
/* vecteurs */
if ((pgl_x=(double*)malloc(sizeof(double)*taille_x))==NULL)
   fprintf(stderr,"memory allocation failed for pgl_x[]\n");
if ((pos_x=(double*)malloc(sizeof(double)*taille_x))==NULL)
   fprintf(stderr,"memory allocation failed for pos_x[]\n");
if ((pgl_y=(double*)malloc(sizeof(double)*taille_y))==NULL)
   fprintf(stderr,"memory allocation failed for pgl_y[]\n");
if ((pos_y=(double*)malloc(sizeof(double)*taille_y))==NULL)
   fprintf(stderr,"memory allocation failed for pos_y[]\n");
if (NDIM == 3) {
  if ((pgl_z=(double*)malloc(sizeof(double)*taille_z))==NULL)
     fprintf(stderr,"memory allocation failed for pgl_z[]\n");
  if ((pos_z=(double*)malloc(sizeof(double)*taille_z))==NULL)
     fprintf(stderr,"memory allocation failed for pos_z[]\n");
}

/* champs 2D */
if (NDIM == 2){
 if ((champ2D_U=(double**)malloc(sizeof(double*)*taille_x))==NULL)
   fprintf(stderr,"memory allocation failed for champ2D_U[]\n");
 else
   for(i=0;i<taille_x;i++){
   if ((champ2D_U[i]=(double*)malloc(sizeof(double)*taille_y))==NULL)
   fprintf(stderr,"memory allocation failed for champ2D_U[%d][]\n",i);
   }
 if ((champ2D_W=(double**)malloc(sizeof(double*)*taille_x))==NULL)
   fprintf(stderr,"memory allocation failed for champ2D_W[]\n");
 else
   for(i=0;i<taille_x;i++){
   if ((champ2D_W[i]=(double*)malloc(sizeof(double)*taille_y))==NULL)
   fprintf(stderr,"memory allocation failed for champ2D_W[%d][]\n",i);
   }
 if ((champ2D_T=(double**)malloc(sizeof(double*)*taille_x))==NULL)
   fprintf(stderr,"memory allocation failed for champ2D_T[]\n");
 else
   for(i=0;i<taille_x;i++){
   if ((champ2D_T[i]=(double*)malloc(sizeof(double)*taille_y))==NULL)
   fprintf(stderr,"memory allocation failed for champ2D_T[%d][]\n",i);
   }
}

/* champs 3D */
if (NDIM == 3) {
 if ((champ3D_U=(double***)malloc(sizeof(double**)*(Nx+1)))==NULL)
   fprintf(stderr,"memory allocation failed for champ3D_U[]\n");
 else
   for(i=0;i<(Nx+1);i++){
   if ((champ3D_U[i]=(double**)malloc(sizeof(double*)*(Ny+1)))==NULL)
   fprintf(stderr,"memory allocation failed for champ3D_U[%d][]\n",i);
   else
    for(j=0;j<(Ny+1);j++){
    if ((champ3D_U[i][j]=(double*)malloc(sizeof(double)*(Nz+1)))==NULL)
    fprintf(stderr,"memory allocation failed for champ3D_U[%d][%d][]\n",i,j);
    }
   }
 if ((champ3D_V=(double***)malloc(sizeof(double**)*(Nx+1)))==NULL)
   fprintf(stderr,"memory allocation failed for champ3D_V[]\n");
 else
   for(i=0;i<(Nx+1);i++){
   if ((champ3D_V[i]=(double**)malloc(sizeof(double*)*(Ny+1)))==NULL)
   fprintf(stderr,"memory allocation failed for champ3D_V[%d][]\n",i);
   else
    for(j=0;j<(Ny+1);j++){
    if ((champ3D_V[i][j]=(double*)malloc(sizeof(double)*(Nz+1)))==NULL)
    fprintf(stderr,"memory allocation failed for champ3D_V[%d][%d][]\n",i,j);
    }
   }
 if ((champ3D_W=(double***)malloc(sizeof(double**)*(Nx+1)))==NULL)
   fprintf(stderr,"memory allocation failed for champ3D_W[]\n");
 else
   for(i=0;i<(Nx+1);i++){
   if ((champ3D_W[i]=(double**)malloc(sizeof(double*)*(Ny+1)))==NULL)
   fprintf(stderr,"memory allocation failed for champ3D_W[%d][]\n",i);
   else
    for(j=0;j<(Ny+1);j++){
    if ((champ3D_W[i][j]=(double*)malloc(sizeof(double)*(Nz+1)))==NULL)
    fprintf(stderr,"memory allocation failed for champ3D_W[%d][%d][]\n",i,j);
    }
   }
 if ((champ3D_T=(double***)malloc(sizeof(double**)*(Nx+1)))==NULL)
   fprintf(stderr,"memory allocation failed for champ3D_T[]\n");
 else
   for(i=0;i<(Nx+1);i++){
   if ((champ3D_T[i]=(double**)malloc(sizeof(double*)*(Ny+1)))==NULL)
   fprintf(stderr,"memory allocation failed for champ3D_T[%d][]\n",i);
   else
    for(j=0;j<(Ny+1);j++){
    if ((champ3D_T[i][j]=(double*)malloc(sizeof(double)*(Nz+1)))==NULL)
    fprintf(stderr,"memory allocation failed for champ3D_T[%d][%d][]\n",i,j);
    }
   }
}

/* 1.5 Construction des coordonnees */
build_pgl(pgl_x,taille_x);
for(i=0;i<taille_x;i++) {
  pos_x[i]=Ax*pgl_x[i];
}
build_pgl(pgl_y,taille_y);
for(i=0;i<taille_y;i++) {
  pos_y[i]=Ay*pgl_y[i];
}
if (NDIM == 3) {
  build_pgl(pgl_z,taille_z);
  for(i=0;i<taille_z;i++){
    pos_z[i]=Az*pgl_z[i];
  }
}

/* 2. Lecture des donnees (2eme partie, les champs) */
if (NDIM == 2) {
  lire_fort_data2D(infile,Nx,Ny,champ2D_U,champ2D_W,champ2D_T);
}
else {
  lire_fort_data3D(infile,Nx,Ny,Nz,champ3D_U,champ3D_V,
                   champ3D_W,champ3D_T);
}

/* 3. Ecriture du fichier tecplot */
fout=fopen(outfile,"w");

if (NDIM == 2) {
  fprintf(fout,"TITLE=\"CHAMPS DANS LE DOMAINE 2D\"\n");
  fprintf(fout,"VARIABLES=\"X\",\"Y\",\"U\",\"W\",\"T\"\n");
  fprintf(fout,"ZONE T=\"2D\",I=%d J=%d\n",taille_x,taille_y);
  /* ecriture des coordonnees et champs */
  for(j=0;j<taille_y;j++){
    for(i=0;i<taille_x;i++){
    fprintf(fout,"%15.8E %15.8E %15.8E %15.8E %15.8E\n",
      pos_x[i],pos_y[j],
      champ2D_U[i][j],champ2D_W[i][j],champ2D_T[i][j]);
    }
  }
}
else {/* 3D datasets */
  fprintf(fout,"TITLE=\"CHAMPS DANS LE DOMAINE 3D\"\n");
  fprintf(fout,"VARIABLES=\"X\",\"Y\",\"Z\",\"U\",\"V\",\"W\",\"T\"\n");
  fprintf(fout,"ZONE T=\"3D\",I=%d J=%d K=%d\n",
               taille_x,taille_y,taille_z);
  /* ecriture des coordonnees et champs */
  for(k=0;k<taille_z;k++){
    for(j=0;j<taille_y;j++){
      for(i=0;i<taille_x;i++){
      fprintf(fout,"%15.8E %15.8E %15.8E %15.8E %15.8E %15.8E %15.8E\n",
        pos_x[i],pos_y[j],pos_z[k],
        champ3D_U[i][j][k],champ3D_V[i][j][k],champ3D_W[i][j][k],
        champ3D_T[i][j][k]);
      }
    }
  }
}

fclose(fout);

return 0;
}