#include <stdio.h>
#include <stdlib.h>
#include "math.h"
#include <math.h>
#include <string.h>
// pyl 05/12 
// 0K 01/2014
//
//  cc sv_ressaut.c ;./a.out | gnuplot
//
// 
// Lagree
//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
// Delestre These

 
    double*x=NULL,*h=NULL,*u=NULL,*Q=NULL,*U=NULL,*Z=NULL,*dZ=NULL,*Zp=NULL;
    double t,dt,tmax,dx,Fr,Cf,h0,xbosse,lbosse,abosse,h0,alpha,e; 
    int cas,nx,icl;
     
  
 void  reconsetat(double hl,double hr,double dz,double *hil,double *hir)
{ 	
    
      *hil=fmax(0.,hl-fmax(0.0,dz));
      *hir=fmax(0.,hr-fmax(0.0,-dz));
      
     // *hil=hl ;
     // *hir=hr;   
}


double FR1(double ug,double ud,double hg,double hd)
 { double c;
 //Rusanov 
   c=fmax(fabs(ug)+sqrt(hg),fabs(ud)+sqrt(hd));
   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));
   return (ug*ug*hg + hg*hg/2. + ud*ud*hd + hd*hd/2.)*0.5 - c*(hd*ud-hg*ug)*0.5;
 }
double FR3(double ug,double ud,double hg,double hd)
 {
   return fmax(fabs(ug)+sqrt(hg),fabs(ud)+sqrt(hd));
 } 
double F(double ug,double ud,double hg,double hd,double *f1,double *f2,double *cfl)
 { double c;
 //Rusanov 
   c=fmax(fabs(ug)+sqrt(hg),fabs(ud)+sqrt(hd));
   *f1=(hg*ug+hd*ud)*0.5-c*(hd-hg)*0.5;
   *f2=(ug*ug*hg + hg*hg/2. + ud*ud*hd + hd*hd/2.)*0.5 - c*(hd*ud-hg*ug)*0.5;
   *cfl=fabs(c);
 } 
 
double FHLL1(double ug,double ud,double hg,double hd)
 { double cfl,c1,c2,c3,f1,f2;
 //HLL 
     c1=fmin(ug -sqrt(hg),ud-sqrt(hd));
     c2=fmax(ug +sqrt(hg),ud+sqrt(hd));
     // c2=fmax(fabs(ug)+sqrt(hg),fabs(ud)+sqrt(hd));
     // c1=-c2;
      if (c1>=0.){
      f1=hg*ug;
      f2=hg*ug*ug+hg*hg/2;
      cfl=fabs(c2);}
      else{ 
      	if ((c1<0.)&&(0.<c2)) {
      f1=(c2*hg*ug-c1*hd*ud)/(c2-c1)+c1*c2*(hd-hg)/(c2-c1);
      f2=(c2*(hg*ug*ug+hg*hg/2)  -c1*(hd*ud*ud+hd*hd/2))/(c2-c1) +c1*c2*(hd*ud-hg*ug)/(c2-c1);
      cfl=fmax(fabs(c1),fabs(c2));
      	} else {
      f1=hd*ud;
      f2=hd*ud*ud+hd*hd/2.;
      cfl=fabs(c1);} }    
   return f1;
 }
 
double FHLL2(double ug,double ud,double hg,double hd)
 { double cfl,c1,c2,c3,f1,f2;
 //HLL 
     c1=fmin(ug -sqrt(hg),ud-sqrt(hd));
     c2=fmax(ug +sqrt(hg),ud+sqrt(hd));
     //c2=fmax(fabs(ug)+sqrt(hg),fabs(ud)+sqrt(hd));
     // c1=-c2;
 
      if (c1>=0.){
      f1=hg*ug;
      f2=hg*ug*ug+hg*hg/2;
      cfl=fabs(c2);}
      else{ 
      	if ((c1<0.)&&(0.<c2)) {
      f1=(c2*hg*ug-c1*hd*ud)/(c2-c1)+c1*c2*(hd-hg)/(c2-c1);
      f2=(c2*(hg*ug*ug+hg*hg/2)  -c1*(hd*ud*ud+hd*hd/2))/(c2-c1) +c1*c2*(hd*ud-hg*ug)/(c2-c1);
      cfl=fmax(fabs(c1),fabs(c2));
      	} else {
      f1=hd*ud;
      f2=hd*ud*ud+hd*hd/2.;
      cfl=fabs(c1);} }    
   return f2;
 }
double FHLL3(double ug,double ud,double hg,double hd)
 { double cfl,c1,c2,c3,f1,f2;
 //HLL 
     c1=fmin(ug -sqrt(hg),ud-sqrt(hd));
     c2=fmax(ug +sqrt(hg),ud+sqrt(hd));
     //     c2=fmax(fabs(ug)+sqrt(hg),fabs(ud)+sqrt(hd));
     // c1=-c2;
if (c1>=0.){
      cfl=fabs(c2);}
      else{ 
      	if ((c1<0.)&&(0.<c2)) {
        cfl=fmax(fabs(c1),fabs(c2));
      	} else {
       cfl=fabs(c1);} }     
   return cfl;
} 

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

