// This is not my code.
// I have no idea how it works.
// It probably requires paid software.
// I'm planning my own implementation in Python.
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <cmath>
#include "TECIO.h"

#ifndef NULL
#define NULL 0
#endif

enum FileType { FULL = 0, GRID = 1, SOLUTION = 2 };

#define BD 20
#define BX 50
#define BY 10
#define NX 2000
#define NY 200
#define KT 10
#define AI -2./3.
#define DT 0.00001
#define RX 100
#define RY 40

#define IY 40
#define IX 80

#define Re 400
#define GR 0.0001
#define T 6000

double dxx[NX], dyy[NY], xi[NX][NY], yi[NX][NY], a[NX*NY], b[NX*NY], c[NX*NY];

double dx[NX][NY], dy[NX][NY];




void PRINT_PLT(void) {
	double SolTime;
	INTEGER4 Debug, I, J, III, DIsDouble, VIsDouble, IMax, JMax, KMax, ZoneType, StrandID, ParentZn, IsBlock;
	INTEGER4 ICellMax, JCellMax, KCellMax, NFConns, FNMode, ShrConn, FileType;

	Debug = 1;
	VIsDouble = 1;
	DIsDouble = 1;
	IMax = NY;
	JMax = NX;
	KMax = 1;
	ZoneType = 0;      /* Ordered */
	SolTime = 360.0;
	StrandID = 0;     /* StaticZone */
	ParentZn = 0;      /* No Parent */
	IsBlock = 1;      /* Block */
	ICellMax = 0;
	JCellMax = 0;
	KCellMax = 0;
	NFConns = 0;
	FNMode = 0;
	ShrConn = 0;
	FileType = FULL;

	/*
	* Open the file and write the tecplot datafile
	* header information
	*/
	I = TECINI112((char*)"SIMPLE DATASET",
		(char*)"X Y U V",
		(char*)"UV.plt",
		(char*)".",
		&FileType,
		&Debug,
		&VIsDouble);

	/*
	* Write the zone header information.
	*/
	I = TECZNE112((char*)"Simple Zone",
		&ZoneType,
		&IMax,
		&JMax,
		&KMax,
		&ICellMax,
		&JCellMax,
		&KCellMax,
		&SolTime,
		&StrandID,
		&ParentZn,
		&IsBlock,
		&NFConns,
		&FNMode,
		0,              /* TotalNumFaceNodes */
		0,              /* NumConnectedBoundaryFaces */
		0,              /* TotalNumBoundaryConnections */
		NULL,           /* PassiveVarList */
		NULL,           /* ValueLocation = Nodal */
		NULL,           /* SharVarFromZone */
		&ShrConn);

	/*
	* Write out the field data.
	*/
	III = IMax * JMax;
	I = TECDAT112(&III, &dx[0][0], &DIsDouble);
	I = TECDAT112(&III, &dy[0][0], &DIsDouble);
	I = TECDAT112(&III, &xi[0][0], &DIsDouble);
	I = TECDAT112(&III, &yi[0][0], &DIsDouble);

	I = TECEND112();


}


void PRINT_VECT_VTK(double U[NX][NY], double V[NX][NY], double *x, double *y) {
	int i, j;

	FILE *vect_field;

	/* File grid.vtk*/
	if ((vect_field = fopen("a.vtk", "w")) == NULL)
		printf("The grid file could not be open\n");
	else{
		fprintf(vect_field, "%d\n", NX);
		fprintf(vect_field, "%d\n", NY);

		for (i = 0; i <= NX - 1; i++)
			fprintf(vect_field, "%f\n", x[i]);

		for (j = 0; j <= NY - 1; j++)
			fprintf(vect_field, "%f\n", y[j]);

		for (i = 0; i <= NY - 1; i++)
			for (j = 0; j <= NX - 1; j++){
			fprintf(vect_field, "%f\t %f\t\n", U[j][i], V[j][i]);
			}
	}
	fclose(vect_field);

}

