using System;
using System.Collections;

namespace Science.Mathematics.SpecialFunction
{
	/// <summary>
	/// Hypergeometric
	/// </summary>
	public class Hypergeometric
	{
		public Hypergeometric()
		{
		}
		private static Complex aa,bb,cc,z0,dz;
	    public static Complex Hypergeometric2F1(Complex a, Complex b, Complex c, Complex z)
		{
			double EPS = 1.0e-6;
			
			Complex ans;
			Complex[] y = new Complex[2];
			double[] yy= new Double[4];

			if (z.Real*z.Real+z.Imaginary*z.Imaginary <= 0.25) 
			{
				hypser(a,b,c,z);
				ans = series;
				y[1] = deriv;
				return ans;
			}
			else if (z.Real < 0.0) z0=new Complex(-0.5,0.0);
			else if (z.Real <= 1.0) z0=new Complex(0.5,0.0);
			else z0=new Complex(0.0,z.Imaginary >= 0.0 ? 0.5 : -0.5);
			aa=a;
			bb=b;
			cc=c;
			dz=(z - z0);
			hypser(aa,bb,cc,z0);
			y[0]=series;
			y[1]=deriv;
			yy[0]=y[0].Real;
			yy[1]=y[0].Imaginary;
			yy[2]=y[1].Real;
			yy[3]=y[1].Imaginary;
			Function der = new Function(hypdrv);
			odeint(yy,4,0.0,1.0,EPS,0.1,0.0001,der);
			y[0]=new Complex(yy[0],yy[1]);
			return y[0];
		} 

		
		private static double[] hypdrv(double s, double[] yy)
		{
			Complex ONE = new Complex(1.0,0.0);
        	double[] dyyds = new Double[4];
			Complex z;
			Complex[] y = new Complex[2];
			Complex[] dyds = new Complex[2];
			y[0]=new Complex(yy[0],yy[1]);
			y[1]=new Complex(yy[2],yy[3]);
			z=(z0 + (s * dz));
			dyds[0]=(y[1] * dz);
			dyds[1]=((((aa*bb)*y[0]) - 
				((cc - (((aa + bb) + ONE) * z)) * y[1])) *
				(dz / (z * (ONE - z))));
			dyyds[0]=dyds[0].Real;
			dyyds[1]=dyds[0].Imaginary;
			dyyds[2]=dyds[1].Real;
			dyyds[3]=dyds[1].Imaginary;
			return dyyds;
		}  

		private static Complex series, deriv;
		private static void hypser(Complex a, Complex b,
			Complex c, Complex z)
		{
			Complex ONE = new Complex(1.0,0.0);
			int n;
			Complex aa,bb,cc,fac,temp;

			deriv = new Complex(0.0,0.0);
			fac=new Complex(1.0,0.0);
			temp=fac;
			aa=a;
			bb=b;
			cc=c;
			for (n=1;n<=1000;n++) 
			{
				fac=(fac * (aa * (bb / cc)));
				deriv.Real+=fac.Real;
				deriv.Imaginary+=fac.Imaginary;
				fac=(fac * ((1.0/n)* z));
				series=(temp + fac);
				if (series.Real == temp.Real && series.Imaginary == temp.Imaginary) return;
				temp= series;
				aa=(aa + ONE);
				bb=(bb + ONE);
				cc=(cc + ONE);

			}
            try { throw new Science.Error(); }
            catch (Science.Error e)
            {
                e.Write("convergence failure in hypser");
            }
		}   		

		private static ArrayList xp = new ArrayList(1);
		private static ArrayList[] yp;
		private static double dxsav=0.0;
		private static int nok, nbad;
		
