#ifndef _INCLUDED_L3_DIFF_H
#define _INCLUDED_L3_DIFF_H

// *Finite differences methods* Copyright (C) Krzysztof Bosak, 1999-12-02...1999-12-03
// All rights reserved.
// kbosak@box43.pl
// http://www.kbosak.prv.pl

#include "l3_array.h"

template<class real>
class Lattice
{
	int _steps;
	int _stepspow2;
	real _lostvalue;
	basearray<int> _indexy;
	basearray<real> _lattice;
public:
	Lattice(int steps, bool clear)
			: _steps(steps),
			_stepspow2(steps*steps),
			_indexy(steps),
			_lattice(steps*steps)
	{
		for(int i=0, iy=0; i<steps; i++, iy+=steps)
		{
			_indexy[i]=iy;
		}
		if(clear==true)
		{
			Clear();
		}
	}
	inline int Size() const
	{
		return _steps;
	}
	inline real& operator()(int x, int y)
	{// SOR needs wide array boundaries...
		assert(x>=-2 && x<_steps+2);
		assert(y>=-2 && y<_steps+2);
		if(x<0 || x>=_steps)
		{
			return _lostvalue=0;
		}
		if(y<0 || y>=_steps)
		{
			return _lostvalue=0;
		}
		return _lattice[x+_indexy[y]];
	}
	inline real operator()(int x, int y) const
	{// SOR needs wide array boundaries...
		assert(x>=-2 && x<_steps+2);
		assert(y>=-2 && y<_steps+2);
		if(x<0 || x>=_steps)
		{
			return 0;
		}
		if(y<0 || y>=_steps)
		{
			return 0;
		}
		return _lattice[x+_indexy[y]];
	}
	inline real& Value(int x, int y)
	{ // Get value by reference
		assert(x>=0 && x<_steps);
		assert(y>=0 && y<_steps);
		return _lattice[x+_indexy[y]];
	}
	inline real Value(int x, int y) const
	{ // Get value
		assert(x>=0 && x<_steps);
		assert(y>=0 && y<_steps);
		return _lattice[x+_indexy[y]];
	}
	inline real& operator[](int i)
	{
		return _lattice[i];
	}
	inline real operator[](int i) const
	{
		return _lattice[i];
	}
	inline void Clear()
	{
		for(int j=0; j<_stepspow2; j++)
		{
			_lattice[j]=0;
		}
	}
	inline void Swap(Lattice& lat2)
	{
		assert(_steps==lat2._steps);
		_lattice.swap(lat2._lattice);
	}
	inline void GetData(Lattice& lat)
	{
		_lattice=lat._lattice;
	}
	void Deriv1x(Lattice& lat, real h) const
	{
		const real c1=1/(12*h);
		const real c2=8/(12*h);
		const int sm1=_steps-1;
		const int sm2=_steps-2;
		int y, iy;
		// Interior
		for(y=0, iy=0; y<_steps; y++, iy+=_steps)
		{
			for(int x=2; x<sm2; x++)
			{
				const int xiy=x+iy;
				const real *ptr=&_lattice[xiy];
				lat[xiy]=c1*(ptr[-2]-ptr[2])+c2*(ptr[1]-ptr[-1]);// O(4) centered
			}
		}
		const real c3=3/(2*h);
		const real c4=4/(2*h);
		const real c5=1/(2*h);
		// Borders
		for(y=0, iy=0; y<_steps; y++, iy+=_steps)
		{
			int xiy=1+iy;
			const real *ptr=&_lattice[xiy];
			lat[xiy]=c5*(ptr[1]-ptr[-1]);// O(2) centered

			ptr=&_lattice[iy];
			lat[iy]=-c3*ptr[0]+c4*ptr[1]-c5*ptr[2];// O(2) decentered

			xiy=sm2+iy;
			ptr=&_lattice[xiy];
			lat[xiy]=c5*(ptr[1]-ptr[-1]);// O(2) centered

			xiy=sm1+iy;
			ptr=&_lattice[xiy];
			lat[xiy]=c3*ptr[0]-c4*ptr[-1]+c5*ptr[-2];// O(2) decentered
		}
	}
	void Deriv1y(Lattice& lat, real h) const
	{
		const real c1=1/(12*h);
		const real c2=8/(12*h);
		const int sm1=_steps-1;
		const int sm2=_steps-2;
		const int steps2=_steps<<1;
		int y, iy;

		// Interior
		for(y=2, iy=steps2; y<sm2; y++, iy+=_steps)
		{
			for(int x=0; x<_steps; x++)
			{
				const int xiy=x+iy;
				const real *ptr=&_lattice[xiy];
				lat[xiy]=c1*(ptr[-steps2]-ptr[steps2])+c2*(ptr[_steps]-ptr[-_steps]);// O(4) centered
			}
		}

		const real c3=3/(2*h);
		const real c4=4/(2*h);
		const real c5=1/(2*h);
		// Borders
		for(int x=0; x<_steps; x++)
		{
			int xiy=x+_steps;
			const real *ptr=&_lattice[xiy];
			lat[xiy]=c5*(ptr[_steps]-ptr[-_steps]);// O(2) centered

			ptr=&_lattice[x];
			lat[x]=-c3*ptr[0]+c4*ptr[_steps]-c5*ptr[steps2];// O(2) decentered

			xiy=x+sm2*_steps;
			ptr=&_lattice[xiy];
			lat[xiy]=c5*(ptr[_steps]-ptr[-_steps]);// O(2) centered

			xiy=x+sm1*_steps;
			ptr=&_lattice[xiy];
			lat[xiy]=c3*ptr[0]-c4*ptr[-_steps]+c5*ptr[-steps2];// O(2) decentered
		}
	}
	void LaplacianNC(Lattice& lat, const real h) const
	{
		const real c1=1/(h*h);
		const real c2=-4/(h*h);
		const int sm1=_steps-1;
		int y, iy;

		// Interior
		for(y=1, iy=_steps; y<sm1; y++, iy+=_steps)
		{
			for(int x=1; x<sm1; x++)
			{
				const int xiy=x+iy;
				const real *ptr=&_lattice[xiy];
				lat[xiy]=c2*ptr[0]+c1*(ptr[-1]+ptr[1]+ptr[-_steps]+ptr[_steps]);// O(2) centered
			}
		}

		// Corners
		lat[0]=c2*_lattice[0]+c1*(_lattice[1]+_lattice[_steps]);// upper left
		lat[sm1]=c2*_lattice[sm1]+c1*(_lattice[sm1-1]+_lattice[sm1+_steps]);// upper right
		const int ylow=sm1*_steps;
		lat[ylow]=c2*_lattice[ylow]+c1*(_lattice[ylow+1]+_lattice[ylow-_steps]);// lower left
		const int ylowsm1=ylow+sm1;
		lat[ylowsm1]=c2*_lattice[ylowsm1]+c1*(_lattice[ylowsm1-1]+_lattice[ylowsm1-_steps]);// lower right

		// Upper and lower borders
		for(int x=1; x<sm1; x++)
		{
			const real *ptr=&_lattice[x];// y=0
			lat[x]=c2*ptr[0]+c1*(ptr[-1]+ptr[1]+ptr[_steps]);// O(2) centered

			const int xiy=x+ylow;
			ptr=&_lattice[xiy];// y=sm1
			lat[xiy]=c2*ptr[0]+c1*(ptr[-1]+ptr[1]+ptr[-_steps]);// O(2) centered
		}

		// Left and right borders
		for(y=1, iy=_steps; y<sm1; y++, iy+=_steps)
		{
			const real *ptr=&_lattice[iy];// x=0
			lat[iy]=c2*ptr[0]+c1*(ptr[1]+ptr[-_steps]+ptr[_steps]);// O(2) centered

			const int xiy=sm1+iy;
			ptr=&_lattice[xiy];// x=sm1
			lat[xiy]=c2*ptr[0]+c1*(ptr[-1]+ptr[-_steps]+ptr[_steps]);// O(2) centered
		}
	}
};

#endif //_INCLUDED_L3_DIFF_H