void PRINT_VTK(double U[NX][NY], double V[NX][NY], double *x, double *y) {
	int i, j;

	FILE *vect_field;


	/* File grid.vtk*/
	if ((vect_field = fopen("b.vtk", "w")) == NULL)
		printf("The grid file could not be open\n");
	else{
		fprintf(vect_field, "# vtk DataFile Version 2.0\n");
		fprintf(vect_field, "Sample rectilinear grid\n");
		fprintf(vect_field, "ASCII\n");
		fprintf(vect_field, "DATASET RECTILINEAR_GRID\n");
		fprintf(vect_field, "DIMENSIONS %d %d \n", NX, NY);
		fprintf(vect_field, "X_COORDINATES %d double\n", NX);
		for (i = 0; i <= NX - 1; i++)
			fprintf(vect_field, "%f\n", x[i]);
		fprintf(vect_field, "Y_COORDINATES %d double\n", NY);
		for (j = 0; j <= NY - 1; j++)
			fprintf(vect_field, "%f\n", y[j]);
//		fprintf(vect_field, "Z_COORDINATES %d double\n", nz);
//		fprintf(vect_field, "%f\n", z);
		fprintf(vect_field, "POINT_DATA %d\n", (NX)*(NY));
		fprintf(vect_field, "VECTORS vectors double\n");
		for (i = 0; i <= NY - 1; i++)
			for (j = 0; j <= NX - 1; j++){
				fprintf(vect_field, "%f\t %f\t %f\n", U[j][i], V[j][i], 0.0);
			}
	}
	fclose(vect_field);

}


