using System;

namespace Science.Mathematics.SpecialFunction
{
	/// <summary>
	/// Elliptic
	/// </summary>
	public class Elliptic
	{
		public Elliptic()
		{
		}	
		public static double CarlsonFirstKind(double x, double y, double z)
		{
			double ERRTOL = 0.08;
			double TINY = 1.5e-38;
			double BIG = 3.0e37;
			double THIRD = 1.0/3.0;
			double C1 = 1.0/24.0;
			double C2 = 0.1;
			double C3 = 3.0/44.0;
			double C4 = 1.0/14.0;

			double alamb,ave,delx,dely,delz,e2,e3,sqrtx,sqrty,sqrtz,xt,yt,zt;

			if (Math.Min(Math.Min(x,y),z) < 0.0 || Math.Min(Math.Min(x+y,x+z),y+z) < TINY ||
				Math.Max(Math.Max(x,y),z) > BIG)
                try { throw new Science.Error(); }
                    catch (Science.Error e)
                    {
                        e.Write("invalid arguments in rf");
                    }
			xt=x;
			yt=y;
			zt=z;
			do 
			{
				sqrtx=Math.Sqrt(xt);
				sqrty=Math.Sqrt(yt);
				sqrtz=Math.Sqrt(zt);
				alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz;
				xt=0.25*(xt+alamb);
				yt=0.25*(yt+alamb);
				zt=0.25*(zt+alamb);
				ave=THIRD*(xt+yt+zt);
				delx=(ave-xt)/ave;
				dely=(ave-yt)/ave;
				delz=(ave-zt)/ave;
			} while (Math.Max(Math.Max(Math.Abs(delx),Math.Abs(dely)),Math.Abs(delz)) > ERRTOL);
			e2=delx*dely-delz*delz;
			e3=delx*dely*delz;
			return (1.0+(C1*e2-C2-C3*e3)*e2+C4*e3)/Math.Sqrt(ave);
		}

