#include <iostream.h>
#include <fstream.h>
#include <stream.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
//************************************************
//**** A function to be used to decide constants *
//************************************************
int Integer(double  a){ int i; 
     if(a>0) { for(i=0;i<a||i==a;++i); return(i-1); }
     if(a<0) { for(i=0;i>a;--i); return(i); } return(0); }
const double pi= 3.14159                     ;// time t<T
double imax(double a,double b){
  if(a>b)
     return(a);
  return(b);
}
//*************************************
//**** Constants to be used in FDM ****
//*************************************
const double T= 70;//10.1*log(10)                      ;// maximum time t<T
const double into=1/log(10)           ;
           //to control output.
           // into number of printouts for each time interval of length 1
const double L1=-1 ;// Left end of the interval
const double L2= 1 ;// Right end of the interval
const double epsilon =0.02;   
const double h=0.01           ;                        // mesh size is fixed
const double lambda = 6/log(10)    ;//Numerical Speed and time step k depend
const double k=h/lambda             ;//on the maximum of the solution. 
const double interval=L2-L1           ;
const int skip=1;//Integer(imax(interval/h/200,1.0)); //to control out put
const int   mesh_max=Integer(interval/h);

//**************************************************************
//**** I have all the head files and constants in "head.h". ****
//**** In the main function we have I/O's and recursive     ****
//**** calls of FDM elementary functions.                   ****
//**************************************************************
#include "head.h2"

//*****************************************************************************
//**** Different Kinds of Schemes
//**** "conse.h"  - traditional conseting method
//**** "homo.h"   - method based on higher dimensional homogeneous
//****              conservation law
//*****************************************************************************

#include "schemes.h2"

//************************************************************
//**** In the main function we havee I/O's and recursive *****
//**** calls of FMD elementary functions.                *****
//************************************************************
//**** Variables defined in "head.h"                      **** 
//**** space domain (L1,L2)       h:space mesh size       ****
//**** k:time mesh size           S:time domain (0,S)     ****
//**** mesh_max:Number of grids                           ****
//************************************************************
int main(){
int n,j,m,i;
double U1[mesh_max][2],U4[mesh_max][2];
double A[mesh_max],B[mesh_max],U[mesh_max];
double R[mesh_max][2],t,x,l,r,integral,oldmax,newmax;
double tmp=-1,d1,d2,d3; 

    //following two initial functions are in head.h2
initialize_E(U1,1,-1);              // Initial Condition for numerical solution
initialize_ROLL(U4);                // Roll wave for comparison

  t=0;
  for(n=0;t<T;++n){                             //Loop for time t
      t=t+k;

    //following loop have to start from j=1 since -1%10=-1, not 9
    // the reason I use operator "%" is to control the boundary.
    for(j=1 ;j<mesh_max+1;++j){               //Loop for space variables
            //"_split" implies (3.7)
            //"_non1" implies  (3.4)
            //"_non2" implies  (3.10)
          U1[j%mesh_max][(n+1)%2]=godunov_non2(h,k,U1[(j-1)%mesh_max][n%2],U1[j%mesh_max][n%2],U1[(j+1)%mesh_max][n%2]);
    } // end of for(j=1 ;j<mesh_max+1;++j)


//following is for out-put control
    for(m=0;m<into*T;++m)
      if(m<(t+0.000001)*into&&(t+0.000001)*into<m+k*into){//checking output time
//following is to print our the solution
//         for(j=0;j<mesh_max;j=j+skip)
//           cout<<form("%f ",U1[j][0]);
//         cout<<endl;

//following is to print M(t) and L^1 convergence
//functions M(U1), L1norm(U1,U4) are in schemes.h2
      cout<<form("\\+&%f  &%e   &%e \\cr",t,M(U1),L1norm(U1,U4))<<endl;
    } //end of if(m<(t+0.000001)*into&&(t+0.000001)*into<m+k*into){
  }  //end of for(n=0;t<T;++n){

  cout<<form("%c",7); // to notify the run is done
return(1);
}