using System;

namespace Science.Mathematics.SpecialFunction
{
	/// <summary>
	/// OrthogonalPolynomial
	/// </summary>
	public class OrthogonalPolynomial
	{
		public OrthogonalPolynomial()
		{
		}
		public static double LegendreP(int l, int m, double x)
		{
			double fact,pll=0.0,pmm,pmmp1,somx2;
			int i,ll;

			if (m < 0 || m > l || Math.Abs(x) > 1.0)
                try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Bad arguments in routine plgndr");
                }
			pmm=1.0;
			if (m > 0) 
			{
				somx2=Math.Sqrt((1.0-x)*(1.0+x));
				fact=1.0;
				for (i=1;i<=m;i++) 
				{
					pmm *= -fact*somx2;
					fact += 2.0;
				}
			}
			if (l == m)
				return pmm;
			else 
			{
				pmmp1=x*(2*m+1)*pmm;
				if (l == (m+1))
					return pmmp1;
				else 
				{
					for (ll=m+2;ll<=l;ll++) 
					{
						pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m);
						pmm=pmmp1;
						pmmp1=pll;
					}
					return pll;
				}
			}
		}
		public static Complex SphericalHarmonicY(int l, int m, double theta, double phi)
		{
			double coef = Math.Sqrt((2.0*(double)l+1.0)/(4.0*Math.PI)
				*GammaRelated.Factorial1(l-m)/GammaRelated.Factorial1(l+m));
			double ct = Math.Cos(theta);
            Science.Mathematics.Complex k = new Science.Mathematics.Complex();
            k.Abs = 1.0;
            k.Arg = (double)m * phi;
			return (coef*LegendreP(l,m,ct))*k;
		}

	/*	public static double ChebyshevT(double x)
		{
			return x;
		}
		public static double ChebyshevU(double x)
		{
			return x;
		}
		public static double GegenbauerC(double x)
		{
			return x;
		}
		public static double HermiteH(double x)
		{
			return x;
		}
		public static double JacobiP(double x)
		{
			return x;
		}
		public static double LaguerreL(double x)
		{
			return x;
		}  */
	}
}