using System;
using System.Collections.Generic;

namespace Science.Mathematics.PartialDifferentialEquation
{
	/// <summary>
	/// PoissonEquationTwoDimension
	/// </summary>
	public class PoissonEquation2D
	{
		public delegate double Source(double x, double y);

		public PoissonEquation2D()
		{
		}
	
		private double delta, rj;
		public double MeshSize
		{
			get{return delta;}
			set{delta=value;}
		}
		public double JacobiRadius
		{
			get{return rj;}
			set{rj=value;}
		}

		private Domain2D dtd;
		public Domain2D Domain
		{
			set{dtd=value;}
		}

		private Source source;
		public Source SourceTerm
		{
			set{source=value;}
		}

        private List<double> al = new List<double>();
		public void SetBoundaryValue(double xFrom, double yFrom, 
			double xTo, double yTo, double boundaryValue)
		{
			al.Add(xFrom);
			al.Add(xTo);
			al.Add(yFrom);
			al.Add(yTo);
			al.Add(boundaryValue);
		}


		private bool[,] isItInDomain;
		private double[,] bv;
		public void Solve()
		{
			double x,y;
			int nbv = al.Count/5;
			double xf = dtd.LowerBoundOfX;
			double xt = dtd.UpperBoundOfX;
			double yf = dtd.LowerBoundOfY;
			double yt = dtd.UpperBoundOfY;
			int jx = (int) ((xt - xf) / delta) + 1;
			int jy = (int) ((yt - yf) / delta) + 1;
			double[,] a = new double[jx,jy];
			double[,] b = new double[jx,jy];
			double[,] c = new double[jx,jy];
			double[,] d = new double[jx,jy];
			double[,] e = new double[jx,jy];
			double[,] f = new double[jx,jy];
			isItInDomain = new bool[jx,jy];
			bv = new double[jx,jy];
			solution = new double[jx,jy];

			for(int k = 0; k < jx; k++)
			{
				x = xf + k*delta;
				for(int m = 0; m < jy; m++)
				{
					y = yf + m*delta;
					a[k,m]=1.0;
					b[k,m]=1.0;
					c[k,m]=1.0;
					d[k,m]=1.0;
					e[k,m]=-4.0;
					f[k,m]=source(x,y)*delta*delta;
					isItInDomain[k,m]=dtd.Check(x,y);
					if(!isItInDomain[k,m])
					{
						for(int kk = 0; kk < nbv; kk++)
						{
							if(Convert.ToDouble(al[5*kk]) < x && x < Convert.ToDouble(al[5*kk+1]) 
								&& Convert.ToDouble(al[5*kk+2]) < y && y < Convert.ToDouble(al[5*kk+3]))
								bv[k,m] = Convert.ToDouble(al[5*kk+4]);
						}
					}
				}
			}
			for(int k = 0; k < jx; k++)
			{
				for(int m = 0; m < jy; m++)
				{
					solution[k,m]=bv[k,m];
				}
			}
			sor(a,b,c,d,e,f,solution,rj);
		}

		private double[,] solution;
		public double[,] Solution
		{
			get{return solution;}
		}
	

		// numerical recipes

		private int MAXITS = 1000;
		private double EPS = 1.0e-5;

		private void sor(double[,] a, double[,] b, double[,] c, double[,] d, double[,] e,
			double[,] f, double[,] u, double rjac)
		{
			int jmax = a.GetLength(0);
			int ipass,j,jsw,l,lsw,n;
			double anorm, omega=1.0,resid;
			//	double anormf=0.0;

			//	for (j=1;j<jmax-1;j++)
			//		for (l=1;l<jmax-1;l++)
			//			anormf += Math.Abs(f[j,l]);
			for (n=1;n<=MAXITS;n++) 
			{
				anorm=0.0;
				jsw=1;
				for (ipass=1;ipass<=2;ipass++) 
				{
					lsw=jsw;
					for (j=1;j<jmax-1;j++) 
					{
						for (l=lsw;l<jmax-1;l+=2) 
						{
							if(isItInDomain[j,l])
							{
								resid=a[j,l]*u[j+1,l]
									+b[j,l]*u[j-1,l]
									+c[j,l]*u[j,l+1]
									+d[j,l]*u[j,l-1]
									+e[j,l]*u[j,l]
									-f[j,l];
								anorm += Math.Abs(resid);
								u[j,l] -= omega*resid/e[j,l];
							}
							else u[j,l] = bv[j,l];
						}
						lsw=3-lsw;
					}
					jsw=3-jsw;
					omega=(n == 1 && ipass == 1 ? 1.0/(1.0-0.5*rjac*rjac) :
						1.0/(1.0-0.25*rjac*rjac*omega));
				}
				if (anorm < EPS) return; //*anormf) return;
			}
            try { throw new Science.Error(); }
            catch (Science.Error er)
            {
                er.Write("MAXITS exceeded");
            }		
		}
	}
}