#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdbool.h>

#include <CImg.h>

#include <omp.h>

#define DELTA_X 1.0
#define DELTA_Y 1.0

#define WIDTH 100
#define HEIGHT 100
#define ITERATIONS 2000
#define EPSILON 0.00

#define OMEGA 1.0


float grid[2][WIDTH][HEIGHT];
bool fixed[WIDTH][HEIGHT];


int main(int argc, char ** argv)
{
	float beta2 = pow(DELTA_X / DELTA_Y, 2);
	//printf("Beta2 is %f\n", beta2);
	omp_set_nested(1);
	
	for (int i = 0; i < WIDTH; ++i)
	{
		for (int j = 0; j < HEIGHT; ++j)
		{
			grid[0][i][j] = 0;
			fixed[i][j] = (i == 0 || j == 0 || i == WIDTH-1 || j == HEIGHT-1);
		}
	}

	// Add the electrodes
	

	int n = 1; int p = 0;
	int k = 0;
	float maxChange = 0.0;
	for (k = 0; k < ITERATIONS; ++k)
	{
		maxChange = 0.0;
		#pragma omp parallel for
		for (int i = 0; i < WIDTH; ++i)
		{
			#pragma omp parallel for
			for (int j = 0; j < WIDTH; ++j)
			{
				if (fixed[i][j] == true)
				{
					grid[n][i][j] = grid[p][i][j];
					continue;
				}
				
				grid[n][i][j] = (1.0 - OMEGA) * grid[p][i][j];
			

				grid[n][i][j] += OMEGA * (grid[n][i-1][j] + grid[p][i+1][j] + beta2 * (grid[n][i][j-1] + grid[p][i][j+1]));
				
				grid[n][i][j] = grid[n][i][j] / (2.0 * (1.0 + beta2));

				float change = (grid[n][i][j] > grid[p][i][j]) ? grid[p][i][j] - grid[n][i][j] : grid[p][i][j] - grid[n][i][j];
				if (grid[p][i][j] != 0.00)
					change = change / grid[p][i][j];
				else
					change = 1.0;
				if (change >= maxChange)
					maxChange = change;
				
				

			}
		}
		
		//fprintf(stdout, "# Complete iteration %d; max change was %f\n", k, maxChange);

		
	

		char file[100];
		sprintf(file, "%d.dat", k);
		FILE * out = fopen(file, "w");

		fprintf(out, "# Grid dimensions %d %d\n", WIDTH, HEIGHT);
		fprintf(out, "# dx = %f, dy = %f, beta2 = %f\n", DELTA_X, DELTA_Y, beta2);
		fprintf(out, "# Ran for %d iterations\n", k);
		fprintf(out, "# Largest change is %f\n", maxChange);
		for (int i = 0; i < WIDTH; ++i)
		{
			for (int j = 0; j < HEIGHT; ++j)
			{
				fprintf(out, "%f\t%f\t%f\n", i*DELTA_X, j*DELTA_Y, grid[n][i][j]);
			}
			fprintf(out, "\n");
		}
		fprintf(out, "# DONE\n");
		if (out != stdout && out != stderr) fclose(out);

		if (maxChange < EPSILON)
			break;

		n = (n == 1) ? 0 : 1;
		p = (p == 1) ? 0 : 1;
	}

}

//EOF