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

//A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows
// Emmanuel Audusse, Franccois Bouchut, Marie-Odile Bristeau, Rupert Klein and Benoıt Perthame
// biblio : Delestre These
// version init mai 2012 version mars 2013
// PYL
    double*x=NULL,*h=NULL,*u=NULL,*Q=NULL;
    double t,dt,tmax,dx,Z; 
    int nx;
 
double FR1(double ug,double ud,double hg,double hd)
 { double c;
 //Rusanov 
   c=fmax(fabs(ug)+sqrt(hg),fabs(ud)+sqrt(hd));
 // c=1*dx/2/dt;
   return (hg*ug+hd*ud)*0.5-c*(hd-hg)*0.5;
 }
double FR2(double ug,double ud,double hg,double hd)
 { double c;
 //Rusanov 
   c=fmax(fabs(ug)+sqrt(hg),fabs(ud)+sqrt(hd));
 //  c=1*dx/2/dt;
   return (ug*ug*hg + hg*hg/2. + ud*ud*hd + hd*hd/2.)*0.5 - c*(hd*ud-hg*ug)*0.5;
 }
  
/*     -------------------------------------------------    */        
int main (int argc, const char *argv[]) {
    int  i,it=0;
    FILE *g;
    double*fp=NULL,*fd=NULL,*un=NULL,*hn=NULL;
    double  q,y1,y2;
    
// parametres --------------------------------------------------------------  
dt=0.01;
tmax=10;
dx=0.25;
nx=160;
t=0;	  
	
  fprintf(stderr,"  ---------------------                         \n");	
  fprintf(stderr,"                     <--                        \n");   
  fprintf(stderr,"                        |--->                   \n");   
  fprintf(stderr,"  ----------------------------------------------\n");   
/* ------------------------------------------------------------------------*/      
      x= (double*)calloc(nx+1,sizeof(double));
      h= (double*)calloc(nx+1,sizeof(double));
      Q= (double*)calloc(nx+1,sizeof(double));
      u= (double*)calloc(nx+1,sizeof(double));
      fp=(double*)calloc(nx+1,sizeof(double));
      fd=(double*)calloc(nx+1,sizeof(double)); 
      un=(double*)calloc(nx+1,sizeof(double));
      hn=(double*)calloc(nx+1,sizeof(double));

      
 fprintf(stderr,"tmax= %lf \n",tmax);
      
// initialisation cond init ----------------------------      
  for(i=0;i<=nx;i++)
    {  x[i]=(i-nx/2)*dx;
       
   //dambreak
	 Z=0;
     h[i]=1*( x[i]<0);
     u[i]=0;
     Q[i]=u[i]*h[i];
     }
       
// initialisation du fichier de sortie ----------------------      
    
      g = fopen("solxhQZ.OUT", "w");
      fclose(g); 
                  
    while(t<tmax){   // boucle en temps
      t=t+dt;
      it=it+1;      
 
   for(i=1;i<=nx;i++)
    {    
    	  fp[i]=FR1(u[i-1],u[i],h[i-1],h[i]);
    	  fd[i]=FR2(u[i-1],u[i],h[i-1],h[i]);  
    }
   
   for(i=1;i<nx;i++)
    {
      hn[i]=h[i]- dt*(fp[i+1]-fp[i])/dx;   //conservation de la masse
  if(h[i]>0.){                             //conservation qunatité de mouvement
      q=h[i]*u[i]-dt*(fd[i+1]-fd[i])/dx ;  
      un[i]=q/hn[i];}
      else{
      un[i]=0.;}
    }
    
    // flux nul en entree sortie 
        hn[0]=hn[1];
		un[0]=un[1];	
		hn[nx]=hn[nx-1];
		un[nx]=un[nx-1];	
   
   //swap		   
  for(i=0;i<=nx;i++)
    { 
      h[i]=hn[i];
      u[i]=un[i];
      Q[i]=hn[i]*un[i];
    }

  if(it%100==0){
  	/* Saving the fields */	
	g = fopen("solxhQt.OUT", "a");
	for (i=0; i<=nx;i++)
	{      fprintf(g,"%lf %lf %lf %lf \n",x[i],h[i],Q[i],t);}
	 fprintf(g,"\n");
	 fclose(g);
 
  
    printf("t=%lf\n",t);
    
 // enlever le commentaire si on utilise gnuplot
 //./a.out | gnuplot   
 //

/*
	printf("set xlabel \"x\" ; set title \"t=%lf \"\n",t);
     y1=-.1,y2=1.25;	
	printf("set label \"+\" at 0,.296296296296296 \n"); 
	printf("h(x)=(((x-0)<-t)+((x-0)>-t)*(2./3*(1-(x-0)/(2*t)))**2)*(((x-0)<2*t)) \n");

	 
	printf("p[%lf:%lf][%lf:%lf] h(x) t'h exact','-' u 1:2 t'Q'  w l,'-'  t'h' w lp,'-'  t'Z' w l linec -1\n ",
			   -nx*dx/2,nx*dx/2,y1,y2); 
		for (i=0; i<=nx;i++)
		{
			printf("%lf %lf \n",x[i],Q[i]);}
		printf("e \n");
		for (i=0; i<=nx;i++)
		{
			printf("%lf %lf \n",x[i],h[i]+Z);}
		printf("e \n");
		for (i=0; i<=nx;i++)
		{
			printf("%lf %lf \n",x[i],Z);}
		printf("e \n");
*/
   
      }
     }  
     free(x);
     free(fp);
     free(fd); 
     free(un);
     free(hn);  
    return 0;
}