#include <stdio.h>
#include <stdlib.h>
#include "math.h"
#include <math.h> 
// cc onde3.c ; ./a.out | gnuplot
//
// OK MARCHE
//
// attention les differences finies c est pas bien
//
// solution de  du/dt = - deta/dx;  deta/dt = -d(hu)/dx  

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





    
int main (int argc, const char *argv[]) {
	
    int  i,nx,it;
    char file;    
    double*x=NULL,*u=NULL,*uo=NULL,*h0=NULL,*eta=NULL,*etao=NULL;
    double dt,dx,x0,L,t,tmax; 
    FILE *g;
    
 
     nx = 800;        
     L = 40;
     x0=-10;
     dx = L/nx;
     dt = 0.05; 
     t = 0;
     tmax=45; 

   	
     
        
  
      
      x= (double*)calloc(nx+1,sizeof(double));
      u= (double*)calloc(nx+1,sizeof(double));
      eta=(double*)calloc(nx+1,sizeof(double));
      uo= (double*)calloc(nx+1,sizeof(double));
      etao=(double*)calloc(nx+1,sizeof(double));
      h0= (double*)calloc(nx+1,sizeof(double));



     fprintf(stderr," Coucou solution de  du/dt = - dh/dx;  dh/dt = -du/dx   dx=%lf  dt/dx=%lf\n",dx,dt/dx);

// init 
     for(i=0;i<=nx;i++)
     { u[i]=0;
       x[i]= x0 +i*dx;
       uo[i]=0;;
       etao[i]=.01*exp(-(x[i])*(x[i]));         
       u[i]=uo[i];
       eta[i]=etao[i];
       h0[i]=1  -  ((x[i]-5)/(x0 +nx*dx))*(x[i]>5);    
     } 
    

 
	
	do{  
		t=t+dt;	
		it++;
		
		// resolution simple                  
		for(i=1;i<nx-1;i++)
             {  
               eta[i] = etao[i]   - dt*(h0[i+1]*uo[i+1]-h0[i-1]*uo[i-1])/2/dx ;               
              } 
        eta[0]=eta[1];
        eta[nx-1]=eta[nx-2];
                
        for(i=0;i<nx;i++) etao[i]=eta[i];
           
           
        for(i=1;i<nx-1;i++)
             {  
               u[i] = uo[i]   - dt*(etao[i+1]-etao[i-1])/2/dx ;               
              } 
        //sort à gauche
        u[0]=u[1];
        // bord à droite
        u[nx-1]=0;
        
       
        for(i=0;i<nx;i++) uo[i]=u[i];  
               
                   	
 
          
    if((it%1)==0){ 
    g = fopen("sol.OUT", "w");
	for (i=0; i<=nx;i++)
	{      fprintf(g,"%lf %lf %lf \n",x[i],u[i],t);}
	fprintf(g,"\n");
	fclose(g);
	
	fprintf(stderr,"  t=%lf\n",t);
	
	
	printf("t=%lf;  f(x)=0.01*exp(-(x-t)*(x-t))/2;p[%lf:%lf][%lf:%lf] f(x) t'exacte',",t,   x[1],x[nx],-.01,.01);
	printf("   '-' u 1:2 t'u'  w l, '-' u 1:2 t'eta' w lp, '-' u 1:2 t'fond' w lp\n ");
	    
		for (i=0; i<=nx;i++)
		{
		printf("%lf %lf \n",x[i],u[i]);}
		printf("e \n");
		 for (i=0; i<=nx;i++)
		{
		printf("%lf %lf \n",x[i],eta[i]);}
		printf("e \n");
		
		 for (i=0; i<=nx;i++)
		{
		printf("%lf %lf \n",x[i],-h0[i]/100);}
		printf("e \n");

    }
	
	}while(t<tmax);
	
	 
}