    	private static void odeint(double[] ystart, int nvar, double x1, double x2, double eps, double h1,
			double hmin, Function derivs)
		{
			int MAXSTP = 10000;
			double TINY = 1.0e-30;

			yp = new ArrayList[nvar];
			for(int j = 0; j < nvar; j++) yp[j] = new ArrayList(1);

			int nstp,i;
			double xsav,x,h;
			double[] yscal,y,dydx;
			yscal=new Double[nvar];
			y=new Double[nvar];
			dydx=new Double[nvar];
			x=x1;
			h=Math.Sign(x2-x1)*h1;
			nok = nbad = 0;
			for (i=0;i<nvar;i++) y[i]=ystart[i];
			xsav=x-dxsav*2.0;
			for (nstp=1;nstp<=MAXSTP;nstp++) 
			{
				dydx = derivs(x,y);
				for (i=0;i<nvar;i++)
					yscal[i]=Math.Abs(y[i])+Math.Abs(dydx[i]*h)+TINY;
				if (Math.Abs(x-xsav) > Math.Abs(dxsav)) 
				{
					xp.Add(x);
					for (i=0;i<nvar;i++) yp[i].Add(y[i]);
					xsav=x;
				}
				if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x;
				rkqs(y,dydx,nvar,x,h,eps,yscal,derivs);
				x = xout;
				if (hdid == h) ++nok; else ++nbad;
				if ((x-x2)*(x2-x1) >= 0.0) 
				{
					for (i=0;i<nvar;i++) ystart[i]=y[i];
					xp.Add(x);
					for (i=0;i<nvar;i++) yp[i].Add(y[i]);
					return;
				}
				if (Math.Abs(hnext) <= hmin)
                    try { throw new Science.Error(); }
                    catch (Science.Error e)
                    {
                        e.Write("Step size too small in odeint");
                    } 
				h=hnext;
			}
            try { throw new Science.Error(); }
            catch (Science.Error e)
            {
                e.Write("Too many steps in routine odeint");
            } 
		}



		private static double xout, hdid, hnext;
		private static void rkqs(double[] y, double[] dydx, int n, double x, double htry, double eps,
			double[] yscal, Function derivs)
		{
			double SAFETY = 0.9;
	    	double PGROW = -0.2;
	    	double PSHRNK = -0.25;
	    	double ERRCON = 1.89e-4;

			int i;
			double errmax,h,xnew;
			double[] yerr = new Double[n];
			double[] ytemp = new Double[n];
			h=htry;
			for (;;) 
			{
				rkck(y,dydx,n,x,h,derivs);
				ytemp = yout;
				errmax=0.0;
				for (i=0;i<n;i++) errmax=Math.Max(errmax,Math.Abs(yerr[i]/yscal[i]));
				errmax /= eps;
				if (errmax > 1.0) 
				{
					h=SAFETY*h*Math.Pow(errmax,PSHRNK);
					if (h < 0.1*h) h *= 0.1;
					xnew=x+h;
                    if (xnew == x) try { throw new Science.Error(); }
                        catch (Science.Error e)
                        {
                            e.Write("stepsize underflow in rkqs");
                        } 
					continue;
				} 
				else 
				{
					if (errmax > ERRCON) hnext=SAFETY*h*Math.Pow(errmax,PGROW);
					else hnext=5.0*h;
					xout = x + (hdid=h);
					for (i=0;i<n;i++) y[i]=ytemp[i];
					break;
				}
			}
		}



