#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <iomanip>
#include <fstream>

const int nx=29;
const int ny=29;
const int nz=21;

float lx=9.9, ly=4.9, lz=2.0;
double dx=0.1, dy=0.1, dz=0.1, dt=0.001;
double cx=7.0, cy=4.0, h=0.8, r1=0.5, r2=0.2;
double potx[nz][ny][nx], poty[nz][ny][nx], potz[nz][ny][nx];
double u[nz][ny][nx], v[nz][ny][nx], w[nz][ny][nx];
double curlx[nz][ny][nx], curly[nz][ny][nx], curlz[nz][ny][nx];
double tempx[nz][ny][nx], tempy[nz][ny][nx], tempz[nz][ny][nx];
double nu=0.2, err=0.001;
int op=40;      //replace by adequate checking of differences!!!

using namespace std;

int main()
{
    int nx= int(lx/dx+1);
    int ny= int(ly/dy+1);
    int nz= int(lz/dz+1);

    for(int k=0;k<nz;k++)
        for(int j=0;j<ny;j++)
            for(int i=0;i<nx;i++) {
                double x= i * dx;
                double y= j * dy;
                double z= k *dz;
                u[k][j][i]=0.0;
                v[k][j][i]=0.0;
                w[k][j][i]=0.0;
                potx[k][j][i]=0.0;
                poty[k][j][i]=0.0;
                potz[k][j][i]=0.0;
                if(y==0)
                    v[k][j][i]=20.0*z*z;

                curlx[k][j][i]=0.0;
                curly[k][j][i]=0.0;
                curlz[k][j][i]=0.0;
                if((y==0)&&(k<nz-1)) {
                    curlx[k][j][i]=-(v[k+1][j][i]-v[k][j][i])/dz;
                    curlz[k][j][i]=-v[k][j][i]/dx;
                }
            }

    int p=0;
    while(p<op) {

        for(int k=1;k<nz-1;k++)
            for(int j=1;j<ny-1;j++)
                for(int i=1;i<nx-1;i++) {
                    tempx[k][j][i]=curlx[k][j][i]-u[k][j][i]*dt*(curlx[k][j][i+1]-curlx[k][j][i])/dx-v[k][j][i]*dt*(curlx[k][j+1][i]-curlx[k][j][i])/dy-w[k][j][i]*dt*(curlx[k+1][j][i]-curlx[k][j][i])/dz+
                        curlx[k][j][i]*dt*(u[k][j][i+1]-u[k][j][i])/dx+curly[k][j][i]*dt*(u[k][j+1][i]-u[k][j][i])/dy+curlz[k][j][i]*dt*(u[k+1][j][i]-u[k][j][i])/dz+
                        nu*dt*((curlx[k][j][i+1]-2*curlx[k][j][i]+curlx[k][j][i-1])/(dx*dx)+(curlx[k][j+1][i]-2*curlx[k][j][i]+curlx[k][j-1][i])/(dy*dy)+(curlx[k+1][j][i]-2*curlx[k][j][i]+curlx[k-1][j][i])/(dy*dy));

                    tempy[k][j][i]=curly[k][j][i]-u[k][j][i]*dt*(curly[k][j][i+1]-curly[k][j][i])/dx-v[k][j][i]*dt*(curly[k][j+1][i]-curly[k][j][i])/dy-w[k][j][i]*dt*(curly[k+1][j][i]-curly[k][j][i])/dz+
                        curlx[k][j][i]*dt*(v[k][j][i+1]-v[k][j][i])/dx+curly[k][j][i]*dt*(v[k][j+1][i]-v[k][j][i])/dy+curlz[k][j][i]*dt*(v[k+1][j][i]-v[k][j][i])/dz+
                        nu*dt*((curly[k][j][i+1]-2*curly[k][j][i]+curly[k][j][i-1])/(dx*dx)+(curly[k][j+1][i]-2*curly[k][j][i]+curly[k][j-1][i])/(dy*dy)+(curly[k+1][j][i]-2*curly[k][j][i]+curly[k-1][j][i])/(dy*dy));

                    tempz[k][j][i]=curlz[k][j][i]-u[k][j][i]*dt*(curlz[k][j][i+1]-curlz[k][j][i])/dx-v[k][j][i]*dt*(curlz[k][j+1][i]-curlz[k][j][i])/dy-w[k][j][i]*dt*(curlz[k+1][j][i]-curlz[k][j][i])/dz+
                        curlx[k][j][i]*dt*(w[k][j][i+1]-w[k][j][i])/dx+curly[k][j][i]*dt*(w[k][j+1][i]-w[k][j][i])/dy+curlz[k][j][i]*dt*(w[k+1][j][i]-w[k][j][i])/dz+
                        nu*dt*((curlz[k][j][i+1]-2*curlz[k][j][i]+curlz[k][j][i-1])/(dx*dx)+(curlz[k][j+1][i]-2*curlz[k][j][i]+curlz[k][j-1][i])/(dy*dy)+(curlz[k+1][j][i]-2*curlz[k][j][i]+curlz[k-1][j][i])/(dy*dy));
                }

        for(int k=1;k<nz-1;k++)
            for(int j=1;j<ny-1;j++)
                for(int i=1;i<nx-1;i++) {
                    curlx[k][j][i]=tempx[k][j][i];
                    curly[k][j][i]=tempy[k][j][i];
                    curlz[k][j][i]=tempz[k][j][i];

                    potx[k][j][i]=0.5/(1/(dx*dx)+1/(dy*dy)+1/(dz*dz))*curlx[k][j][i];
                    poty[k][j][i]=0.5/(1/(dx*dx)+1/(dy*dy)+1/(dz*dz))*curly[k][j][i];
                    potz[k][j][i]=0.5/(1/(dx*dx)+1/(dy*dy)+1/(dz*dz))*curlz[k][j][i];
                }

        double diff=err+1.0;
        while( diff > err ) {
            diff=0.0;

            for(int k=1;k<nz-1;k++)
                for(int j=1;j<ny-1;j++)
                    for(int i=1;i<nx-1;i++) {
                        tempx[k][j][i]=0.5/(1/(dx*dx)+1/(dy*dy)+1/(dz*dz))*((potx[k][j][i+1]+potx[k][j][i-1])/(dx*dx)+(potx[k][j+1][i]+potx[k][j-1][i])/(dy*dy)+(potx[k+1][j][i]+potx[k-1][j][i])/(dz*dz)+curlx[k][j][i]);
                        if ( fabs(tempx[k][j][i]-potx[k][j][i]) > diff )
                            diff=fabs(tempx[k][j][i]-potx[k][j][i]);
                    }

            for(int k=1;k<nz-1;k++)
                for(int j=1;j<ny-1;j++)
                    for(int i=1;i<nx-1;i++)
                        potx[k][j][i]=tempx[k][j][i];

        }

        diff=err+1.0;
        while( diff > err ) {
            diff=0.0;
            for(int k=1;k<nz-1;k++) {
                for(int j=1;j<ny-1;j++) {
                    for(int i=1;i<nx-1;i++) {
                        tempy[k][j][i]=0.5/(1/(dx*dx)+1/(dy*dy)+1/(dz*dz))*((poty[k][j][i+1]+poty[k][j][i-1])/(dx*dx)+(poty[k][j+1][i]+poty[k][j-1][i])/(dy*dy)+(poty[k+1][j][i]+poty[k-1][j][i])/(dz*dz)+curly[k][j][i]);
                        if( fabs(tempy[k][j][i]-poty[k][j][i]) > diff )
                            diff = fabs(tempy[k][j][i]-poty[k][j][i]);
                    }
                }
            }

            for(int k=1;k<nz-1;k++)
                for(int j=1;j<ny-1;j++)
                    for(int i=1;i<nx-1;i++)
                        poty[k][j][i] = tempy[k][j][i];
        }

        diff=err+1.0;
        while( diff > err ) {
            diff=0.0;
            for(int k=1;k<nz-1;k++)
                for(int j=1;j<ny-1;j++)
                    for(int i=1;i<nx-1;i++){
                        tempz[k][j][i]=0.5/(1/(dx*dx)+1/(dy*dy)+1/(dz*dz))*((potz[k][j][i+1]+potz[k][j][i-1])/(dx*dx)+(potz[k][j+1][i]+potz[k][j-1][i])/(dy*dy)+(potz[k+1][j][i]+potz[k-1][j][i])/(dz*dz)+curlz[k][j][i]);
                        if( fabs(tempz[k][j][i]-potz[k][j][i]) > diff )
                            diff = fabs(tempz[k][j][i]-potz[k][j][i]);
                    }

            for(int k=1;k<nz-1;k++)
                for(int j=1;j<ny-1;j++)
                    for(int i=1;i<nx-1;i++)
                        potz[k][j][i] = tempz[k][j][i];
        }

        for(int k=1;k<nz-1;k++)
            for(int j=1;j<ny-1;j++)
                for(int i=1;i<nx-1;i++) {
                    u[k][j][i]=(potz[k][j+1][i]-potz[k][j][i])/dy-(poty[k+1][j][i]-poty[k][j][i])/dz;
                    v[k][j][i]=(potx[k+1][j][i]-potx[k][j][i])/dz-(potz[k][j][i+1]-potz[k][j][i])/dx;
                    w[k][j][i]=(poty[k][j][i+1]-poty[k][j][i])/dx-(potx[k][j+1][i]-potx[k][j][i])/dy;
                }

        for(int k=0;k<nz;k++)
            for(int j=0;j<ny;j++)
                for(int i=0;i<nx;i++) {
                    double x= i * dx;
                    double y= j * dy;
                    double z= k *dz;
                    if(z>h)
                        continue;
                    double r=r2+(r1-r2)*(h-z)/h;
                    if((x<cx-r)||(x>cx+r))
                        continue;
                    if((y<cx-sqrt(r*r-(x-cx)*(x-cx)))||(y>cx+sqrt(r*r-(x-cx)*(x-cx))))
                        continue;
                    u[k][j][i]=0.0;
                    v[k][j][i]=0.0;
                    w[k][j][i]=0.0;
                    potx[k][j][i]=0.0;
                    poty[k][j][i]=0.0;
                    potz[k][j][i]=0.0;
                }

        p++;
    }

    ofstream XZ("XZ.txt");
    XZ << nx << "\t" << nz << "\t" << endl;
    for (int i=0; i<nx; i++){
        for (int j=0;j<nz;j++) {
          XZ << fixed << u[j][10][i] << "\t";
          XZ << fixed << w[j][10][i] << "\t";
        }
        XZ<<endl;
    }

    ofstream YZ("YZ.txt");
    YZ << ny << "\t" << nz << "\t" << endl;
    for (int i=0;i<ny;i++) {
        for (int j=0;j<nz;j++) {
            YZ << fixed << v[j][i][10] << "\t";
            YZ << fixed << w[j][i][10] << "\t";
        }
        YZ<<endl;
    }

    ofstream YX("YX.txt");
    YX << ny << "\t" << nx << "\t" << endl;
    for (int i=0;i<ny;i++) {
        for (int j=0;j<nx;j++) {
            YX << fixed << u[5][i][j] << "\t";
            YX << fixed << v[5][i][j] << "\t";
        }
        YX<<endl;
    }

    return 0;
}