#include <stdio.h>    /* Advection. */
#include <stdlib.h>   /* On considère que u=1 (vitesse) */
#include <malloc.h>   /* et que le segment en espace est [0, 1]. */
#include <math.h>


double maximum(double x, double y)
{
  if(x>=y)
    return(x);
  else 
    return(y);
}

double minimum(double x, double y)
{
  if(x<=y)
    return(x);
  else 
    return(y);
}

double minmod(double x, double y)
{
  return(maximum(0,minimum(x,y)));
}

double phiminmodmoi(double r, double cfl)
{
  if(r<0.)
    return(0.);
  else                                  // LE NÔTRE. 
    if(r<cfl/2.)
      return(r*2./cfl);
    else
      return(1.);
}

double phisuperbee(double r, double cfl)
{
  if(r<0.)
    return(0.);
  else
    if(r<0.5)
      return(2*r);
    else                               // SUPERBEE. 
      if(r<1.)      
	return(1.);
      else
	if(r<2.)
	  return(r);
	else
	  return(2.);
}

double phiminmod(double r, double cfl)
{
  if(r<0.)
    return(0.);
  else
    if(r<1.)                            // MINMOD. 
      return(r);
    else
      return(1.);
}

double phimoi(double r, double cfl)
{
  if(r<0.)
    return(0.);
  else                                  // LE NÔTRE. 
    if(r<cfl/(1-cfl))
      return(r*2./cfl);
    else
      return(2./(1.-cfl));
}

double phimoispecial(double deltau, double gradient, double r, double cfl)
{
double s1,s2,s3;

  if(r<0.)
    s2=0.;
  else                                  // LE NÔTRE. 
    if(r<cfl/(1-cfl))
      s2=r*2./cfl;
    else
      s2=2./(1.-cfl);

  s1=2.*gradient/((1-cfl)*(s2+0.0000000000000001)
*sqrt(deltau*deltau+0.0000000001)  );

    if (s1<1.) 
    { s2=s1*s2;
    //  printf("s1= %lf , s2= %lf\n",s1,s2);
    }
  //  s2=s2/2.;
  return(s2);

}

double phi1(double r, double cfl)
{                                     // LAX-WENDROFF (SANS LIMITEUR). 
  return(1.); 
}

double phiWB(double r, double cfl)
{ 
  if(r<=0)
    return(0.);                         // WARMING & BEAM. 
  else
    return(r);
}

double phi0(double r, double cfl)
{                                     // DÉCENTRÉ AMONT. 
  return(0.); 
}

double phimoimodifie(double r, double cfl)
{ 
  if(r<=0.)
    return(0.);
  else 
    if(r>=0.9)
      if(r<=1.1)
	return(1.+100.*(r-1)*(r-1));
      else
	return(phimoi(r,cfl));
    else
      return(phimoi(r,cfl));
} 

double phialeatoire(double r, double cfl)
{ 
  double ale;
  ale=maximum(0.,minimum(2.*r/cfl,2./(1.-cfl)));
  return(rand()*ale/(RAND_MAX+1.));
}

double phi2minmod(double r, double cfl)
{
  if(r<0.)
    return(0.);
  else
    if(r<1.)                            // 2*MINMOD. 
      return(2.*r);
    else
      return(2.);
}


double phiessai(double r, double cfl)
{
  if(r<=1.-cfl)
    return(phimoi(r,cfl));
  else  
    return(phimoi(1./r,cfl));
}




