#ifndef lbmInOut
#define lbmInOut

#include <math.h>
#include "nw2cuda.h"

#undef __FUNC__
#define __FUNC__ "get_int"

int get_int(char int_variable_name[], int *int_variable)
{
	char dummyget[128];
	sprintf(dummyget,"-%s",int_variable_name);
	GetInt(dummyget,int_variable);

	return 0;
}

#undef __FUNC__
#define __FUNC__ "get_Real"

int get_Real(char scalar_variable_name[], double *scalar_variable)
{
	char dummyget[128];
	sprintf(dummyget,"-%s",scalar_variable_name);
	GetDouble(dummyget,scalar_variable);
	return 0;
}

#undef __FUNC__
#define __FUNC__ "lb_input"

int lb_input(int *hnodes, int *hnodes_x, int *hnodes_y, int *hniter, int *hniter_cycle, int *hncycles,int *hlast_iter,
		int *hlnodes, int *hlnodes_x, int *hlnodes_y,
		double *hbbspeed, double *htcritical, double *htfluid, double *hchic2, double *htwochic2, double *htwochi2c4, double *hdtimetemp,
		double *hRTzero, double *htemperature, double *hdelta_time, double *hkappa, double *hrho_zero,
		double *htau, double *hdelta_space, double *htwo_delta_space,double *hdelta_space_three_three,
		const int hnw, const double hlength, const double hone, const double htwo, const double hthree,
		const double hchi)
{
	char snodes_x[]="hnodes_x", snodes_y[]="hnodes_y";
	char stemperature[]="htemperature", sRTzero[]="hRTzero";
	char sdelta_time[]="hdelta_time", skappa[]="hkappa";
	char sniter[]="hniter", sniter_cycle[]="hniter_cycle",
			srho_zero[]="hrho_zero";
	char stau[]="htau", sncycles[]="hncycles";
	char slast_iter[]="hlast_iter";

	get_int(snodes_x, hnodes_x);
	get_int(snodes_y, hnodes_y);
	get_int(sniter, hniter);
	get_int(sniter_cycle, hniter_cycle);
	get_int(sncycles, hncycles);
	get_int(slast_iter, hlast_iter);

	get_Real(sRTzero, hRTzero);
	get_Real(stemperature, htemperature);
	get_Real(sdelta_time, hdelta_time);
	get_Real(skappa, hkappa);
	get_Real(srho_zero, hrho_zero);
	get_Real(stau, htau);

	//*hdelta_space = hlength / ((double)(*hnodes_x)); not working for nodes_x > 128
	*hdelta_space = hlength / ((double) 128);
	*htwo_delta_space = htwo * (*hdelta_space);
	*hdelta_space_three_three = ((*hdelta_space) * (*hdelta_space) * (*hdelta_space)) / hthree;

	*hnodes = (*hnodes_x) * (*hnodes_y);
	*hlnodes_x = (*hnodes_x) + (2 * hnw);
	*hlnodes_y = (*hnodes_y) + (2 * hnw);
	*hlnodes = (*hlnodes_x) * (*hlnodes_y);
	//hnw = 2;

	*htcritical = hone / (*hRTzero);
	*htfluid = (*htemperature) / (*hRTzero);

	*hchic2 = *htfluid;
	*hbbspeed = sqrt((*htfluid) / hchi);

	*htwochic2 = htwo * (*hchic2);
	*htwochi2c4 = htwo * (*hchic2) * (*hchic2);

	*hdtimetemp  = (*hdelta_time) / (*hchic2);

	return SUCCESS;
}

int build_names(char *name,
		const int hnodes_x, const int hnodes_y, const double hdelta_space,
		const double hdelta_time, const double htau, const double hrho_zero,
		const double hRTzero, const double htfluid, const double hkappa)
{
	sprintf(name,
			"_X%3.3d_Y%3.3d_S%6.4f_T%7.5f_Ta%5.3e_R%5.3f_Z%4.2f_T%4.2f_K%8.6f",
			hnodes_x,hnodes_y,hdelta_space,hdelta_time,htau,hrho_zero,hRTzero,htfluid*hRTzero,hkappa);

	return SUCCESS;
}

#undef __FUNC__
#define __FUNC__ "profile"