void inixh(double *a)
{
	int i, j, k;
	for (i = 0; i < NX; i++) {
		for (j = 0; j < NY; j++) {

			k = i *NY + j;
			if (i<RX && j>RY)
				a[k] = 1.;
			else
				a[k] = 0.;
		}
	}
}
void inih(double *a)
{
	int i, j, k;

	for (i = 0; i < NX; i++) {
		for (j = 0; j < NY; j++) {
			k = i *NY + j;
			a[k] = 0.;
		}
	}
}
void USh(double *a, double *b, double *c, double *d)
{
	double u1;
	int i, j, k;
	for (i = 0; i < NX; i++) {
		for (j = 0; j < NY; j++) {
			k = i *NY + j;

			if (i < NX - 2 && i>0 && j < NY - 1 && j>0 && ((i <= RX && j <= RY) != 1)){
				u1 = -((pow(a[k + NY], 2) - pow(a[k - NY], 2))*NX / 2.0f + (a[k + 1] * 0.5*(b[k] + b[k + NY]) - a[k - 1] * 0.5*(b[k - 1] + b[k + NY - 1]))*NY*KT / 2.0);
				u1 = u1 + (1.0 / Re)*((a[k + NY] - 2.0*a[k] + a[k - NY])*NX*NX + (a[k + 1] - 2.0*a[k] + a[k - 1])*NY*NY*KT*KT);
				u1 = a[k] + DT*(u1 - 1.0*NX*(c[k + NY] - c[k]));
				d[k] = u1;
			}
			else
				d[k] = a[k];

		}
	}
}
void VSh(double *a, double *b, double *c, double *d)
{

	double v1;
	int i, j, k;
	for (i = 0; i < NX; i++) {
		for (j = 0; j < NY; j++) {
			k = i *NY + j;

			if (i < NX - 1 && i>0 && j < NY - 2 && j>0 && ((i <= RX && j <= RY) != 1)){
				v1 = -((b[k + NY] * 0.5*(a[k] + a[k + 1]) - b[k - NY] * 0.5*(a[k - NY] + a[k - NY + 1]))*NX / 2.0 + (b[k + 1] * b[k + 1] - b[k - 1] * b[k - 1])*KT*NY / 2.0);
				v1 = v1 + (1.0f / Re)*((b[k + NY] - 2.0f*b[k] + b[k - NY])*NX*NX + (b[k + 1] - 2.0f*b[k] + b[k - 1])*NY*NY*KT*KT);
				v1 = b[k] + DT*(v1 - (1.0f*NY*KT)*(c[k + 1] - c[k]));
				d[k] = v1;
			}
			else
				d[k] = b[k];
		}
	}
}
void sumah(double *a, double *b)
{
	int i, j, k;
	for (i = 0; i < NX; i++) {
		for (j = 0; j < NY; j++) {
			k = i *NY + j;
			b[k] += a[k];
		}
	}
}
void rubnixh(double *a)
{
	int i, j, k;
	for (i = 0; i < NX; i++) {
		for (j = 0; j < NY; j++) {
			k = i *NY + j;

			if (i == 0 && j>RY)
				a[k] = 1.;
			if (i == NX - 2)
				a[k] = a[k - NY];
			if (j == 0)
				a[k] = 0.;
			if (j == NY - 1)
				a[k] = a[k - 1];
			if (j == RY && i < RX)
				a[k] = 0.;
			if (i == RX && j <= RY)
				a[k] = a[k + NY] * AI;
		}
	}
}
void rubniyh(double *a)
{
	int i, j, k;
	for (i = 0; i < NX; i++) {
		for (j = 0; j < NY; j++) {
			k = i *NY + j;

			if (j == 0)
				a[k] = AI*a[k + 1];
			if (j == NY - 2)
				a[k] = a[k - 1];
			if (i == 0)
				a[k] = 0.;
			if (i == NX - 1)
				a[k] = a[k - NY];
			if (j == RY && i <= RX)
				a[k] = AI*a[k + 1];
			if (i == RX && j < RY)
				a[k] = 0.;
		}
	}
}
void rubniph(double *a)
{
	int i, j, k;
	for (i = 0; i < NX; i++) {
		for (j = 0; j < NY; j++) {
			k = i *NY + j;


			if (j == 0 && i != 0 && i != NX - 1)
				a[k] = a[k + 1];
			if (j == NY - 1 && i != 0 && i != NX - 1)
				a[k] = 0.0;
			if (i == 0 && j != 0 && j != NY - 1)
				a[k] = 0;
			if (i == NX - 1 && j != 0 && j != NY - 1)
				a[k] = 0.;
			if (i == RX && j <= RY)
				a[k] = a[k + NY];
			if (i <= RX && j == RY)
				a[k] = a[k + 1];
		}
	}
}
void korekh(double *a, double *b, double *c)
{
	int i, j, k;
	for (i = 0; i < NX; i++) {
		for (j = 0; j < NY; j++) {
			k = i *NY + j;

			if (i>0 && i < NX - 2 && ((i <= RX && j <= RY) != 1)){
				a[k] -= DT*NX*(c[k + NY] - c[k]);
			}

			if (j>0 && j < NY - 2 && ((i <= RX && j <= RY) != 1)){
				b[k] -= DT*NY*KT*(c[k + 1] - c[k]);
			}

		}
	}
}

void poissh(double *a, double *b, double *c, double *d)
{
	int i, j, k;
	double corr;
	for (i = 0; i < NX; i++) {
		for (j = 0; j < NY; j++) {
			k = i *NY + j;
			corr = c[k];
			if (i > 0 && j > 0 && i < NX - 1 && j < NY - 1 && ((i <= RX && j <= RY) != 1)){
				c[k] = -(1.0 / (2 * (DT*NX*NX + DT*NY*NY*KT*KT)))*((-DT*NX*NX)*c[k + NY] + (-DT*NX*NX)*c[k - NY] + (-DT*NY*NY*KT*KT)*c[k + 1] + (-DT*NY*NY*KT*KT)*c[k - 1] + (1.0*NX)*(a[k] - a[k - NY]) + (1.0*NY*KT)*(b[k] - b[k - 1]));
			}
			corr -= c[k];
			if (corr < 0)
				corr *= (-1);
			d[k] = corr;
		}
	}
}



