#include <stdio.h>
#include <new.h>
#include "params.h"

class simuland
  {
    private:
      float dC_dy(int x,int y,int z);
      float dC_dz(int x,int y,int z);
      float far ***C,far ***NextC;
    public:
      float dx,dy,dz,dt;
      int   nx,ny,nz;
      float dC_dx(int x,int y,int z);
      float dC2_dx(int x,int y,int z);
      float dC2_dy(int x,int y,int z);
      float dC2_dz(int x,int y,int z);
      void Start();
      void Step();
      void End();
  };

float simuland::dC_dx(int x,int y,int z)
  {
    if (x>0) return (C[x][y][z]-C[x-1][y][z])/dx;
    else     return (C[x+1][y][z]-C[x][y][z])/dx;
  }

float simuland::dC_dy(int x,int y,int z)
  {
    if (y>0) return (C[x][y][z]-C[x][y-1][z])/dy;
    else     return (C[x][y+1][z]-C[x][y][z])/dy;
  }

float simuland::dC_dz(int x,int y,int z)
  {
    if (z>0) return (C[x][y][z]-C[x][y][z-1])/dz;
    else     return (C[x][y][z+1]-C[x][y][z])/dz;
  }

float simuland::dC2_dx(int x,int y,int z)
  {
    if (x>0) return (dC_dx(x,y,z)-dC_dx(x-1,y,z))/dx;
    else     return (dC_dx(x+1,y,z)-dC_dx(x,y,z))/dx;
  }

float simuland::dC2_dy(int x,int y,int z)
  {
    if (y>0) return (dC_dy(x,y,z)-dC_dy(x,y-1,z))/dy;
    else     return (dC_dy(x,y+1,z)-dC_dy(x,y,z))/dy;
  }

float simuland::dC2_dz(int x,int y,int z)
  {
    if (z>0) return (dC_dz(x,y,z)-dC_dz(x,y,z-1))/dz;
    else     return (dC_dz(x,y,z+1)-dC_dz(x,y,z))/dz;
  }

const float dC=0.5;
const float NumFloat=100000.0;

void simuland::Start()
  {
    int i,j,k;
    t=0.0;
    nx=100;
    ny=10;
    nz=5;
    dx=x/(float)nx;
    dy=y/(float)ny;
    dz=z/(float)nz;
    // Lambda currently left out of dt estimation
    dt=dC/(-(u*dC)/(Rf*dx)+
           (Dx*dC*dC)/(Rf*dx*dx)+
           (Dy*dC*dC)/(Rf*dy*dy)+
           (Dz*dC*dC)/(Rf*dz*dz));
    NextC=new float **[nx];
    C=new float **[nx];
    for (i=0;i<nx;i++)
      {
        NextC[i]=new float *[ny];
        C[i]=new float *[ny];
        for (j=0;j<ny;j++)
          {
            NextC[i][j]=new float[nz];
            C[i][j]=new float[nz];
            for (k=0;k<nz;k++)
              {
                C[i][j][k]=0.0;
                NextC[i][j][k]=0.0;
              }
          }
      }
  }

void simuland::Step()
  {
    int i,j,k;
    for (i=0;i<nx;i++)
      for (j=0;j<ny;j++)
        for (k=0;k<nz;k++)
          NextC[i][j][k]=C[i][j][k]+dt*(
                         Dx*dC2_dx(i,j,k)/Rf
                        +Dy*dC2_dy(i,j,k)/Rf
                        +Dz*dC2_dz(i,j,k)/Rf
                        -Lambda*C[i][j][k]
                        -u*dC_dx(i,j,k)/Rf);
    for (i=0;i<nx;i++)
      for (j=0;j<ny;j++)
        for (k=0;k<nz;k++)
          C[i][j][k]=NextC[i][j][k];
    t+=dt;
  }

void simuland::End()
  {
    int i,j;
    for (i=0;i<nx;i++)
      {
        for (j=0;j<ny;j++)
          {
            delete[] NextC[i][j];
            delete[] C[i][j];
          }
        delete [] NextC[i];
        delete [] C[i];
      }
    delete[] NextC;
    delete[] C;
  }

void main()
  {
    simuland sim;
    sim.u=0.01;
    sim.Dx=100.0; sim.Dy=10.0; sim.Dz=1.0;
    sim.Lambda=0.0001;
    sim.Rf=1;
    sim.x=100.0; sim.y=10.0; sim.z=1.0;
    sim.Start();
    while (sim.t<10000)
      {
        printf("Time %f\r",sim.t);
        sim.Step();
      }
    sim.End();
  }