		public static double CarlsonSecondKind(double x, double y, double z)
		{
			double ERRTOL = 0.05;
			double TINY = 1.0e-25;
			double BIG = 4.5e21;
			double C1 = 3.0/14.0;
			double C2 = 1.0/6.0;
			double C3 = 9.0/22.0;
			double C4 = 3.0/26.0;
			double C5 = 0.25*9.0/22.0;
			double C6 = 1.5*3.0/26.0;

            double alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,sqrtx,sqrty,
				sqrtz,sum,xt,yt,zt;

			if (Math.Min(x,y) < 0.0 || Math.Min(x+y,z) < TINY || Math.Max(Math.Max(x,y),z) > BIG)
                try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("invalid arguments in rd");
                } 
            xt=x;
			yt=y;
			zt=z;
			sum=0.0;
			fac=1.0;
			do 
			{
				sqrtx=Math.Sqrt(xt);
				sqrty=Math.Sqrt(yt);
				sqrtz=Math.Sqrt(zt);
				alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz;
				sum += fac/(sqrtz*(zt+alamb));
				fac=0.25*fac;
				xt=0.25*(xt+alamb);
				yt=0.25*(yt+alamb);
				zt=0.25*(zt+alamb);
				ave=0.2*(xt+yt+3.0*zt);
				delx=(ave-xt)/ave;
				dely=(ave-yt)/ave;
				delz=(ave-zt)/ave;
			} while (Math.Max(Math.Max(Math.Abs(delx),Math.Abs(dely)),Math.Abs(delz)) > ERRTOL);
			ea=delx*dely;
			eb=delz*delz;
			ec=ea-eb;
			ed=ea-6.0*eb;
			ee=ed+ec+ec;
			return 3.0*sum+fac*(1.0+ed*(-C1+C5*ed-C6*delz*ee)
				+delz*(C2*ee+delz*(-C3*ec+delz*C4*ea)))/(ave*Math.Sqrt(ave));
		}

		
    	public static double CarlsonThirdKind(double x, double y, double z, double p)
		{
			double ERRTOL = 0.05;
			double TINY = 2.5e-13;
			double BIG = 9.0e11;
			double C1 = (3.0/14.0);
			double C2 = (1.0/3.0);
			double C3 = (3.0/22.0);
			double C4 = (3.0/26.0);
			double C5 = (0.75*3.0/22.0);
			double C6 = (1.5*3.0/26.0);
			double C7 = (0.5*1.0/3.0);
			double C8 = (6.0/22.0);

			double a=0,alamb,alpha,ans,ave,b=0,beta,delp,delx,dely,delz,ea,eb,ec,
				ed,ee,fac,pt,rcx=0,rho,sqrtx,sqrty,sqrtz,sum,tau,xt,yt,zt;

			if (Math.Min(Math.Min(x,y),z) < 0.0 || Math.Min(Math.Min(x+y,x+z),Math.Min(y+z,Math.Abs(p))) < TINY
				|| Math.Max(Math.Max(x,y),Math.Max(z,Math.Abs(p))) > BIG)
                try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("invalid arguments in rj");
                } 
            sum=0.0;
			fac=1.0;
			if (p > 0.0) 
			{
				xt=x;
				yt=y;
				zt=z;
				pt=p;
			} 
			else 
			{
				xt=Math.Min(Math.Min(x,y),z);
				zt=Math.Max(Math.Max(x,y),z);
				yt=x+y+z-xt-zt;
				a=1.0/(yt-p);
				b=a*(zt-yt)*(yt-xt);
				pt=yt+b;
				rho=xt*zt/yt;
				tau=p*pt/yt;
				rcx=CarlsonDegenerate(rho,tau);
			}
			do 
			{
				sqrtx=Math.Sqrt(xt);
				sqrty=Math.Sqrt(yt);
				sqrtz=Math.Sqrt(zt);
				alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz;
				alpha=(pt*(sqrtx+sqrty+sqrtz)+sqrtx*sqrty*sqrtz)*(pt*(sqrtx+sqrty+sqrtz)+sqrtx*sqrty*sqrtz);
				beta=pt*(pt+alamb)*(pt+alamb);
				sum += fac*CarlsonDegenerate(alpha,beta);
				fac=0.25*fac;
				xt=0.25*(xt+alamb);
				yt=0.25*(yt+alamb);
				zt=0.25*(zt+alamb);
				pt=0.25*(pt+alamb);
				ave=0.2*(xt+yt+zt+pt+pt);
				delx=(ave-xt)/ave;
				dely=(ave-yt)/ave;
				delz=(ave-zt)/ave;
				delp=(ave-pt)/ave;
			} while (Math.Max(Math.Max(Math.Abs(delx),Math.Abs(dely)),
				Math.Max(Math.Abs(delz),Math.Abs(delp))) > ERRTOL);
			ea=delx*(dely+delz)+dely*delz;
			eb=delx*dely*delz;
			ec=delp*delp;
			ed=ea-3.0*ec;
			ee=eb+2.0*delp*(ea-ec);
			ans=3.0*sum+fac*(1.0+ed*(-C1+C5*ed-C6*ee)+eb*(C7+delp*(-C8+delp*C4))
				+delp*ea*(C2-delp*C3)-C2*delp*ec)/(ave*Math.Sqrt(ave));
			if (p <= 0.0) ans=a*(b*ans+3.0*(rcx-CarlsonFirstKind(xt,yt,zt)));
			return ans;
		}

		
    	public static double CarlsonDegenerate(double x, double y)
		{
			double ERRTOL = 0.04;
			double TINY = 1.69e-38;
			double BIG = 3.0e37;
			double COMP1 = 2.236/1.3e-19;
			double COMP2 = 1.69e-38*3.0e37*1.69e-38*3.0e37/25.0;
			double THIRD = (1.0/3.0);
			double C1 = 0.3;
			double C2 = (1.0/7.0);
			double C3 = 0.375;
			double C4 = (9.0/22.0);

			double alamb,ave,s,w,xt,yt;
			if (x < 0.0 || y == 0.0 || (x+Math.Abs(y)) < TINY || (x+Math.Abs(y)) > BIG ||
				(y<-COMP1 && x > 0.0 && x < COMP2))
                try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("invalid arguments in rc");
                } 
			if (y > 0.0) 
			{
				xt=x;
				yt=y;
				w=1.0;
			} 
			else 
			{
				xt=x-y;
				yt = -y;
				w=Math.Sqrt(x)/Math.Sqrt(xt);
			}
			do 
			{
				alamb=2.0*Math.Sqrt(xt)*Math.Sqrt(yt)+yt;
				xt=0.25*(xt+alamb);
				yt=0.25*(yt+alamb);
				ave=THIRD*(xt+yt+yt);
				s=(yt-ave)/ave;
			} while (Math.Abs(s) > ERRTOL);
			return w*(1.0+s*s*(C1+s*(C2+s*(C3+s*C4))))/Math.Sqrt(ave);
		}

		public static double LegendreFirstKindF(double phi, double ak)
		{
			double s;
			s=Math.Sin(phi);
			return s*CarlsonFirstKind(Math.Cos(phi)*Math.Cos(phi),(1.0-s*ak)*(1.0+s*ak),1.0);
		} 

		public static double LegendreSecondKindE(double phi, double ak)
		{
			double cc,q,s;

			s=Math.Sin(phi);
			cc=Math.Cos(phi)*Math.Cos(phi);
			q=(1.0-s*ak)*(1.0+s*ak);
			return s*(CarlsonFirstKind(cc,q,1.0)-(s*ak*s*ak)*CarlsonSecondKind(cc,q,1.0)/3.0);
		}   

		public static double LegendreThirdKindPI(double phi, double en, double ak)
		{
			double cc,enss,q,s;

			s=Math.Sin(phi);
			enss=en*s*s;
			cc=Math.Cos(phi)*Math.Cos(phi);
			q=(1.0-s*ak)*(1.0+s*ak);
			return s*(CarlsonFirstKind(cc,q,1.0)-enss*CarlsonThirdKind(cc,q,1.0,1.0+enss)/3.0);
		} 

		public static double JacobianSn(double u, double k)
		{
			return sncndn(u,1.0-k*k,0);
		}
		public static double JacobianCn(double u, double k)
		{
			return sncndn(u,1.0-k*k,1);
		}
		public static double JacobianDn(double u, double k)
		{
			return sncndn(u,1.0-k*k,2);
		}

		private static double sncndn(double uu, double emmc, int returnwhat)
		{
			double CA = 0.0003;
			double sn, cn, dn;
			double a,b,c=0.0,d=0.0,emc,u;
			double[] em = new Double[13];
			double[] en = new Double[13];
			int i,ii,l=0;
			bool bo = false;

			emc=emmc;
			u=uu;
			if (emc != 0.0) 
			{
				if(emc < 0.0) bo = true;
				if (bo) 
				{
					d=1.0-emc;
					emc /= -1.0/d;
					u *= (d=Math.Sqrt(d));
				}
				a=1.0;
				dn=1.0;
				for (i=1;i<=13;i++) 
				{
					l=i;
					em[i-1]=a;
					en[i-1]=(emc=Math.Sqrt(emc));
					c=0.5*(a+emc);
					if (Math.Abs(a-emc) <= CA*a) break;
					emc *= a;
					a=c;
				}
				u *= c;
				sn=Math.Sin(u);
				cn=Math.Cos(u);
				if (sn != 0.0) 
				{
					a=cn/sn;
					c *= a;
					for (ii=l;ii>=1;ii--) 
					{
						b=em[ii-1];
						a *= c;
						c *= dn;
						dn=(en[ii-1]+a)/(b+a);
						a=c/b;
					}
					a=1.0/Math.Sqrt(c*c+1.0);
					sn=(sn >= 0.0 ? a : -a);
					cn=c*sn;
				}
				if (bo) 
				{
					a=dn;
					dn=cn;
					cn=a;
					sn /= d;
				}
			} 
			else 
			{
				cn=1.0/Math.Cosh(u);
				dn=cn;
				sn=Math.Tanh(u);
			}
			
			if(returnwhat==0) return sn;
			else if(returnwhat==1) return cn;
			else return dn;
		} 


	/*	public static double DedekindEta(double x)
		{
			return x;
		}
		public static double EllipticExp(double x)
		{
			return x;
		}
		public static double EllipticTheta(double x)
		{
			return x;
		}
		public static double EllipticThetaPrime(double x)
		{
			return x;
		}
		public static double InverseEllipticNomeQ(double x)
		{
			return x;
		}
		public static double InverseJacobiSn(double x)
		{
			return x;
		}
		public static double InverseWeiserstrassP(double x)
		{
			return x;
		}
		public static double JacobiAmplitude(double x)
		{
			return x;
		}
		public static double KleinInvariantJ(double x)
		{
			return x;
		}
		public static double ModularLambda(double x)
		{
			return x;
		}
		public static double WeierstrassHalfPeriods(double x)
		{
			return x;
		}
		public static double WeierstrassP(double x)
		{
			return x;
		}
		public static double WeierstrassPPrime(double x)
		{
			return x;
		}
		public static double WeierstrassSigma(double x)
		{
			return x;
		}
		public static double WeierstrassZeta(double x)
		{
			return x;
		}  */		
	}
}