#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "header.h"
#include "mpi.h"

/*global variables:
 * X...........extent of the finest lattice in x-direction
 * Y...........extent of the finest lattice in y-direction
 * XL..........extent of the coarsest lattice in x-direction
 * YL..........extent of the coarsest lattice in y-direction
 * L...........size of the box, set to max(XL-1,YL-1)[later on]
 * NP..........# iterations used to solve the laplace equation or precision
 * NPHI........# of different potentials set at the start
 * NGROUND.....# of different grounded areas
 * outdir[]....output directory
 * method[]....method specified in the infile [jacobi / red-black / gauss-seidel]
 * exerc[].....exercise-label
 * *phi0.......pointer to the start-values of the potential
 * **ground....grounded areas.=>ground[#grounded areas][coordinate]
 *             example: ground[0][x1],ground[0][x2],ground[0][y1],ground[0][y2];
 *  */
int X, Y, XL, YL, L, F, NPHI, NGROUND, p=0;
char outdir[1000], method[100], exerc[100];
pot *phi0;
int **ground;
double NP;

/*
 * reads the input parameters and checks them for validity.
 *
 *  variables are explained in readme.in
 *
 * */
void read_infile(char *in) {
	FILE *file;
	int n;
	check_file(in, "r");
	file = fopen(in, "r");
	fscanf(file, "%s\n", exerc);
	fscanf(file, "%s\n", outdir);
	fscanf(file, "%s\n", method);
	if(!(strcmp(method, "gauss-seidel")==0 || strcmp(method, "jacobi")==0 || strcmp(method, "red-black")==0)){
		printf("%s is not a valid method.\n", method);
		abort();
	}/*abort if the method is unknown.*/

	fscanf(file, "X %d\n", &X);
	fscanf(file, "Y %d\n", &Y);
	fscanf(file, "%*s %d\n", &F);
	fscanf(file, "%*s %d\n", &XL);
	fscanf(file, "%*s %d\n", &YL);
	fscanf(file, "%*s %lf\n", &NP);
	fscanf(file, "%*s %d\n", &NPHI);
	phi0 = malloc(NPHI * sizeof(pot));
	for(n=0; n<NPHI; n++){
		fscanf(file, "%d %d %d %d %lf\n", &phi0[n].pos[0], &phi0[n].pos[1],&phi0[n].pos[2],&phi0[n].pos[3],&phi0[n].val);
		if(phi0[n].pos[0]<0 || phi0[n].pos[1]>XL || phi0[n].pos[2]<0 || phi0[n].pos[3]>YL){
			printf("X=0..%d Y=0..%d\tpotential outside boundary of the small grid\n", XL-1, YL-1);
			abort();
		}/*abort if the potential is outside of the boundaries of the coarsest lattice*/
	}
	fscanf(file, "%*s %d\n", &NGROUND);
	if(NGROUND>0 && (strcmp(method, "red-black")==0 || strcmp(method, "gauss-seidel")==0)){
		printf("Grounded areas are not implemented for the %s method.\n", method);
		abort();
	}
	ground=malloc(NGROUND*sizeof(int *));
	for(n=0; n<NGROUND; n++){
		ground[n]=malloc(4*sizeof(int));
		fscanf(file, "%d %d %d %d\n", &ground[n][0],&ground[n][1],&ground[n][2],&ground[n][3]);
		if(ground[n][0]<0 || ground[n][1]>=XL || ground[n][2]<0 || ground[n][3]>=YL){
			printf("X=0..%d Y=0..%d\tgrounded area outside or on the boundary\n", XL, YL);
			abort();
		}/*grounded areas must be inside of the box, they may not overlap with the boundary*/
	}
	fscanf(file, "plot %d", &p);
	fclose(file);

	/*prints the infile in the terminal for checks.*/
	printf("%s\n", exerc);
	printf("%s\n", outdir);
	printf("%s\n", method);
	printf("X %d\n", X);
	printf("Y %d\n", Y);
	printf("F %d\n", F);
	printf("XL %d\n", XL);
	printf("YL %d\n", YL);
	printf("NP %.2e\n", NP);
	printf("npotential %d\n", NPHI);
	for(n=0; n<NPHI; n++){
		printf("phi[%d]=% .2e at %d %d %d %d\n", n, phi0[n].val, phi0[n].pos[0], phi0[n].pos[1], phi0[n].pos[2],phi0[n].pos[3]);
	}
	printf("nground %d\n", NGROUND);
	for(n=0; n<NGROUND; n++){
		printf("ground[%d] at %d %d %d %d\n", n, ground[n][0],ground[n][1],ground[n][2],ground[n][3]);
	}
}



