using System;

namespace Science.Mathematics.SpecialFunction
{
	/// <summary>
	/// Elliptic Integral Function 
	/// </summary>
	public class Integral
	{
		public Integral()
		{
		}
		public static double FresnelS(double x)
		{
			return frenel(x,1);
		}
		public static double FresnelC(double x)
		{
			return frenel(x,0);
		}
		private static double frenel(double x, int returnwhat)
		{
			double EPS = 6.0e-8;
			int MAXIT = 100;
			double FPMIN = 1.0e-30;
			double XMIN = 1.5; 
			double PI = Math.PI;
			double PIBY2 = Math.PI/2.0;
			bool TRUE = true;
			Complex ONE = new Complex(1.0,0.0);

			double s, c;
			int k,n;
			bool odd;
			double a,ax,fact,pix2,sign,sum,sumc,sums,term,test;
			Complex b,cc,d,h,del,cs;

			ax=Math.Abs(x);
			if (ax < Math.Sqrt(FPMIN)) 
			{
				s=0.0;
				c=ax;
			} 
			else if (ax <= XMIN) 
			{
				sum=sums=0.0;
				sumc=ax;
				sign=1.0;
				fact=PIBY2*ax*ax;
				odd=TRUE;
				term=ax;
				n=3;
				for (k=1;k<=MAXIT;k++) 
				{
					term *= fact/k;
					sum += sign*term/n;
					test=Math.Abs(sum)*EPS;
					if (odd) 
					{
						sign = -sign;
						sums=sum;
						sum=sumc;
					} 
					else 
					{
						sumc=sum;
						sum=sums;
					}
					if (term < test) break;
					odd=!odd;
					n += 2;
				}
                if (k > MAXIT) try { throw new Science.Error(); }
                    catch (Science.Error e)
                    {
                        e.Write("series failed in frenel");
                    } 
				s=sums;
				c=sumc;
			} 
			else 
			{
				pix2=PI*ax*ax;
				b=new Complex(1.0,-pix2);
				cc=new Complex(1.0/FPMIN,0.0);
				Complex t;
				d=h=(ONE / b);
				n = -1;
				for (k=2;k<=MAXIT;k++) 
				{
					n += 2;
					a = -n*(n+1);
					t = new Complex(4.0,0.0);
					b=(b + t);
					d=(ONE / ((a*d) + b));
					t = new Complex(a,0.0);
					cc=(b + (t / cc));
					del=(cc * d);
					h=(h * del);
					if (Math.Abs(del.Real-1.0)+Math.Abs(del.Imaginary) < EPS) break;
				}
                if (k > MAXIT) try { throw new Science.Error(); }
                    catch (Science.Error e)
                    {
                        e.Write("cf failed in frenel");
                    } 
				t = new Complex(ax,-ax);
				h=(t * h);
				t = new Complex(0.5,0.5);
				Complex tt = 
					new Complex(Math.Cos(0.5*pix2),Math.Sin(0.5*pix2));
				cs=(t * (ONE - (tt * h)));
				c=cs.Real;
				s=cs.Imaginary;
			}
			if (x < 0.0) 
			{
				c = -c;
				s = -s;
			}
			if(returnwhat==1) return s;
			else return c;
		}