int profile(const char *name, const double *hgf,
		const int hnodes, const int hnodes_x, const int hnodes_y,const int hkn,
		const int hniter, const double hlength,
		const double hzero, const double htwo, const double hdelta_space)
{
	int ix, iy, valxv, ixv;
	double ntot, maxval;
	double *x_position;
	double *nprofile;

	char filename[128];
	FILE *xhandle;

	nprofile = (double*)calloc(hnodes_x,sizeof(double));
	CHK_NULL_RETURN(nprofile);

	x_position = (double*)calloc(hnodes_x,sizeof(double));
	CHK_NULL_RETURN(x_position);

	ntot = hzero;

	for(iy=0; iy < hnodes_y; iy++)
	{
		for(ix=0; ix < hnodes_x; ix++)
		{
			ntot += hgf[hkn * hnodes + iy * hnodes_x + ix];
			//aznloc[ix][iy] =  hgf[hkn*hnodes + iy * hnodes_x + ix];
		}
	}

	for(ix=0; ix < hnodes_x; ix++)
	{
		x_position[ix] = hdelta_space * ((double) (ix+1)) - hdelta_space / htwo;
		nprofile[ix] = hgf[hkn * hnodes + (hnodes_y / 2) * hnodes_x + ix];

	}

	printf("iter=%10.10d  ntot=%25.20e\n",hniter,ntot);

	sprintf(filename,"BP%s.%09d",name, hniter);
	xhandle = fopen(filename,"w"); 
	CHK_NULL_RETURN(xhandle);

	for(ix=0; ix< hnodes_x; ix++)
	{
		fprintf(xhandle,"%15.6e %15.6e \n",
				x_position[ix]-hlength/htwo, nprofile[ix]);
	}

	fclose(xhandle);

	sprintf(filename,"BV%s.%09d",name, hniter);
	xhandle = fopen(filename,"w"); 
	CHK_NULL_RETURN(xhandle);

	fprintf(xhandle, "P2\n%4d%5d\n255\n", hnodes_x, hnodes_y);

	maxval = 2.500;
	ixv = 0;

	for(ix=0; ix < hnodes_x; ix++)
	{
		for(iy=0; iy < hnodes_y; iy++)
		{
			ixv++;
			valxv = (int) floor ((hgf[hkn * hnodes + iy * hnodes_x + ix] * 255.000) / maxval);
			valxv = 256 - valxv;

			if(valxv > 255)
			{
				valxv = 255;
			}
			if(valxv < 0)
			{
				valxv = 0;
			}
			fprintf(xhandle, "%4d", valxv);

			if( ixv == 15)
			{
				ixv = 0;
				fprintf(xhandle,"\n");
			}
		}
	}
	if(ixv)
	{
		fprintf(xhandle,"\n");
	}
	fclose(xhandle);

	/* write vtk files */
	sprintf(filename,"PV%s.%09d.vtk",name,hniter);
	xhandle = fopen(filename,"w");
	CHK_NULL_RETURN(xhandle);

	fprintf(xhandle,"# vtk DataFile Version 3.0\n");
	fprintf(xhandle,"Output for PARAVIEW\n");
	fprintf(xhandle,"ASCII \nDATASET STRUCTURED_GRID\n");
	fprintf(xhandle,"DIMENSIONS %d %d %d\n", hnodes_x,hnodes_y,1);
	fprintf(xhandle,"POINTS %d float\n",(hnodes_x)*(hnodes_y));

	for(ix=0; ix < hnodes_x; ix++)
	{
		for(iy=0; iy < hnodes_y; iy++)
		{
			fprintf(xhandle,"%.4e %.4e %.4e\n",
					hdelta_space * (double) ix, hdelta_space * (double) iy, 0.0);
		}
	}
	fprintf(xhandle,"POINT_DATA %d\n", (hnodes_x)*(hnodes_y));
	fprintf(xhandle,"SCALARS rho float 1 \nLOOKUP_TABLE default\n");
	for(ix=0; ix < hnodes_x; ix++)
	{
		for(iy=0; iy < hnodes_y; iy++)
		{
			fprintf(xhandle,"%.4e\n",hgf[hkn * hnodes + iy * hnodes_x + ix]);
		}
	}
	fclose(xhandle);
	/* end write vtk files */

	free(x_position);
	free(nprofile);

	return SUCCESS;
}


#endif