//**************************************************************************
//****  Following is a sample program  named 'sample.C'.
//****  Compile the program  (Command is : g++ sample.C -lm )
//****  and run it  (Command is : a.out  )
//****  Then the result will be stored in a file named 'outputfilename'.
//**************************************************************************
#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;
const double X2=3000;
const double Y1=0;
const double Y2=3000;
const double days=1;
const double T=days*24;                        // final hour
const int MAX=121;                            // MAX=(L2-L1)/dx+1
const double nv=10.0;

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

//**************************************************************************
//**** Initial data setup
//**************************************************************************
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]=0.0;
      c[j][l][1]=0.0;
    }
  }
}
//**************************************************************************
//**** Random function which returns a number between 0 and 2PI
//**************************************************************************
double myrandomnumber(){
  return(((random()%1000000)*PI)/500000.0);
}
//**************************************************************************
//**** Wind Velocity
//**************************************************************************
double xocean(double beta){
  return(5+5*cos(beta));
}
double xland(double beta){
  return(2+2*cos(beta));
}
double U(double x,double y,double beta){
//   if(y>3000) return(0);
   if(y>2300){
     if(x<y-700) return(xland(beta));
     else        return(xocean(beta));
   } 
   if(y>2000){
     if(x<y-700||2300<x&&x<2500) return(xland(beta));
     else                        return(xocean(beta));
   }
   if(y>1600){
     if(x<700||1100<x&&x<1400||2300<x&&x<2500) return(xland(beta));
     else return(xocean(beta));
   }
   if(y>1300){
     if(x<700||1100<x&&x<1400||1600<x&&x<2500) return(xland(beta));
     else return(xocean(beta));
   }
   if(y>1200){
     if(x<700) return(xland(beta));
     else      return(xocean(beta));
   }
   if(y>1000){
     if(x<700||1400<x&&x<1600) return(xland(beta));
     else      return(xocean(beta));
   }
   if(y>600){
     if(x<700) return(xland(beta));
     else      return(xocean(beta));
   }
   if(y>0){
     if(x<y+100) return(xland(beta));
     else      return(xocean(beta));
   }
//   else return(0);
}
double yocean(double beta){
  return(5*sin(beta));
}
double yland(double beta){
  return(2*sin(beta));
}
double V(double x,double y,double beta){
//   if(y>3000) return(0);
   if(y>2300){
     if(x<y-700) return(yland(beta));
     else        return(yocean(beta));
   }
   if(y>2000){
     if(x<y-700||2300<x&&x<2500) return(yland(beta));
     else                        return(yocean(beta));
   }
   if(y>1600){
     if(x<700||1100<x&&x<1400||2300<x&&x<2500) return(yland(beta));
     else return(yocean(beta));
   }
   if(y>1300){
     if(x<700||1100<x&&x<1400||1600<x&&x<2500) return(yland(beta));
     else return(yocean(beta));
   }
   if(y>1200){
     if(x<700) return(yland(beta));
     else      return(yocean(beta));
   }
   if(y>1000){
     if(x<700||1400<x&&x<1600) return(yland(beta));
     else      return(yocean(beta));
   }
   if(y>600){
     if(x<700) return(yland(beta));
     else      return(yocean(beta));
   }
   if(y>0){
     if(x<y+100) return(yland(beta));
     else      return(yocean(beta));
   }
//   else return(0);
}

//**************************************************************************
//**** Diffusion
//**************************************************************************
double diffusion(double a,double b,double c){
   return((c-2*b+a)/dx/dx);
}
//**************************************************************************
//**** Polution Source
//**************************************************************************
double g(double x,double y){
   if(690-100<x&&x<690&&1000<y&&y<2000) return(100.0);
   else return(0.0);
}
//************************************************************
//**** In the main function we have I/O's and recursive  *****
//**** calls of FDM 
//************************************************************
void main(){
  filebuf fi;
  int j,m,l,next,now;
  double x,y,beta[MAX][MAX],beta0;
  double c[MAX][MAX][2];


  initialize(c);

  next=0;
  for(m=0;m*dt<T;++m){
    beta0=myrandomnumber();
    now=next;
    next=(now+1)%2;
    for(j=1;j<MAX-1;++j){
      for(l=1;l<MAX-1;++l){
        beta[j][l]=myrandomnumber();
      }
    }
    for(j=1;j<MAX-1;++j){
      x=X1+j*dx;
      for(l=1;l<MAX-1;++l){
        y=Y1+l*dy;
        c[j][l][next]=(c[j-1][l][now]+c[j+1][l][now]+c[j-1][l-1][now]+c[j][l-1][now]+c[j+1][l-1][now]+c[j-1][l+1][now]+c[j][l+1][now]+c[j+1][l+1][now])/8
                -0.5*(dt/dx)*(U(x+dx,y,beta[j+1][l])*c[j+1][l][now]
                             -U(x-dx,y,beta[j-1][l])*c[j-1][l][now])
                -0.5*(dt/dy)*(V(x,y+dy,beta[j][l+1])*c[j][l+1][now]
                             -V(x,y-dy,beta[j][l-1])*c[j][l-1][now])
                +g(x+dx/2,y+dy/2)*dt                                   ;
      }
    }
  } 
  ofstream output("output", ios::out);
  for(l=0;l<MAX;++l){
    for(j=0;j<MAX;++j){
       output<<form("%f ",c[j][l][0]);
    }
    output<<endl;
  }
cout<<form("%c",7);
}