		private delegate double[] Function(double y, double[] x);
		private static double[] yout, yerr;
		private static void rkck(double[] y, double[] dydx, int n, double x, double h, 
			Function derivs)
		{
			yout = new Double[n];
			yerr = new Double[n];
			int i;
			double a2=0.2,a3=0.3,a4=0.6,a5=1.0,a6=0.875,b21=0.2,
				b31=3.0/40.0,b32=9.0/40.0,b41=0.3,b42 = -0.9,b43=1.2,
				b51 = -11.0/54.0, b52=2.5,b53 = -70.0/27.0,b54=35.0/27.0,
				b61=1631.0/55296.0,b62=175.0/512.0,b63=575.0/13824.0,
				b64=44275.0/110592.0,b65=253.0/4096.0,c1=37.0/378.0,
				c3=250.0/621.0,c4=125.0/594.0,c6=512.0/1771.0,
				dc5 = -277.0/14336.0;
			double dc1=c1-2825.0/27648.0,dc3=c3-18575.0/48384.0,
				dc4=c4-13525.0/55296.0,dc6=c6-0.25;
			double[] ak2 = new Double[n];
			double[] ak3 = new Double[n];
			double[] ak4 = new Double[n];
			double[] ak5 = new Double[n];
			double[] ak6 = new Double[n];
			double[] ytemp = new Double[n];
			
			for (i=0;i<n;i++)
				ytemp[i]=y[i]+b21*h*dydx[i];
			ak2 = derivs(x+a2*h,ytemp);
			for (i=0;i<n;i++)
				ytemp[i]=y[i]+h*(b31*dydx[i]+b32*ak2[i]);
			ak3 = derivs(x+a3*h,ytemp);
			for (i=0;i<n;i++)
				ytemp[i]=y[i]+h*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);
			ak4 = derivs(x+a4*h,ytemp);
			for (i=0;i<n;i++)
				ytemp[i]=y[i]+h*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]);
			ak5 = derivs(x+a5*h,ytemp);
			for (i=0;i<n;i++)
				ytemp[i]=y[i]+h*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]);
			ak6 = derivs(x+a6*h,ytemp);
			for (i=0;i<n;i++)
				yout[i]=y[i]+h*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);
			for (i=0;i<n;i++)
				yerr[i]=h*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]);
		}  


	/*	public static double AppellF1(double x)
		{
			return x;
		}
		public static double ChebyshevT(double x)
		{
			return x;
		}
		public static double ChebyshevU(double x)
		{
			return x;
		}
		public static double CoshIntegral(double x)
		{
			return x;
		}
		public static double CosIntegral(double x)
		{
			return x;
		}
		public static double EllipticE(double x)
		{
			return x;
		}
		public static double EllipticK(double x)
		{
			return x;
		}
		public static double Erf(double x)
		{
			return x;
		}
		public static double Erfc(double x)
		{
			return x;
		}
		public static double Erfi(double x)
		{
			return x;
		}
		public static double ExpIntegralE(double x)
		{
			return x;
		}
		public static double ExpIntegralEi(double x)
		{
			return x;
		}
		public static double GegenbauerC(double x)
		{
			return x;
		}
		public static double HermiteH(double x)
		{
			return x;
		}
		public static double Hypergeometric0F1(double x)
		{
			return x;
		}
		public static double Hypergeometric0F1Regularized(double x)
		{
			return x;
		}
		public static double Hypergeometric1F1(double x)
		{
			return x;
		}
		public static double Hypergeometric1F1Regularized(double x)
		{
			return x;
		}
		public static double Hypergeometric2F1Regularized(double x)
		{
			return x;
		}
		public static double HypergeometricPFQ(double x)
		{
			return x;
		}
		public static double HypergeometricPFQRegularized(double x)
		{
			return x;
		}
		public static double HypergeometricU(double x)
		{
			return x;
		}
		public static double InverseBetaRegularized(double x)
		{
			return x;
		}
		public static double InverseErf(double x)
		{
			return x;
		}
		public static double InverseErfc(double x)
		{
			return x;
		}
		public static double InverseGammaRegularized(double x)
		{
			return x;
		}
		public static double JacobiP(double x)
		{
			return x;
		}
		public static double LaguerreL(double x)
		{
			return x;
		}
		public static double LegendreP(double x)
		{
			return x;
		}
		public static double LegendreQ(double x)
		{
			return x;
		}
		public static double LogIntegral(double x)
		{
			return x;
		}
		public static double Meijer(double x)
		{
			return x;
		}
		public static double ProductLog(double x)
		{
			return x;
		}
		public static double SinhIntegral(double x)
		{
			return x;
		}
		public static double SinIntegral(double x)
		{
			return x;
		}
		public static double ShpericalHarmonicY(double x)
		{
			return x;
		}
		public static double StruveH(double x)
		{
			return x;
		}
		public static double StruveL(double x)
		{
			return x;
		}  */
	}
}