void ecrit()
{    
	extern double  dx,t,dt,tmax; 
  	extern int nx;
	int i;
	double y1,y2,xf,hf;
	extern double *x,*h,*Z,*dZ,*Q,Fr,e;
	FILE *g;
	extern int cas;
 
	/* Saving the fields */ 

	 // sortie GNUPLOT	
	
	
	printf("t=%lf\n",t);
	printf("set xlabel \"x\" ; set title \"t=%lf  e=%lf \"\n",t,e); 
 
	y1=0,y2=(-1+sqrt(1+ 8*Fr*Fr))/2-1; ;	 
	printf("h(x)=(%lf)*(x>0)+1 \n",y2); 
	y2=fmax(2,y2+1.5);
 	 
	printf("p[%lf:%lf][%lf:%lf] h(x) t'h exact','-' u 1:2 t'Q'  w l,'-'  t'h' w d,'-'  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[i]);}
		printf("e \n");
		for (i=0; i<=nx;i++)
		{
			printf("%lf %lf \n",x[i],Z[i]);}
		printf("e \n");
	
	
	
} 
     
    
int main (int argc, const char *argv[]) {
    int  i,it=0; 
    char file;
    double*fp=NULL,*fd=NULL,*hg=NULL,*hd=NULL,*hid=NULL,*hig=NULL,*ug=NULL,*ud=NULL,*un=NULL,*hn=NULL;
    double cfl2,q,f1,f2,sign;
       
 
      dt=0.015;
     tmax=10;
     dx=0.025;
     nx=1000;
     Fr=1.5;
     
    
    
     
 

// alloc 
     
      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));
      Z= (double*)calloc(nx+1,sizeof(double));
      dZ=(double*)calloc(nx+1,sizeof(double));
      Zp=(double*)calloc(nx+1,sizeof(double)); 
      fp=(double*)calloc(nx+1,sizeof(double));
      fd=(double*)calloc(nx+1,sizeof(double)); 
      hg=(double*)calloc(nx+1,sizeof(double));
      hd=(double*)calloc(nx+1,sizeof(double)); 
      hig=(double*)calloc(nx+1,sizeof(double));
      hid=(double*)calloc(nx+1,sizeof(double)); 
      ug=(double*)calloc(nx+1,sizeof(double));
      ud=(double*)calloc(nx+1,sizeof(double)); 
      un=(double*)calloc(nx+1,sizeof(double));
      hn=(double*)calloc(nx+1,sizeof(double));
      
      
    // initialisation            
       for(i=0;i<=nx;i++)
    {  x[i]=(i-nx/2)*dx;  
	 Z[i]=0;
     h[i]=1+(((-1 + sqrt(1+8*Fr*Fr))/2)-1)*(1+tanh(x[i]))/2;
     Q[i]=Fr;
     u[i]=Q[i]/h[i];
	     }        
     
 for(i=0;i<nx;i++)
    {  dZ[i]= Z[i+1]-Z[i];	}
    dZ[nx]=0;




     
      
      
    while(t<=tmax){
      t=t+dt;
      it=it+1;
      
// variables pour la reconstruction
       for(i=0;i<=nx;i++)
    { 
      hd[i]=h[i];
      ud[i]=u[i];
      hg[i]=h[i];
      ug[i]=u[i];}

      /*


c-----la reconstruction "hydrostatique"
c-----construction d une section a gauche a l interface A_G(i-1/2) notee Ad(i-1)
c-----construction d une section a droite a l interface A_D(i-1/2) notee Ag(i)
      call reconsetat(Ar(i-1),Al(i),delRaA0(i-1),Ad(i-1),Ag(i))
 */ 
 for(i=1;i<=nx;i++)
    { 
       reconsetat(hd[i-1],hg[i],dZ[i-1],&f1,&f2);
       hid[i-1]=f1;
       hig[i]=f2;
         
      }
         
 	FILE *g = fopen("sol.OUT", "w");
	//for (i=0; i<=nx;i++) fprintf(g,"%lf %lf %lf %lf %lf \n",x[i],h[i],Q[i],Z[i],hig[i],hid[i]);
	for (i=0; i<=nx;i++) fprintf(g,"%lf %lf %lf %lf \n",x[i],h[i],Q[i],Z[i]);
	fclose(g);	
	
// construction de flux intermediaires a l'interface : i-1/2  fp:premier  fd:deuxieme

   for(i=1;i<=nx;i++)
    {    
    	  fp[i]=FHLL1(ud[i-1],ug[i],hid[i-1],hig[i]);
    	  fd[i]=FHLL2(ud[i-1],ug[i],hid[i-1],hig[i]); 
          cfl2 =FHLL3(ud[i-1],ug[i],hid[i-1],hig[i]); 
         // F(ud[i-1],ug[i],hd[i-1],hg[i],&fp[i],&fd[i],&cfl2);
    }
  
  
   if(cfl2>dx/dt){ dt = dx/cfl2/2;}
    
// boucle principale
//   U=(h,hu)
//  	Unouveau = Uancien  + dt (Flux gauche-Flux droit)
		
   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 
       -dt*( hig[i]*hig[i]/2 - hg[i]*hg[i]/2
            + hd[i]*hd[i]/2- hid[i]*hid[i]/2)/dx;
          
      
	  // frottement     
	  if(q>0)  sign= 1; 
      if(q<0)  sign=-1;
      if((q!=0)&&(Cf>0)) q = q/(1+dt*Cf/h[i]/h[i]);
	  
      un[i]=q/hn[i];}
      else{
      un[i]=0.;}
    }
		
 
			hn[0]=hn[1];
			un[0]=un[1];
		
		hn[nx]=hn[nx-1];
		un[nx]=un[nx-1];	
		//u[nx]=u[nx-1]-sqrt(h[nx-1]) + sqrt(h[nx]);
		//; 
		e=0;
 
  for(i=0;i<=nx;i++)
    { // e=e+(h[i]-hn[i])*(h[i]-hn[i]);
	   e=e + hn[i]*dx;
      h[i]=hn[i];
      u[i]=un[i];
      Q[i]=hn[i]*un[i];
    }


  if(it%150==0)ecrit();
		
      }
      
     free(x);
     free(fp);
     free(fd);
     free(hd);
     free(hg); 
    return 0;
}


//http://www.ontko.com/pub/rayo/cs35/pointers.html