		public static double CosineCi(double x)
		{
			return cisi(x,0);
		}
		public static double SineSi(double x)
		{
			return cisi(x,1);
		}
		private static double cisi(double x, int returnwhat)
		{
			double EPS = 6.0e-8;
			double EULER = 0.57721566;
			int MAXIT = 100;
			double PIBY2 = 1.5707963;
			double FPMIN = 1.0e-30;
			double TMIN = 2.0;
			bool TRUE = true;
			Complex ONE = new Complex(1.0,0.0);

			double ci, si;
			int i,k;
			bool odd;
			double a,err,fact,sign,sum,sumc,sums,t,term;
			Complex h,b,c,d,del,tc;

			t=Math.Abs(x);
			if (t == 0.0) 
			{
				si=0.0;
				ci = -1.0/FPMIN;
				if(returnwhat==1) return si;
				else return ci;
			}
			if (t > TMIN) 
			{
				b=new Complex(1.0,t);
				c=new Complex(1.0/FPMIN,0.0);
				d=h=(ONE / b);
				for (i=2;i<=MAXIT;i++) 
				{
					a = -(i-1)*(i-1);
					tc = new Complex(2.0,0.0);
					b=(b + tc);
					d=(ONE / ((a*d) + b));
					tc = new Complex(a,0.0);
					c=(b + (tc / c));
					del=(c * d);
					h=(h * del);
					if (Math.Abs(del.Real-1.0)+Math.Abs(del.Imaginary) < EPS) break;
				}
                if (i > MAXIT) try { throw new Science.Error(); }
                    catch (Science.Error e)
                    {
                        e.Write("cf failed in cisi");
                    }
				tc = new Complex(Math.Cos(t),-Math.Sin(t));
				h=(tc * h);
				ci = -h.Real;
				si=PIBY2+h.Imaginary;
			} 
			else 
			{
				if (t < Math.Sqrt(FPMIN)) 
				{
					sumc=0.0;
					sums=t;
				} 
				else 
				{
					sum=sums=sumc=0.0;
					sign=fact=1.0;
					odd=TRUE;
					for (k=1;k<=MAXIT;k++) 
					{
						fact *= t/k;
						term=fact/k;
						sum += sign*term;
						err=term/Math.Abs(sum);
						if (odd) 
						{
							sign = -sign;
							sums=sum;
							sum=sumc;
						} 
						else 
						{
							sumc=sum;
							sum=sums;
						}
						if (err < EPS) break;
						odd=!odd;
					}
                    if (k > MAXIT) try { throw new Science.Error(); }
                        catch (Science.Error e)
                        {
                            e.Write("maxits exceeded in cisi");
                        } 
				}
				si=sums;
				ci=sumc+Math.Log(t)+EULER;
			}
			if (x < 0.0) si = -si;
			if(returnwhat==1) return si;
			else return ci;
		}	


		public static double Dawson(double x)
		{
			
			int NMAX = 6;
			double H = 0.4;
			double A1 = 2.0/3.0;
			double A2 = 0.4;
			double A3 = 2.0/7.0;
			double[] c = new Double[7];
			int init = 0;

			int i,n0;
			double d1,d2,e1,e2,sum,x2,xp,xx,ans;

			if (init == 0) 
			{
				init=1;
				for (i=1;i<=NMAX;i++) c[i]=Math.Exp(-((2.0*i-1.0)*H)*((2.0*i-1.0)*H));
			}
			if (Math.Abs(x) < 0.2) 
			{
				x2=x*x;
				ans=x*(1.0-A1*x2*(1.0-A2*x2*(1.0-A3*x2)));
			} 
			else 
			{
				xx=Math.Abs(x);
				n0=2*(int)(0.5*xx/H+0.5);
				xp=xx-n0*H;
				e1=Math.Exp(2.0*xp*H);
				e2=e1*e1;
				d1=n0+1;
				d2=d1-2.0;
				sum=0.0;
				for (i=1;i<=NMAX;i++,d1+=2.0,d2-=2.0,e1*=e2)
					sum += c[i]*(e1/d1+1.0/(d2*e1));
				ans=0.5641895835*Math.Exp(-xp*xp)*Math.Sign(x)*sum;
			}
			return ans;
		} 

	/*	public static double EllipticE(double x)
		{
			return x;
		}
		public static double EllipticF(double x)
		{
			return x;
		}
		public static double EllipticK(double x)
		{
			return x;
		}
		public static double EllipticLog(double x)
		{
			return x;
		}
		public static double EllipticNomeQ(double x)
		{
			return x;
		}
		public static double EllipticPi(double x)
		{
			return x;
		}
		public static double JacobiZeta(double x)
		{
			return x;
		}  */
	}
}