using System;

namespace Science.Mathematics
{
	/// <summary>
	/// The complex numbers are the field of numbers of 
	/// the form x+iy, where x and y are real numbers 
	/// and i is the imaginary unit equal to the square root 
	/// of -1, \sqrt{-1}. 
	/// </summary>
    public class Complex
    {
        private double re, im, abs, arg;
        public double Real
        {
            get { return re; }
            set { re = value; FindAngle(); }
        }
        public double Imaginary
        {
            get { return im; }
            set { im = value; FindAngle(); }
        }
        public double Abs
        {
            get { return abs; }
            set { abs = value; re = abs * Math.Cos(Arg); im = abs * Math.Sin(Arg); }
        }
        public double Arg
        {
            get { return arg; }
            set { arg = value; re = abs * Math.Cos(Arg); im = abs * Math.Sin(Arg); }
        }
        public Complex()
        {
        }
        public Complex(double x, double y)
        {
            re = x;
            im = y;
            FindAngle();
        }
        private void FindAngle()
        {
            abs = Math.Sqrt(re * re + im * im);
            if (re > 0.0 & im > 0.0)
                arg = Math.Atan(im / re);
            else if (re > 0.0 & im < 0.0)
                arg = 2.0 * Math.PI + Math.Atan(im / re);
            else if (re < 0.0 & im > 0.0)
                arg = Math.PI + Math.Atan(im / re);
            else arg = Math.PI + Math.Atan(im / re);
        }

        public static Complex operator +(Complex a, Complex b)
        {
            Complex c = new Complex();
            c.Real = a.Real + b.Real;
            c.Imaginary = a.Imaginary + b.Imaginary;
            return c;
        }
        public static Complex operator -(Complex a, Complex b)
        {
            Complex c = new Complex();
            c.Real = a.Real - b.Real;
            c.Imaginary = a.Imaginary - b.Imaginary;
            return c;
        }
        public static Complex operator -(Complex a)
        {
            Complex c = new Complex();
            c.Real = - a.Real;
            c.Imaginary = - a.Imaginary;
            return c;
        }
        public static Complex operator *(Complex a, Complex b)
        {
            Complex c = new Complex();
            c.Real = a.Real * b.Real - a.Imaginary * b.Imaginary;
            c.Imaginary = a.Imaginary * b.Real + a.Real * b.Imaginary;
            return c;
        }
        public static Complex operator /(Complex a, Complex b)
        {
            Complex c = new Complex();
            double div = b.Real * b.Real + b.Imaginary * b.Imaginary;
            c.Real = (a.Real * b.Real + a.Imaginary * b.Imaginary) / div;
            c.Imaginary = (a.Imaginary * b.Real - a.Real * b.Imaginary) / div;
            return c;
        }
        public static Complex operator +(Complex a, double b)
        {
            Complex c = new Complex();
            c.Real = a.Real + b;
            c.Imaginary = a.Imaginary;
            return c;
        }
        public static Complex operator +(double a, Complex b)
        {
            Complex c = new Complex();
            c.Real = a + b.Real;
            c.Imaginary = b.Imaginary;
            return c;
        }
        public static Complex operator -(Complex a, double b)
        {
            Complex c = new Complex();
            c.Real = a.Real - b;
            c.Imaginary = a.Imaginary;
            return c;
        }
        public static Complex operator -(double a, Complex b)
        {
            Complex c = new Complex();
            c.Real = a - b.Real;
            c.Imaginary = - b.Imaginary;
            return c;
        }
        public static Complex operator *(double x, Complex a)
        {
            Complex c = new Complex();
            c.Real = x * a.Real;
            c.Imaginary = x * a.Imaginary;
            return c;
        }
        public static Complex operator *(Complex a, double x)
        {
            Complex c = new Complex();
            c.Real = x * a.Real;
            c.Imaginary = x * a.Imaginary;
            return c;
        }
        public static Complex operator /(double x, Complex a)
        {
            Complex xx = new Complex(x, 0.0);
            Complex c = xx / a;
            return c;
        }
        public static Complex operator /(Complex a, double x)
        {
            Complex xx = new Complex(x, 0.0);
            Complex c = a / xx;
            return c;
        }
        public static Complex Conjugate(Complex z)
        {
            Complex c = new Complex();
            c.Real = z.Real;
            c.Imaginary = -z.Imaginary;
            return c;
        }
        public static double Norm(Complex z)
        {
            double x, y, ans, temp;
            x = Math.Abs(z.Real);
            y = Math.Abs(z.Imaginary);
            if (x == 0.0)
                ans = y;
            else if (y == 0.0)
                ans = x;
            else if (x > y)
            {
                temp = y / x;
                ans = x * Math.Sqrt(1.0 + temp * temp);
            }
            else
            {
                temp = x / y;
                ans = y * Math.Sqrt(1.0 + temp * temp);
            }
            return ans;
        }
        public static Complex Sqrt(Complex z)
        {
            Complex c = new Complex(0.0, 0.0);
            double x, y, w, r;
            if ((z.Real == 0.0) && (z.Imaginary == 0.0))
            {
                c.Real = 0.0;
                c.Imaginary = 0.0;
                return c;
            }
            else
            {
                x = Math.Abs(z.Real);
                y = Math.Abs(z.Imaginary);
                if (x >= y)
                {
                    r = y / x;
                    w = Math.Sqrt(x) * Math.Sqrt(0.5 * (1.0 + Math.Sqrt(1.0 + r * r)));
                }
                else
                {
                    r = x / y;
                    w = Math.Sqrt(y) * Math.Sqrt(0.5 * (r + Math.Sqrt(1.0 + r * r)));
                }
                if (z.Real >= 0.0)
                {
                    c.Real = w;
                    c.Imaginary = z.Imaginary / (2.0 * w);
                }
                else
                {
                    c.Imaginary = (z.Imaginary >= 0) ? w : -w;
                    c.Real = z.Imaginary / (2.0 * c.Imaginary);
                }
                return c;
            }
        }

        public static Complex I
        {
            get { Complex i = new Complex(0.0, 1.0); return i; }
        }

        public void Chop()
        {
            if (Math.Abs(this.Real) < 0.0000000001) this.Real = 0.0;
            if (Math.Abs(this.Imaginary) < 0.0000000001) this.Imaginary = 0.0;      
        }

        public override string ToString()
        {
            string res = this.Real.ToString();
            if (this.Imaginary == 0.0) res += "";
            else if (this.Imaginary > 0.0) 
                res += " +I "+this.Imaginary.ToString() + "";
            else
                res += " -I " + (-this.Imaginary).ToString() + "";
            return res;
        } 
    }
}