using System;
using Science.Mathematics;
using L=Science.Mathematics.OrdinaryDifferentialEquation;

namespace ScienceTest.MathematicsTest.OrdinaryDifferentialEquationTest
{
	/// <summary>
	/// SpheroidalHarmonics
	/// </summary>
	public class SpheroidalHarmonics
	{
		public SpheroidalHarmonics()
		{
		}
		private string result;
		public string Result
		{
			get{return result;}
		}

		private int m=2,n=2;   // input m and n
		private double c2=0.1; // input c^2

		private double dx=0.0001,gmma=1.0;  

		public void Compute()
		{
            int q1 = n;
            for (int i = 1; i <= m; i++) gmma *= -0.5 * (n + i) * (q1--) / (double)i;

			double start = -1.0 + dx;
			double end = 0.0;
            Function.ToLastType<double, double[], double[]> load
                = new Function.ToLastType<double, double[], double[]>(Load);
            Function.ToLastType<double, double[], double[], double[]> score
                = new Function.ToLastType<double, double[], double[], double[]>(Score);
			L.BoundaryCondition bc = new L.BoundaryCondition();
            bc.Load = load;
            bc.Score = score;
            bc.Start = start;
            bc.End = end;
    
			L.EquationWithBoundaryCondition eq 
				= new L.EquationWithBoundaryCondition();
            Function.ToLastType<double, double[], double[]> func
                = new Function.ToLastType<double, double[], double[]>(Derivs);
            eq.FirstDerivativeFunction = func;
            eq.Condition = bc;
            eq.HowMany = 20;

			double[] initialParameter = new Double[1];
			initialParameter[0]=n*(n+1)-m*(m+1)+c2/2.0;
			eq.InitialGuess = initialParameter;
            
			eq.Solve();

            for (int k = 0; k < eq.IndependentVariable.Count; k++)
                result += k + "\t" + Convert.ToString(eq.IndependentVariable[k]) + "\t" +
					Convert.ToString(eq.Solution[0][k])+"\t"+
					Convert.ToString(eq.Solution[1][k])+"\r\n";
				result += Convert.ToString(eq.Solution[2][0])+"\r\n";
		}		
	
		private double[] Load(double start, double[] p)
		{
			double[] y = new Double[3];
			double y1 = ((n-m)%2 == 1 ? -gmma : gmma);
			y[2] = p[0];
			y[1] = -(y[2]-c2)*y1/(2.0*(m+1));
			y[0] = y1+y[1]*dx;
			return y;
		}
		private double[] Score(double end, double[] v, double[] p)
		{
			double[] f = new Double[p.Length];
			f[0]=((n-m)%2 == 1 ? v[0] : v[1]);
			return f;
		}
		private double[] Derivs(double x, double[] y)
		{
			double[] dydx = new Double[y.Length];
			dydx[0]=y[1];
			dydx[1]=(2.0*x*(m+1.0)*y[1]-(y[2]-c2*x*x)*y[0])/(1.0-x*x);
			dydx[2]=0.0;
			return dydx;
		} 
	}
}
/*
0	-0.9999	3.0000042866843	0.042866843029173
1	-0.949905	3.00209419580067	0.0407415338111994
2	-0.89991	3.00407797644603	0.0386170783529632
3	-0.849915	3.00595546398091	0.0364894135296038
4	-0.79992	3.00772650235836	0.0343587156256431
5	-0.749925	3.00939094435339	0.0322251612732701
6	-0.69993	3.01094865158012	0.030088927425051
7	-0.649935	3.01239949450683	0.0279501913415
8	-0.59994	3.01374335247117	0.0258091305650124
9	-0.549945	3.01498011369385	0.0236659229007263
10	-0.49995	3.01610967529088	0.0215207463999522
11	-0.449955	3.01713194328569	0.0193737793358981
12	-0.39996	3.01804683262076	0.0172252001803616
13	-0.349965	3.01885426716412	0.015075187597363
14	-0.29997	3.01955417972089	0.0129239204096349
15	-0.249975	3.02014651204036	0.0107715775811334
16	-0.19998	3.02063121482215	0.00861833819860957
17	-0.149985	3.02100824772194	0.00646438145066932
18	-0.0999899999999998	3.02127757935603	0.00430988660738155
19	-0.0499949999999998	3.021439187305	0.00215503299991733
20	0	3.0214930581161	6.72205346941013E-18
0.0142663139416539
*/