void main()
{
  FILE* fp = NULL;
  FILE* l2 = NULL;
  FILE* fpini   = NULL;    /* Fichiers de donnees */
FILE* fichierperio = NULL;

  
  const int kb=40000;  /* deltax=1/ka */
  
char nomfichierperio[255];



  double rho[kb],rhoini[kb],rhotemp[kb],fluxrho[kb],
    deltax,CFL,massetotale,r1,r,
    deltataleatoire,deltatmax,t;
  int k,n,ka;       /* k indice d'espace, n indice d'itération. */
  int dk,uk,ukk,k1,schema,condini;
  char   dummy[255];
  double deltat;
  double tempsfinal;
  double gradient;
  double scal,scal2,normespec;

  double normel1,normel2,normelinfty;
  /**************************************************************************/
  /*                  Initialisation des variables par défaut.              */
  /**************************************************************************/
  
  /* Ouverture du fichier de données. */
  
  fpini=fopen("transport.ini","r");
  

  
  /* Nombre de points du maillage. */
  fgets(dummy,255,fpini);
  fscanf(fpini,"%d\n",&ka);
  
 
  /* Temps final. */
  fgets(dummy,255,fpini);
  fscanf(fpini,"%lf\n",&tempsfinal);
  
  /* Condition CFL. */
  fgets(dummy,255,fpini);
  fscanf(fpini,"%lf\n",&CFL);
  
  /* Schéma. */
  fgets(dummy,255,fpini);
  fscanf(fpini,"%d\n",&schema);
  

  /* Condition initiale     . */
  fgets(dummy,255,fpini);
  fscanf(fpini,"%d\n",&condini);


  
  /* Fermeture du fichier de donnees */
  fclose(fpini);



  //  printf("longueur max de gradient \n");
  // scanf("%lf",&gradient);
  // printf("longueur max de gradient = %lf \n", gradient);

  /**************************************************************************/
  /*                          Choix des données.                            */
  /**************************************************************************/
  
  do
    {
      printf("\n\n");
      printf("*************************************************************************\n");
      printf("** Code de calcul  Advection                                           **\n");
      printf("*************************************************************************\n");
      printf("    a) Nombre de points du maillage                            : %d\n",ka);
      printf("    b) Temps final                                             : %20.15f\n",tempsfinal);
      printf("    c) Condition CFL                                           : %20.15f\n",CFL);
      printf("    d) Schema (1=up,2=minmod,3=SB,4=UB,5=LW,6=BW, 7=03)        : %d\n",schema);
      printf("    e) Initialisation (1=Rie,2=Rect,3=Ram,4=Dir,5=Gaus,6=Cos) : %d\n",condini);
      printf("*************************************************************************\n");
      printf("  Lancer le code............ [ESPACE]\n");
      printf("  Modifier les parametres... [lettre]\n");

  
  
        printf("  \nChoix --> ");
      gets(dummy);
      
      switch(dummy[0])
	{case 'a':
  printf("  Nombre de mailles : ");
  gets(dummy);
  sscanf(dummy,"%d",&ka);
  break;



  case 'b':
  printf("  Temps final : ");
  gets(dummy);
  sscanf(dummy,"%lf",&tempsfinal);
  break;

  case 'c':
  printf("  Condition CFL : ");
  gets(dummy);
  sscanf(dummy,"%lf",&CFL);
  break;

  case 'd':
  printf("  Schema  (1=up,2=minmod,3=SB,4=UB,5=LW) : ");
  gets(dummy);
  sscanf(dummy,"%d",&schema);
  break;

  case 'e':
  printf("  Condition initiale  (1=Riemann, 2=Rectangle, 3=Rampe, 4=Dirac, 5=Gaussienne), 6=Cosinus : ");
  gets(dummy);
  sscanf(dummy,"%d",&condini);
  break;
	}
    }
  while(dummy[0]!=' ');
  
  /* Sauvegarde du fichier d'initialisation par défaut. */
  fpini=fopen("transport.ini","w");
  

  
  /* Nombre de points du maillage. */
  fprintf(fpini,"Nombre de points du maillage\n%d\n",ka);
  
 
    /* Temps final. */
  fprintf(fpini,"Temps final\n%20.15lf\n",tempsfinal);

  
  /* Condition CFL. */
  fprintf(fpini,"CFL\n%20.15lf\n",CFL);
  
  /* Schéma. */
  fprintf(fpini,"Schéma\n%d\n",schema);
  
 /* Condition iniaitale. */
  fprintf(fpini,"Condition initiale\n%d\n",condini);
  

  
  /* Fermeture du fichier d'initialisation par défaut. */
  fclose(fpini);




  deltax=1./(1.*ka);
  deltat=CFL*deltax;
  

  gradient=gradient*deltax;
   


  n=0;
  t=0.;
  deltatmax=deltat;
  
  
  if (condini==1)
    {
      for(k=0;k<ka/2;k++)
	rho[k]=1.0;    // Condition initiale pour le problème de Riemann.  
      for(k=ka/2;k<ka;k++)
	rho[k]=0.;
    } 
  
  
  
  else if (condini==2)
    {
      for(k=0;k<2*ka/10;k++)
	rho[k]=0.;     
      for(k=2*ka/10;k<4*ka/10;k++)  // Rectangle  
	rho[k]=1.0;  
      for(k=4*ka/10;k<ka;k++)  
	rho[k]=0.; 
      
    }
  
  
  else if (condini==3)
    {
      for(k=0;k<3*ka/10;k++)
	rho[k]=1.-10.*k/(3.*ka);    // Condition initiale rampe.  
      for(k=3*ka/10;k<7*ka/10;k++)
	rho[k]=0.; 
      for(k=7*ka/10;k<ka;k++)
	rho[k]=(k-7.*ka/10)*10/(3.*ka); 
      
    }
  
  
  
  else if (condini==4)
    {
      for(k=0;k<ka;k++)   
	rho[k]=0.;    // Condition initiale Dirac.
      rho[ka/2]=1./deltax; 
      rho[ka/2+1]=1.;
      rho[ka/2+2]=2.;
      rho[ka/2+3]=2.;
    }
  
  
  else if (condini==5)
    {
      for(k=0;k<ka;k++)           // Condition initiale gaussienne. 
	rho[k]=exp(-50.*(1.*k/ka-1./2.)*(1.*k/ka-1./2.));
    }
  
  else  if (condini==6)
    {
      for(k=0;k<ka;k++)
      rho[k]=-cos(k*2.*3.14159/ka); // Condition initiale cosinus.
    }
 else  if (condini==8)
    {
      for(k=0;k<ka;k++)
	rho[k]=pow(4.*k*1./ka*(1.-k*1./ka),.05); // Condition initiale cosinus.
    }
  else
    {
      for(k=0;k<ka;k++)                // Condition ponctuelle 
	rho[k]=0.;
      rho[ka/2]=1.;
    }
  
  
  
   for (k=0;k<ka;k++)
     {
       rhoini[k]=rho[k];   //stockage condition initiale
     }
  
  
      l2=fopen("normes","w");
  
  
  srand(1);
  


  normespec=1.;

  while(t+deltat<tempsfinal)
    
    {
      

      
      for(k=0;k<ka;k++)
	{dk=(k-2+2*ka)%ka;
	uk=(k-1+2*ka)%ka;
	ukk=(k-2+2*ka)%ka;
	k1=(k+1+2*ka)%ka;
	
	
	// printf("Test %i %i %i %i\n",dk,uk,k,k1);

        r1=(rho[k]-rho[uk])/(rho[k1]-rho[k]);
	
	
	
        r=(rho[uk]-rho[dk])/(rho[k]-rho[uk]);
	
	
	if(schema==1)
	  rhotemp[k]=rho[k]-CFL*(rho[k]-rho[uk])
	    -CFL*(1-CFL)/2.
	    *(phi0(r1,CFL)*(rho[k1]-rho[k])
	      -phi0(r,CFL)*(rho[k]-rho[uk]));
        else if(schema==2)
	  rhotemp[k]=rho[k]-CFL*(rho[k]-rho[uk])
	    -CFL*(1-CFL)/2.
	    *(phiminmod(r1,1)*(rho[k1]-rho[k])
	      -phiminmod(r,1)*(rho[k]-rho[uk]));
        else if(schema==3)
	  rhotemp[k]=rho[k]-CFL*(rho[k]-rho[uk])
	    -CFL*(1-CFL)/2.
	    *(phisuperbee(r1,CFL)*(rho[k1]-rho[k])
	      -phisuperbee(r,CFL)*(rho[k]-rho[uk]));
        else if(schema==4)
	  rhotemp[k]=rho[k]-CFL*(rho[k]-rho[uk])
	    -CFL*(1-CFL)/2.
	    *(phimoi(r1,CFL)*(rho[k1]-rho[k])
	      -phimoi(r,CFL)*(rho[k]-rho[uk]));
        else if(schema==5)
	  rhotemp[k]=rho[k]-CFL*(rho[k]-rho[uk])
	    -CFL*(1-CFL)/2.
	    *(phi1(rho[k1]-rho[k],CFL)*(rho[k1]-rho[k])
	     -phi1(rho[k]-rho[uk],CFL)*(rho[k]-rho[uk]));
        else if(schema==6)
	  rhotemp[k]=rho[k]-CFL*(rho[k]-rho[uk])
	    -CFL*(1-CFL)/2.
	    *(phi1(rho[k1]-rho[k],CFL)*(rho[k]-rho[uk])
	     -phi1(rho[k]-rho[uk],CFL)*(rho[uk]-rho[ukk]));
        else if(schema==7)
	  {
	  scal=(1.+CFL)/3.;

	  rhotemp[k]=rho[k]-CFL*(rho[k]-rho[uk])
	    -CFL*(1-CFL)/2.
	    	    *(phi1(rho[k1]-rho[k],CFL)*((1.-scal)*(rho[k1]-rho[k])+scal*(rho[k]-rho[uk]))
	     -phi1(rho[k]-rho[uk],CFL)*((1.-scal)*(rho[k]-rho[uk])+scal*(rho[uk]-rho[ukk]))   )
	    ;
	  }
	else
	  // if(schema==10)
	  // rhotemp[k]=rho[k]-CFL*(rho[k]-rho[uk])
	  // -CFL*(1-CFL)/2.
	  // *(phiminmodmoi(r1,CFL)*(rho[k1]-rho[k])
	  //   -phiminmodmoi(r,CFL)*(rho[k]-rho[uk]));
	  //	else
	  rhotemp[k]=rho[k]-CFL*(rho[k]-rho[uk])
	    -CFL*(1-CFL)/2.
	    *(phi1(r1,CFL)*(rho[k1]-rho[k])
	      -phi1(r,CFL)*(rho[k]-rho[uk]));
	
	}
      for(k=0;k<ka;k++)
	{rho[k]=rhotemp[k];
	}

      scal=0.;
      scal2=0.;
      normel1=0.;
      normel2=0.;
      normelinfty=0.;

      for (k=0;k<ka;k++)
	{
	  scal=scal+deltax * rho[k]*(rho[k]);
	  scal2=scal2+deltax * sqrt( rho[k]*(rho[k]));
	  normel1=normel1+deltax*sqrt( (rho[k]-rhoini[k])*(rho[k]-rhoini[k]));
	  normel2=normel2+deltax*( (rho[k]-rhoini[k])*(rho[k]-rhoini[k]));
	  normelinfty=maximum(normelinfty, 
			  sqrt( (rho[k]-rhoini[k])*(rho[k]-rhoini[k])));
	}
      normel2=sqrt(normel2);
      fprintf(l2,"%f %f %f %f %f %f %f %f \n",t,scal2/deltax,scal/deltax,normel1,normel2,normelinfty,t*(log(normel1)-log(normespec))/deltat,log(10.));
      
      normespec=normel1;
  


        t+=deltat;
    n++;



    //////////////////////////////
    //  ecriture fichier   film  /
    /////////////////////////////	  
    //if(temps - tempsfinal*numero/(nbresul-1) >= 0.)

    /*
	    {
	    sprintf(nomfichierperio,"film.%05d",n);
	    if ( !(fichierperio=fopen(nomfichierperio,"w")) )  
	      {printf("Probleme d'ouverture ou de creation du fichier ");
	      printf("%s\n",nomfichierperio);
	      exit(-1);
	      }
	    
	    for(k=0;k<ka;k++)
	      fprintf(fichierperio,"%.20f %.20f\n",1.*k/ka,rho[k]);
	    fprintf(fichierperio,"\n");
	    fclose(fichierperio);
	        }
    */
	  /////////////////////
	  ///   fin film /////
	  ///////////////////

    }







  fclose(l2);

  
      fp=fopen("solution","w");

  for(k=0;k<ka;k++)
    fprintf(fp,"%f  %20.15lf %20.15lf\n",1.*k/ka,rho[k],rhoini[k]);

 
  printf("Terminé au bout de %i itérations\n",n);
  printf("Temps final %f \n", t);
  printf("Delta t %f, Delta x %f \n", deltat,deltax);
  printf("Erreur par rapport a la solution initiale en  norme l1,l2,linfini :  %20.15lf %20.15lf %20.15lf \n",
         normel1,normel2,normelinfty);

  fclose(fp);
  
}