int main(int arg, char *args[]) {
	int i, x, y, Xstep, Ystep, ix, iy, step;
	double *phi, *phibuf, *range, x1, x2, y1, y2, buf;
	char out[1000], title[1000], xlabel[100], ylabel[100], plot[1000], dif[1000];
	static int my_rank;
	FILE *file;
	plotinfo *pi;

	if(arg!=2){
		printf("incorrect number of input parameters.\n");
		abort();
	}

	read_infile(args[1]);/*get input data and test it*/

	/*initialization of the MPI-process, required to measure processor-time*/
	MPI_Init(&arg, &args);

	phi=malloc(XL*YL*sizeof(double));
	pi=malloc(3*sizeof(plotinfo));
	range=malloc(6*sizeof(double));
	for(i=0; i<6; i++)
		range[i]=0.;

	/*initialize potential*/
	for(x=0; x<XL; x++){
		for(y=0; y<YL; y++){
			phi[y*XL+x]=0.;
		}
	}
	for(i=0; i<NPHI; i++){
		printf("nphi[%d]\n", i);
		for(x=phi0[i].pos[0]; x<=phi0[i].pos[1]; x++){
			for(y=phi0[i].pos[2]; y<=phi0[i].pos[3]; y++){
				phi[y*XL+x]=phi0[i].val;
				printf("%d %d %d %f\n", i, x, y, phi[y*XL+x]);
			}
		}
	}

	Xstep=XL;/*lattice extent of the current step*/
	Ystep=YL;
	step=1;
	while(Xstep<=X && Ystep<=Y){
		printf("solving lattice: %dx%d\n", Xstep, Ystep);
		sprintf(dif, "%s/dif_X%d_Y%d_%s_%s_m%d.dat", outdir, Xstep, Ystep,exerc, method, step);
		sprintf(out, "%s/%s_%s.log", outdir, exerc, method);
		/*compute the solution of the laplace equation using the functions defined in functions.c*/
		if(0==strcmp(method, "jacobi"))
			jacobi(Xstep, Ystep, NP, phi, NPHI, phi0, NGROUND, ground, dif, out);
		else if(0==strcmp(method, "red-black"))
			red_black(Xstep, Ystep, NP, phi, NPHI, phi0, dif, out);
		else if(0==strcmp(method, "gauss-seidel"))
			gauss_seidel(Xstep, Ystep, NP, phi, NPHI, phi0, dif, out);

		/*define the plotranges*/
		if(Xstep>Ystep)
			L=Xstep-1;
		else
			L=Ystep-1;
		range[0]=0, range[1]=(double)(Xstep-1)/L,range[2]=0, range[3]=(double)(Ystep-1)/L;/*plotranges*/
		/*generate the output for the potential*/
		sprintf(out, "%s/phi_X%d_Y%d_%s_%s.dat", outdir, Xstep, Ystep,exerc, method);
		printf("%s\n", out);
		check_file(out, "w");
		file=fopen(out, "w");
		fprintf(file, "#Plotdata for:\n");
		fprintf(file, "#%s\n", exerc);
		fprintf(file, "#%s\n", outdir);
		fprintf(file, "#%s\n", method);
		fprintf(file, "#X %d\n", X);
		fprintf(file, "#Y %d\n", Y);
		fprintf(file, "#NP %.2e\n", NP);
		fprintf(file, "#npotential %d\n", NPHI);
		for(i=0; i<NPHI; i++){
			fprintf(file, "#phi[%d]=% .2e at %d %d %d %d\n", i, phi0[i].val,
					phi0[i].pos[0], phi0[i].pos[1], phi0[i].pos[2],phi0[i].pos[3]);
		}
		fprintf(file, "#nground %d\n", NGROUND);
		for(i=0; i<NGROUND; i++){
			fprintf(file, "#ground[%d] at %d %d %d %d\n", i,
					ground[i][0],ground[i][1],ground[i][2],ground[i][3]);
		}
		for(x=0; x<Xstep; x++){
			for(y=0; y<Ystep; y++){
				fprintf(file, "%4f %4f % .6e\n",(double)x/L, (double)y/L, phi[y * Xstep + x]);
			}
			fprintf(file, "\n");/*to create "set pm3d map"
			each matrix line has to be separated by the next by an additional line.*/
		}
		fclose(file);
		if(p==1){
			sprintf(pi[0].color, "red");
			sprintf(pi[0].data, "phi_X%d_Y%d_%s_%s.dat", Xstep, Ystep,exerc, method);
			sprintf(pi[0].label, "$\\Phi(NP=%.2e)$", NP);
			pi[0].lw=1;
			pi[0].pt=1;
			pi[0].x=1;
			pi[0].y=2;
			sprintf(xlabel, "x[a]");
			sprintf(ylabel, "y[a]");
			sprintf(title, "$\\Phi(x,y)$ at NP=%.2e and $m=%d$", NP, step);
			sprintf(plot, "potential_X%d_Y%d_%s_%s", Xstep, Ystep,exerc, method);
			plot_data_color(outdir, plot, title, xlabel, ylabel, range, Xstep, Ystep, 1, pi);
		}
		/*generate the output for the electric field*/
		sprintf(out, "%s/E_X%d_Y%d_%s_%s.dat", outdir, Xstep, Ystep,exerc, method);
		printf("%s\n", out);
		check_file(out, "w");
		file=fopen(out, "w");
		for(x=0; x<Xstep-1; x++){
			for(y=0; y<Ystep-1; y++){
				/*discretized gradient of phi:*/
				x1=(phi[y*Xstep+x]+phi[(y+1)*Xstep+x])/2;
				x2=(phi[y*Xstep+x+1]+phi[(y+1)*Xstep+x+1])/2;
				y1=(phi[y*Xstep+x]+phi[y*Xstep+x+1])/2;
				y2=(phi[(y+1)*Xstep+x]+phi[(y+1)*Xstep+x+1])/2;
				buf=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))/2;
				fprintf(file, "%f %f %f %f\n", (double)(2*x+1)/2,(double)(2*y+1)/2,	-(x2-x1)/2/buf, -(y2-y1)/2/buf);
			}
			fprintf(file, "\n");
		}
		fclose(file);
		if(p==1){
			sprintf(pi[0].color, "black");
			sprintf(pi[0].data, "E_X%d_Y%d_%s_%s.dat", Xstep, Ystep,exerc, method);
			sprintf(pi[0].label, "$\\vec{E}$");
			pi[0].lw=3;
			pi[0].x=1;
			pi[0].y=2;
			pi[0].u=3;
			pi[0].v=4;
			sprintf(xlabel, "x[a]");
			sprintf(ylabel, "y[a]");
			sprintf(title, "$\\vec{E}(x,y)=-\\nabla\\Phi$ at $m=%d$", step);
			sprintf(plot, "E_X%d_Y%d_%s_%s", Xstep, Ystep,exerc, method);
			range[0]=0., range[1]=(double)Xstep, range[2]=0., range[3]=(double)Ystep;
			plot_data_vecfield(outdir, plot, title, xlabel, ylabel, range, Xstep, Ystep, 1, pi);
		}
		/*generate the output for the magnetic field*/
		sprintf(out, "%s/B_X%d_Y%d_%s_%s.dat", outdir, Xstep, Ystep,exerc, method);
		printf("%s\n", out);
		check_file(out, "w");
		file=fopen(out, "w");
		for(x=0; x<Xstep; x++){
			for(y=0; y<Ystep; y++){
				if((x>0 && x<Xstep-1) && (y>0 && y<Ystep-1)){
					/*discretized \nabla x \phi:*/
					x1=(phi[(y+1)*Xstep+x]-phi[(y-1)*Xstep+x])/2.;
					y1=-(phi[y*Xstep+x+1]-phi[y*Xstep+x-1])/2.;
					buf=sqrt((x1*x1+y1*y1));
					fprintf(file, "%f %f %f %f\n", (double)x,(double)y,	x1/buf, y1/buf);
				}
				else{
					fprintf(file, "%f %f %f %f\n", (double)x,(double)y,	0., 0.);
				}
			}
			fprintf(file, "\n");
		}
		fclose(file);
		if(p==1){
			sprintf(pi[0].color, "black");
			sprintf(pi[0].data, "B_X%d_Y%d_%s_%s.dat", Xstep, Ystep,exerc, method);
			sprintf(pi[0].label, "$\\vec{E}$");
			pi[0].lw=3;
			pi[0].x=1;
			pi[0].y=2;
			pi[0].u=3;
			pi[0].v=4;
			sprintf(xlabel, "x[a]");
			sprintf(ylabel, "y[a]");
			sprintf(title, "$\\vec{B}(x,y)=\\nabla\\times A$ at $m=%d$", step);
			sprintf(plot, "B_X%d_Y%d_%s_%s", Xstep, Ystep,exerc, method);
			range[0]=0., range[1]=(double)Xstep, range[2]=0., range[3]=(double)Ystep;
			plot_data_vecfield(outdir, plot, title, xlabel, ylabel, range, Xstep, Ystep, 1, pi);
			sprintf(out, "rm %s/*.dat", outdir);
			printf("%s\n", out);
			system(out);
		}
