using System;

namespace Science.Mathematics.IntegralEquation
{
	/// <summary>
	/// InhomogeneousFredholmEquationSecondKind
	/// </summary>
	public class InhomogeneousFredholmEquationSecondKind
	{
        public InhomogeneousFredholmEquationSecondKind()
        {
        }
        private Function.ToLastType<double,double,double> kp;
        private Function.ToLastType<double,double> gp;
        public Function.ToLastType<double,double,double> Kernel
        {
            set { kp = value; }
        }
        public Function.ToLastType<double,double> KernelFunction
        {
            set { gp = value; }
        }
				
		private int n;
		public int NumberOfPoints
		{
			get{return n;}
			set{n=value;}
		}
		private double from, to;
		public double From
		{
			get{return from;}
			set{from=value;}
		}
		public double To
		{
			get{return to;}
			set{to=value;}
		}

        private Fred2 obj = new Fred2();
		public void Solve()
		{            
			obj.fred2(from, to, n, gp, kp);
		}
		private double at;
		public double At
		{
			get{return at;}
			set{at=value;}
		}
        public double Solution
		{
			get
			{
                return obj.fredin(at);
			}
		}		
	}


    // private class
    class Fred2
    {
        public Fred2()
        {
        }
        Function.ToLastType<double, double> g;
        Function.ToLastType<double, double, double> ak;
        double a, b;
        int n;
        public double[] t, f, w;
        public void fred2(double aa, double bb, int nn,
            Function.ToLastType<double, double> gg,
            Function.ToLastType<double, double, double> akk)
        {
            a = aa;
            b = bb;
            n = nn;
            g = gg;
            ak = akk;
            t = new double[n];
            f = new double[n];
            w = new double[n];
            double[,] omk = new Double[n, n];
            Gauleg gau = new Gauleg();
            gau.gauleg(a, b, t, w);
            double lam;
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    if (i == j) lam = 1.0;
                    else lam = 0.0;
                    omk[i, j] = lam - ak(t[i], t[j]) * w[j];
                }
                f[i] = g(t[i]);
            }
            LUdcmp lud = new LUdcmp(omk);
            lud.solve(f, f);
        }

        public double fredin(double x)
        {
            double sum = 0.0;

            for (int i = 0; i < n; i++) sum += ak(x, t[i]) * w[i] * f[i];
            return g(x) + sum;
        }
    }

    class Gauleg
    {
        public Gauleg()
        {
        }
        public void gauleg(double x1, double x2, double[] x, double[] w)
        {
            double EPS = 1.0e-14;
            int n = x.Length;
            int m, j, i;
            double z1, z, xm, xl, pp, p3, p2, p1;

            m = (n + 1) / 2;
            xm = 0.5 * (x2 + x1);
            xl = 0.5 * (x2 - x1);
            for (i = 0; i < m; i++)
            {
                z = Math.Cos(3.141592654 * (i + 1.0 - 0.25) / (n + 0.5));
                do
                {
                    p1 = 1.0;
                    p2 = 0.0;
                    for (j = 0; j < n; j++)
                    {
                        p3 = p2;
                        p2 = p1;
                        p1 = ((2.0 * (j + 1.0) - 1.0) * z * p2 - (j + 1.0 - 1.0) * p3) / (j + 1.0);
                    }
                    pp = n * (z * p1 - p2) / (z * z - 1.0);
                    z1 = z;
                    z = z1 - p1 / pp;
                } while (Math.Abs(z - z1) > EPS);
                x[i] = xm - xl * z;
                x[n - i - 1] = xm + xl * z;
                w[i] = 2.0 * xl / ((1.0 - z * z) * pp * pp);
                w[n - i - 1] = w[i];
            }
        }
    }


    class LUdcmp
    {
        int n;
        private double[,] lu;
        private int[] indx;
        private double d;
        private double[,] aref;

        public double[,] LU
        {
            get { return lu; }  // output
        }
        public int[] RowPermutation // output
        {
            get { return indx; }
        }
        public double Sign // output
        {
            get { return d; }
        }
        public LUdcmp(double[,] a)   // a input
        {
            n = a.GetLength(0);
            lu = new double[n, n];
            for (int ir = 0; ir < n; ir++)
                for (int ic = 0; ic < n; ic++) lu[ir, ic] = a[ir, ic];
            aref = a;
            indx = new int[n];
            double TINY = 1.0e-40;
            int i, imax = 0, j, k;
            double big, temp;
            double[] vv = new Double[n];
            d = 1.0;
            for (i = 0; i < n; i++)
            {
                big = 0.0;
                for (j = 0; j < n; j++)
                    if ((temp = Math.Abs(lu[i, j])) > big) big = temp;
                if (big == 0.0) try { throw new Science.Error(); }
                    catch (Science.Error e)
                    {
                        e.Write("Singular matrix in LUdcmp");
                    } 
                vv[i] = 1.0 / big;
            }
            for (k = 0; k < n; k++)
            {
                big = 0.0;
                for (i = k; i < n; i++)
                {
                    temp = vv[i] * Math.Abs(lu[i, k]);
                    if (temp > big)
                    {
                        big = temp;
                        imax = i;
                    }
                }
                if (k != imax)
                {
                    for (j = 0; j < n; j++)
                    {
                        temp = lu[imax, j];
                        lu[imax, j] = lu[k, j];
                        lu[k, j] = temp;
                    }
                    d = -d;
                    vv[imax] = vv[k];
                }
                indx[k] = imax;
                if (lu[k, k] == 0.0) lu[k, k] = TINY;
                for (i = k + 1; i < n; i++)
                {
                    temp = lu[i, k] /= lu[k, k];
                    for (j = k + 1; j < n; j++) lu[i, j] -= temp * lu[k, j];
                }
            }
        }

        public void solve(double[] b, double[] x) // b input and x output
        {
            int i, ii = 0, ip, j;
            double sum;
            if (b.Length != n || x.Length != n) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Bad sizes");
                } 
            for (i = 0; i < n; i++) x[i] = b[i];
            for (i = 0; i < n; i++)
            {
                ip = indx[i];
                sum = b[ip];
                x[ip] = x[i];
                if (ii != 0) for (j = ii - 1; j < i; j++) sum -= lu[i, j] * x[j];
                else if (sum != 0.0) ii = i + 1;
                x[i] = sum;
            }
            for (i = n - 1; i >= 0; i--)
            {
                sum = x[i];
                for (j = i + 1; j < n; j++) sum -= lu[i, j] * x[j];
                x[i] = sum / lu[i, i];
            }
        }

        public void solve(double[,] b, double[,] x) // b input and x output
        {
            int i, j, m = b.GetLength(1);
            if (b.GetLength(0) != n || x.GetLength(0) != n
                || b.GetLength(1) != x.GetLength(1)) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Bad sizes");
                } 
            double[] xx = new double[n];
            for (j = 0; j < m; j++)
            {
                for (i = 0; i < n; i++) xx[i] = b[i, j];
                solve(xx, xx);
                for (i = 0; i < n; i++) x[i, j] = xx[i];
            }
        }

        public double[,] inverse()
        {
            double[,] ainv = new double[n, n];
            for (int i = 0; i < n; i++) ainv[i, i] = 1.0;
            solve(ainv, ainv);
            return ainv;
        }

        public double det()
        {
            double dd = d;
            for (int i = 0; i < n; i++) dd *= lu[i, i];
            return dd;
        }

        public void mprove(double[] b, double[] x)
        {
            int j, i;
            double sdp;
            double[] r = new Double[n];

            for (i = 0; i < n; i++)
            {
                sdp = -b[i];
                for (j = 0; j < n; j++) sdp += aref[i, j] * x[j];
                r[i] = sdp;
            }
            solve(r, r);
            for (i = 0; i < n; i++) x[i] -= r[i];
        }
    }
}