//********************************
//****  Elementary Functions  ****
//********************************
double abs(double x) {
  if(x<0)
    return (-x);
  return(x);
}
double  max(double  a,double  b){
  if(a>b)
    return(a);
  return(b);
}

double  min(double  a,double  b){
  if(a>b)
    return(b);
  return(a);
}

double  max(double U[],int n){
  int j;
  double m=U[0];

  for(j=1;j<mesh_max;++j)
  if(abs(U[j])>m)
    m=abs(U[j]);

  return(m);
}

double  max(double U[][2],int n){
  int j;
  double m=U[0][n];

  for(j=1;j<mesh_max;++j)
  if(abs(U[j][n])>m)
    m=abs(U[j][n]);

  return(m);
}

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);
}

//*****************************
//**** Universal Constants ****
//*****************************

const double PI=3.14159;

//*************************************
//**** Constants to be used in FDM ****
//*************************************





//**********************************************
//**** Functions to give Initial conditions ****
//**** Returns Nothing. Chnages the array   ****
//**********************************************
void initialize_ROLL(double U[][2]){
  int j;
  double x;

  x=L1-h;
  for(j=0;j<mesh_max;++j){
    x=x+h;
   
    if     (x< L1+h/2) U[j][0]=0;
    else if(x<   -h/2) U[j][0]=x-L1 ;
    else if(x<    h/2) U[j][0]= 0;
    else if(x< L2-h/2) U[j][0]=x-L2;
    else               U[j][0]=0;
  }
  for(j=0;j<mesh_max;++j) U[j][1]=U[j][0];
}
 
void initialize_E(double U[][2],double a, double b){
  int j; 
  double x=L1-h,s;
 
  for(j=0;j<mesh_max;++j){
    x=x+h;
    
    if     (x<L1+h/2)        U[j][0]=0;
    else if(x< (L1+L2)/2-h/2) U[j][0]=a;
    else if(x< (L1+L2)/2+h/2) U[j][0]=0;
    else if(x< L2-h/2)        U[j][0]=b;
    else                      U[j][0]=0;
  }
}

void initialize_E(double U[][2],double a, double b,double c,double d){
  int j; 
  double x=L1-h,s;
 
  for(j=0;j<mesh_max;++j){
    x=x+h;
    
    if     (x<L1+h/2) U[j][0]=0;
    else if(x<L1/2-h/2) U[j][0]=a;
    else if(x<L1/2+h/2) U[j][0]=0;
    else if(x<  0-h/2) U[j][0]=b;
    else if(x<  0+h/2) U[j][0]=0;
    else if(x<L2/2-h/2) U[j][0]=c;
    else if(x<L2/2+h/2) U[j][0]=0;
    else if(x<L2-h/2) U[j][0]=d;
    else if(x<L2+h/2) U[j][0]=0;
    else               U[j][0]=0;
  }
}