/*allocate memory for a temporary array and store the current potential*/
		phibuf=malloc(Xstep*Ystep*sizeof(double));
		for(x=0; x<Xstep; x++){
			for(y=0; y<Ystep; y++){
				phibuf[y*Xstep+x]=phi[y*Xstep+x];
			}
		}
		free(phi);/*free the memory used by the current step's potential*/
		Xstep*=F;/*increase the lattice-extent to the next level*/
		Ystep*=F;
		/*allocate memory for the next step's potential and store the old
		 * potential, such that one lattice site of the previous step is stored in
		 * F*F lattice sites at the appropriate coordinates.
		 * Example: F=2: x-> xx
		 * 					 xx*/
		phi=malloc(Xstep*Ystep*sizeof(double));
		for(x=0, ix=0; x<Xstep; x++){
			for(y=0, iy=0; y<Ystep; y++){
				phi[y*Xstep+x]=phibuf[iy*(int)Xstep/F+ix];
				if((y+1)%F==0)
					iy++;
			}
			if((x+1)%F==0)
				ix++;
		}
		free(phibuf);/*free memory used by the temporary potential*/
		/*
		 * as the lattice size is increased the predefined potentials and grounded
		 * areas must grow by the same factor:
		 * */
		for(i=0; i<NPHI; i++){
			for(x=0; x<4; x++){
				phi0[i].pos[x]*=F;
			}
		}
		for(i=0; i<NGROUND; i++){
			for(x=0; x<4; x++){
				ground[i][x]*=F;
			}
		}
		step++;
	}/*end loop over lattice-volumes*/

	return 0;
}