//**************************************************************************
//****  Following is a sample program  named 'sample2.C'.
//****  Compile the program  (Command is : g++ sample2.C -lm )
//****  and run it  (Command is : a.out  )
//****  Then the result will be stored in a file named 'output'.
//****  This program is to solve Example 2.5.2, Friedman
//**************************************************************************

#include <iostream.h>
#include <fstream.h>
#include <stream.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>

//**************************************************************************
//**** Constants to be used in FDM
//**************************************************************************
const double X1=0.0;
const double X2=15;
const double Y1=-15;
const double Y2=0;
const double T=5.0;                            // final time
const int MAX=41;                            // MAX=(L2-L1)/dx+1
const double nv=1.5;

const double dx=(X2-X1)/(MAX-1);                // dx
const double dy=(Y2-Y1)/(MAX-1);                // dy
const double dt=dx/nv;                          // dt

//**************************************************************************
//**** Initial data setup
//**************************************************************************
double initial_function(double x,double y){
  double r;
  r=sqrt((x-5)*(x-5)+(y+10)*(y+10));
  
  if(r>4) return(0.0);
  else    return(50*(1+cos(PI*r/4.0)));
}
void initialize(double c[][MAX][2]){
  int j,l;
  double x,y;

  for(j=0;j<MAX;++j){
    x=X1+j*dx;
    for(l=0;l<MAX;++l){
      y=Y1+l*dy;
      c[j][l][0]=initial_function(x,y);
      c[j][l][1]=0.0;
    }
  }
}

//************************************************************
//**** In the main function we have I/O's and recursive  *****
//**** calls of FDM 
//************************************************************
void main(){
  filebuf fi;
  int j,l,m;
  double c[MAX][MAX][2];
  double x,y;

  ofstream output("output", ios::out);

  initialize(c);

  for(m=0;m*dt<T;++m){
      for(j=1;j<MAX-1;++j){
        x=X1+j*dx;
        for(l=1;l<MAX-1;++l){
          y=Y1+l*dx;
          c[j][l][(m+1)%2]=c[j][l][m%2]
                     +y/sqrt(x*x+y*y)*(dt/dx)*(c[j][l][m%2]-c[j-1][l][m%2])
                     -x/sqrt(x*x+y*y)*(dt/dy)*(c[j][l][m%2]-c[j][l-1][m%2]);
        }
    }
  }

  for(j=0;j<MAX-1;++j){
    for(l=0;l<MAX-1;++l){
       output<<form("%f ",c[j][l][0]);
    }
    output<<endl;
  }
}