#define CRLF "\n"    
/*  novembre 2004        */
/* 26 octobre 2002       */
/*                       */ 
/* jv 2002               */
/* Sept 2001             */
/* 4 3 1                  */
/* fev 2001 
   06/11/2000                                                        */
/* 01 /2001                                                          */
/* Nouvelle version 12/05/2000                                       */
/* OK TF issue de KS                                                 */
/* 26/09/95   vers 19 Sept 95                                        */
/* 26/09/95     OK 26/10/96                                          */
/*                                                                   */
/*                                                                   */
/* LIRE LE LIVRE NUMERICAL RECIPIES                                  */
/* sur www.lmm.jussieu.fr/INFO2/info.html: source de Num Rec.        */
/*   
/*  COMPILATION dans le fichier Fait.txt:
cc -O3 -ffast-math dfour1.c  fftTD_b.c -lm -o ff.out;ff.out
                                                                     */
/* utilisation de Numˇrical Recipies (en Fortran!) p 390!            */                                                                 

/* dfour1.c                                                          */
/* transformee de Fourier en DOUBLE precision                        */
/*            j=N-1                                                  */
/*             ___                                                   */
/*      ^      \          I k (2 Pi j/N)                             */
/*      u   =  /    u    e                                           */
/*       k     ---   j                                               */
/*             j=0                                                   */
/*                                                                   */
/* remarquer que les valeurs sont prises en xj=(2 Pi j/N)            */ 
/*                                                                   */
/* La TF inverse s'ecrit alors                                       */
/*            k=N-1                                                  */
/*             ___                                                   */
/*             \    ^     -I k (2 Pi j/N)                            */
/*      u   =  /    u    e                                           */
/*       j     ---   k                                               */
/*             k=0                                                   */
/*                                                                   */
/* La derivation s'ecrit alors                                       */
/*            k=N-1                                                  */
/*             ___                                                   */
/*     d       \    ^            -I k (2 Pi j/N)                     */
/*     -- u =  /    u  (- I k)  e                                    */
/*     dx  j   ---   k                                               */
/*             k=0                                                   */
/*                                                                   */
/*                                                                   */
/*         ON MULTIPLIE PAR (-Ik)  POUR DERIVER                      */
/*                                                                   */
/*      ACHTUNG !!!!!                                                */
/*            les variations sont entre -Pi et Pi                    */
/*                     => correction des dˇrivˇes                    */
/*                    d'ou le facteur    (Pi/Lx)                     */
/* ================================================================  */ 
 

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

/*---------------------------------------------------------*/
double* u=NULL,* x=NULL,* yw=NULL,* ywp=NULL;
double* TFf=NULL,* TFfp=NULL,* TFp=NULL,*TFz=NULL,*Esp=NULL;
double* p=NULL;
double* tau=NULL,* TFtau=NULL;
double* Tzr=NULL,* Tzi=NULL;
/*---------------------------------------------------------*/
double zr,zi,a,a13;
double Pi;
double Lx,Lb,h0,dx;
int it,nt,reg,type;
int nx ;
#define Power pow
int npix,npiy;
FILE  *g;

  
/*---------------------------------------------------------*/
 