int main()
{

	clock_t t2, t1;
	FILE *fp, *fp1, *fp2;

	double *a_h, *a_h1, *b_h, *b_h1, *c_h, *c_h1, *corr_h;

	double  dx1, dy1;
	int j, i, t, k;

	fp = fopen("x.dat", "w");
	fp1 = fopen("y.dat", "w");
	fp2 = fopen("tlak.dat", "w");
	t1 = clock();

	a_h = (double*)malloc(sizeof(double)*NX*NY);
	a_h1 = (double*)malloc(sizeof(double)*NX*NY);
	b_h = (double*)malloc(sizeof(double)*NX*NY);
	b_h1 = (double*)malloc(sizeof(double)*NX*NY);
	c_h = (double*)malloc(sizeof(double)*NX*NY);
	c_h1 = (double*)malloc(sizeof(double)*NX*NY);
	corr_h = (double*)malloc(sizeof(double)*NX*NY);

	inixh(a_h);
	inixh(a_h1);
	inih(b_h);
	inih(b_h1);
	inih(c_h);
	inih(c_h1);


	printf("Start");
	for (t = 0; t<T; t++)
	{

		USh(a_h, b_h, c_h, a_h1);
		VSh(a_h, b_h, c_h, b_h1);
		rubnixh(a_h1);
		rubniyh(b_h1);

		for (i = 0; i<30; i++)
			poissh(a_h1, b_h1, c_h1, corr_h);

		rubniph(c_h1);
		korekh(a_h1, b_h1, c_h1);
		rubnixh(a_h1);
		rubniyh(b_h1);
		sumah(c_h1, c_h);
		inih(c_h1);
		rubniph(c_h);
		USh(a_h1, b_h, c_h, a_h);
		VSh(a_h, b_h1, c_h, b_h);
		rubnixh(a_h);
		rubniyh(b_h);

		for (i = 0; i<30; i++)
	    poissh(a_h, b_h, c_h1, corr_h);


		rubniph(c_h1);
		korekh(a_h, b_h, c_h1);
		rubnixh(a_h);
		rubniyh(b_h);
		sumah(c_h1, c_h);
		inih(c_h1);
		rubniph(c_h);


	}
	for (i = 0; i<NX; i++)
		for (j = 0; j<NY; j++) if (i != NX - 1)
			fprintf(fp, "\n%d %d %f", i, j, a_h[i*NY + j]);
	for (i = 0; i<NX; i++)
		for (j = 0; j<NY; j++) if (j != NY - 1)
			fprintf(fp1, "\n%d %d %f", i, j, b_h[i*NY + j]);
	for (i = 0; i<NX; i++)
		for (j = 0; j<NY; j++)
			fprintf(fp2, "\n%d %d %f", i, j, c_h[i*NY + j]);
	t2 = clock();
	printf(" %f\n", ((double)(t2 - t1)) / CLOCKS_PER_SEC);

	dx1 = 1.0*NX - 1;
	dy1 = KT*NY - 1;
	for (i = 0; i<NX; i++)
		dxx[i] = 1.0*i / dx1;
	for (i = 0; i<NY; i++)
		dyy[i] = 1.0*i / dy1;
	for (i = 0; i<NX; i++)
	{
		xi[i][0] = 0.0;
		yi[i][0] = 0.0;
		xi[i][NY - 1] = 0.0;
		yi[i][NY - 1] = 0.0;
	}
	for (i = 0; i<NY; i++)
	{
		xi[0][i] = 1.0;
		yi[0][i] = 0.0;
		xi[NX - 1][i] = 0.0;
		yi[NX - 1][i] = 0.0;
	}
	for (i = 1; i <= NX - 2; i++)
		for (j = 1; j <= NY - 1; j++)
		{
		k = (i - 1)*NY + j;

		xi[i][j] = 0.5*(a_h[k] + a_h[k + NY]);
		}
	for (i = 1; i <= NX - 1; i++)
		for (j = 1; j <= NY - 2; j++)
		{
		k = i*NY + j - 1;

		yi[i][j] = 0.5*(b_h[k] + b_h[k + 1]);
		}

	PRINT_VECT_VTK(xi, yi, dxx, dyy);
	PRINT_VTK(xi, yi, dxx, dyy);

	for (j = 0; j < NX; j++)
	for (i = 0; i < NY ; i++)
		 {
			dx[j][i] = (1.0*j / dx1);
			dy[j][i]=  (1.0*i / dy1);
		}



	PRINT_PLT();


	printf("END");


    return 0;
}