#include <stdio.h>
#include <stdlib.h>
#include "math.h"
#include <math.h> 
// cc onde1.c ; ./a.out | gnuplot

// OK MARCHE
//
// attention les differences finies c est pas bien
//
//
// solution de  du/dt = - deta/dx;  deta/dt = -du/dx 


// cas instable; jouer avec dt pour stabiliser

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





    
int main (int argc, const char *argv[]) {
	
    int  i,nx,it;
    char file;    
    double*x=NULL,*u=NULL,*uo=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.4; 
     t = 0;    
     tmax=3;	
     
        
  
      
      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));



     fprintf(stderr," Coucou solution de  du/dt = - deta/dx;  deta/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];
     } 
    

 
	
	do{  
		t=t+dt;	
		it++;
		
		// resolution simple                  
		for(i=1;i<nx-1;i++)
             {  
               eta[i] = etao[i]   - dt*(uo[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 ;               
              } 
        u[0]=0;
        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(" p[%lf:%lf][%lf:%lf]  '-' u 1:2 t'u'  w l linec 3 , '-' u 1:2 t'eta' w lp\n ",
			   x[1],x[nx],-.01,.01); 
		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");

    }
	
	}while(t<tmax);
	
	 
}