#ifndef lbmInit
#define lbmInit

extern  double randomUniformInterval(double a, double b);
//#include "nw2cuda.h"

#undef __FUNC__
#define __FUNC__ "init_base_vectors"


int init_base_vectors(double *hbweight,double *hcweight, double *hbspeed,double *hcspeed_x, double *hcspeed_y,
		double *hcfl, const double hbbspeed, const double hzero, const double hdelta_time, const double hdelta_space)
{

	hbweight[0] = ((double) 2) / ((double) 3);
	hbweight[1] = ((double) 1) / ((double) 6);

	hbspeed[0] = hzero;
	hbspeed[1] = hbbspeed;

	hcweight[0]  =   hbweight[0] * hbweight[0];
	hcspeed_x[0] =   hbspeed[0];
	hcspeed_y[0] =   hbspeed[0];

	hcweight[1]  =   hbweight[1] * hbweight[0];
	hcspeed_x[1] =   hbspeed[1];
	hcspeed_y[1] =   hbspeed[0];

	hcweight[2]  =   hbweight[0] * hbweight[1];
	hcspeed_x[2] =   hbspeed[0];
	hcspeed_y[2] =   hbspeed[1];

	hcweight[3]  =   hbweight[1] * hbweight[0];
	hcspeed_x[3] = - hbspeed[1];
	hcspeed_y[3] =   hbspeed[0];

	hcweight[4]  =   hbweight[0] * hbweight[1];
	hcspeed_x[4] =   hbspeed[0];
	hcspeed_y[4] = - hbspeed[1];

	hcweight[5]  =   hbweight[1] * hbweight[1];
	hcspeed_x[5] =   hbspeed[1];
	hcspeed_y[5] =   hbspeed[1];

	hcweight[6]  =   hbweight[1] * hbweight[1];
	hcspeed_x[6] = - hbspeed[1];
	hcspeed_y[6] =   hbspeed[1];

	hcweight[7]  =   hbweight[1] * hbweight[1];
	hcspeed_x[7] = - hbspeed[1];
	hcspeed_y[7] = - hbspeed[1];

	hcweight[8]  =   hbweight[1] * hbweight[1];
	hcspeed_x[8] =   hbspeed[1];
	hcspeed_y[8] = - hbspeed[1];


	*hcfl = hbspeed[1] * hdelta_time / hdelta_space;

	return 0;
}

#undef __FUNC__
#define __FUNC__ "init_arrays"

int init_arrays(double *hgf,double *hfeq,
		const double *hcspeed_x, const double *hcspeed_y, const double *hcweight,
		const int hnodes, const int hnodes_x, const int hnodes_y, const int hnvelocities,
		const int hkn, const double hzero, const double hone,
		const double hrho_zero,const double hchic2,
		const double htwochic2, const double htwochi2c4)
{
	int ix, iy, k, key_init = 0;
	double interval_start, interval_end;
	double ux, uy, hedotu;

#if 1
	switch(key_init)
	{
	case 0:         /*  homogeneous system  */
	default:

		interval_start = 0.999 * hrho_zero;
		interval_end = 1.001*hrho_zero;

		for (iy = 0; iy < hnodes_y; iy++)
		{
			for (ix=0; ix < hnodes_x; ix++)
			{
				hgf[hkn * hnodes + iy * hnodes_x + ix] = randomUniformInterval(interval_start, interval_end);
				//hgf[hkn * hnodes + iy * hnodes_x + ix] = hrho_zero;
			}
		}
		break;
	case 1:         /*  other  */
		break;
	}

	for(iy=0; iy < hnodes_y; iy++)
	{
		for(ix=0; ix < hnodes_x; ix++)
		{
			ux = uy = hzero;

			hfeq[0] = hone - (ux*ux + uy*uy) / htwochic2;

			for(k=1; k<hnvelocities; k++)
			{
				hedotu  = hcspeed_x[k] * ux + hcspeed_y[k] * uy;

				hfeq[k] = hfeq[0] + hedotu / hchic2 + hedotu * hedotu / htwochi2c4;

				hgf[k * hnodes + iy * hnodes_x + ix] = hfeq[k] * hgf[hkn * hnodes + iy * hnodes_x + ix] * hcweight[k];

			}
			hgf[0 * hnodes + iy * hnodes_x + ix] = hfeq[0] * hgf[hkn * hnodes + iy * hnodes_x + ix] * hcweight[0];
		}
	}
#endif
	return 0;
}


#endif