void trace(void);
void reLit(void);
/*---------------------------------------------------------*/

                                                        
int main(void)
{
int i,itt,nbt;
double fx,fy,qq;
 
/* ouverture des fichiers de sauvegarde  */
 
  
 
printf("     ___                                    \n");
printf("^    \\          I k (2 Pi j/N)             \n");
printf("u  = /    u    e                            \n");
printf(" k   ---    j                               \n");
printf("                                            \n");
printf("                                            \n");
printf("                                            \n");
printf("                                            \n");
printf(" tau = F[y ]                                \n");
printf("          w                                 \n");     	   
 /* lecture du fichier donnees */

/* a modifier ici, lire le nombre de points de la forme 2^n  */ 
Pi=4*atan(1);
 

/* domaine en x et en t                         */ 
reLit();
  /* RAZ   */
 
 
 
printf("  nbre de points en 2^n = 2^%d    \n",nx);
printf(" Lx = %7f (1/2 taille domaine)    \n",Lx);
printf(" Lb = %7f (1/2 longueur bosse)    \n",Lb);
printf(" dx = %7f (calcule 2*Lx/(nx-1))   \n",dx); 
printf(" regime = %d ->  ",reg);
 if (reg==-1){printf("subcritique   F<1 \n");}
 if (reg== 0){printf("tuyau             \n");}
 if (reg== 1){printf("supercritique F>1 \n");}
 if (reg== 2){printf("prof. infinie     \n");}
 if (reg== 3){printf("tuayu p2-p1=A''   \n");}
 if (reg==10){printf("supersonique M>1  \n");}
 if (reg==100){printf("Special Hilbert  \n");}	   
 printf("  type de bosse initiale  ->     ");
 if((type)==0){printf("bruit blanc\n");}
 if((type)==1){printf("gaussienne de largeur Lb\n");}
 if((type)==2){printf("sinus\n");}
 if((type)==3){printf("marche descendante\n");}
 if((type)==4){printf("diamant\n");}
 if((type)==99){printf("reprise\n");}
 



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

/* fin des lectures */



/* declaration */
nbt=nx+15;
u   =(double*)calloc(nbt,sizeof(double)); 
Tzr =(double*)calloc(nbt,sizeof(double)); 
Tzi =(double*)calloc(nbt,sizeof(double)); 
/* u est la fonction dont on cherche la TF */
nbt=2*nx+30;
x     =(double*)calloc(nbt,sizeof(double));
yw    =(double*)calloc(nbt,sizeof(double));
ywp   =(double*)calloc(nbt,sizeof(double));
tau   =(double*)calloc(nbt,sizeof(double));
p     =(double*)calloc(nbt,sizeof(double));
TFf   =(double*)calloc(nbt,sizeof(double));
TFfp  =(double*)calloc(nbt,sizeof(double));
TFp   =(double*)calloc(nbt,sizeof(double));
TFtau =(double*)calloc(nbt,sizeof(double));

/* fin declaration */

 
/* Initialisation conditions initiales        */
 
 
 for (i=0; i<=nx;i++)
 {x[i]=-Lx+ i*dx;}

 
/*   expon               */     
if((type)==1){
  for (i=0; i<=nx;i++)
 {
  yw[i]=  h0*exp(-Pi*pow((x[i])/Lb,2));
       }
 }
/*   cos               */   
if((type)==2){ 
 for (i=0; i<=nx;i++)
  {
  yw[i]= h0*cos( 4*Pi*x[i]/Lb /2) ;
       }   
  }
  
if((type)==3){
   for (i=0; i<nx;i++)
  { yw[i]=  h0/2*(tanh(- (x[i])/Lb  )-1);
 if(x[i]>Lx/2) {yw[i]=  h0/2*(-tanh(- (x[i]-3*Lx/4)/Lb  )-1);;}
 }}
 
 
if((type)==4){
   for (i=0; i<nx;i++)
  { 
   yw[i]= 0;
  if(x[i]>0 )    {yw[i]=  x[i]/(Lb/2)*h0;}
  if(x[i]>(Lb/2)){yw[i]=(-x[i]/(Lb/2) +2)*h0;}
  if(x[i]>Lb)  {yw[i]=0;}
 }} 
 
 
/* reprise             */    
  if((type)==99){ 
  g= fopen("xypt.OUT","r");	       
   for (i=0; i<nx;i++)
  { 
  fscanf(g,"%lf %lf  \n",&fx,&fy,&qq,&qq ); 
    yw[i]=fy;printf("%lf %lf \n",x[i],yw[i]);  
  }
  fclose(g);
   }
/* fin init geometrie */         



 
  
  
  
/* le calcul commence maintenant                      */   
 
 
 
 
   
/* on transvase Re[TFf]=yw; Im[TFf]=0                 */
/* la syntaxe provient de la routine FFT              */ 
  for (i=0; i<=nx;i++)
 {TFf[2*i+1]=yw[i];
  TFf[2*i+2]=0;} 
/* transforme fourier un tableau "TFf" de yw          */
  dfour1(TFf,nx,1);   

  
  hydro(TFf,TFfp,TFp,TFtau,nx,Lx,reg);
  
/* &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*/
/* $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ */  
 

  
   
/* on revient dans l'espace physique, on divise donc par nx? */ 
  dfour1(TFp,nx,-1);
  dfour1(TFtau,nx,-1);
  dfour1(TFfp,nx,-1);
	   
  for (i=1; i<=2*nx+2;i++){TFp[i]=TFp[i]/nx;TFtau[i]=TFtau[i]/nx;TFfp[i]=TFfp[i]/nx;}
 
    
    for (i=0; i<nx;i++)
 {p[i]=TFp[2*i+1] -TFp[ 1]*0.000;
 tau[i]=TFtau[2*i+1];
 ywp[i]=TFfp[2*i+1];
 
 }   

 
/* ================================================================*/ 
/* trace de la fonction                                            */ 

 
 
 trace();
 
 //getchar();
  
 
 return 0;}


/*---------------------------------------------------------*/
 
/*---------------------------------------------------------*/
void trace()
{
 extern int nxt;
 extern double *x,*tau,*p, *yw, *ywp;
static   int i;
    FILE *fa;
  
/* sauve  */
 
 
 
 
/*==================sauvegarde profils instantannes =========================================*/
   fa= fopen("xypt.OUT","w");
   for (i=0; i<nx;i++)
  { 
  fprintf(fa,"%7f %7f %7f %7f  %s",
         x[i],yw[i],p[i],tau[i],CRLF); 
  }
  fclose(fa); 
 
  
 
  
 }
 
 
 
/*---------------------------------------------------------*/ 
void reLit()
 { 
       FILE *g;
       extern int nx,type,reg;
       extern double dxLx,Lb,h0;

	   
g = fopen("FfftTD_b.IN", "r");


fscanf(g,"nx=%i             \n",&nx);
fscanf(g,"Lx=%lf            \n",&Lx);
fscanf(g,"Lb=%lf            \n",&Lb);
fscanf(g,"h0=%lf            \n",&h0);
fscanf(g,"type=%i           \n",&type);
fscanf(g,"reg=%i            \n",&reg);
 
fclose(g);
dx=2*Lx/(nx-1);   
	   }
  
/*---------------------------------------------------------*/