/*
 * Open Source Physics software is free software as described near the bottom of this code file.
 *
 * For additional information and documentation on Open Source Physics please see:
 * <http://www.opensourcephysics.org/>
 */

/*

    Adapted to OSP by Javier E. Hasbun (2009) with added functionality
    from the original JEP - Java Math Expression Parser 2.24 of Nathan Funk.
  <P>
    @Copyright (c) 2009
    This software is to support the Open Source Physics library
    http://www.opensourcephysics.org under the terms of the GNU General Public
    License (GPL) as published by the Free Software Foundation.

 **/
package org.opensourcephysics.numerics;

/**
 * Class description
 *
*/
public class Complex {
  /** the real component */
  private double re;

  /** the imaginary component */
  private double im;

  //------------------------------------------------------------------------
  // Constructors

  /**
   * Default constructor.
   */
  public Complex() {
    re = 0;
    im = 0;
  }

  /**
   * Constructor from a single double value. The complex number is
   * initialized with the real component equal to the parameter, and
   * the imaginary component equal to zero.
   * @param re_in
   */
  public Complex(double re_in) {
    re = re_in;
    im = 0;
  }

  /**
   * Construct from a Number. This constructor uses the doubleValue()
   * method of the parameter to initialize the real component of the
   * complex number. The imaginary component is initialized to zero.
   * @param re_in
   */
  public Complex(Number re_in) {
    re = re_in.doubleValue();
    im = 0;
  }

  /**
   * Copy constructor
   * @param z
   */
  public Complex(Complex z) {
    re = z.re;
    im = z.im;
  }

  /**
   * Initialize the real and imaginary components to the values given
   * by the parameters.
   * @param re_in
   * @param im_in
   */
  public Complex(double re_in, double im_in) {
    re = re_in;
    im = im_in;
  }

  /**
   * Convert Cartesian to polar
   */
  public Complex polar() {
    double r = StrictMath.sqrt(re*re+im*im);
    double a = StrictMath.atan2(im, re);
    return new Complex(r, a);
  }

  /**
   * Convert polar to Cartesian
   */
  public Complex cartesian() {
    return new Complex(re*StrictMath.cos(im), re*StrictMath.sin(im));
  }

  /**
   * Returns the real component of this object
   */
  public double re() {
    return re;
  }

  /**

   * Returns the imaginary component of this object

   */
  public double im() {
    return im;
  }

  /**
   * Copies the values from the parameter object to this object
   */
  public void set(Complex z) {
    re = z.re;
    im = z.im;
  }

  /**
   * Sets the real and imaginary values of the object.
   */
  public void set(double re_in, double im_in) {
    re = re_in;
    im = im_in;
  }

  /**
   * Sets the real component of the object
   */
  public void setRe(double re_in) {
    re = re_in;
  }

  /**
   * Sets the imaginary component of the object
   */
  public void setIm(double im_in) {
    im = im_in;
  }

  //------------------------------------------------------------------------
  // Various functions

  /**
   * Compares this object with the Complex number given as parameter
   * <pre>b</pre>. The <pre>tolerance</pre> parameter is the radius
   * within which the <pre>b</pre> number must lie for the two
   * complex numbers to be considered equal.
   *
   * @return <pre>true</pre> if the complex number are considered equal,
   * <pre>false</pre> otherwise.
   */
  public boolean equals(Complex b, double tolerance) {
    double temp1 = (re-b.re);
    double temp2 = (im-b.im);
    return(temp1*temp1+temp2*temp2)<=tolerance*tolerance;
  }

  /**
   * Returns the value of this complex number as a string in the format:
   * <pre>(real, imaginary)</pre>.
   */
  public String toString() {
    return "("+re+", "+im+")"; //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
  }

  /**
   * Returns the absolute value of the complex number.
   */
  public double abs() {
    double absRe = Math.abs(re);
    double absIm = Math.abs(im);
    if((absRe==0)&&(absIm==0)) {
      return 0;
    } else if(absRe>absIm) {
      double temp = absIm/absRe;
      return absRe*Math.sqrt(1+temp*temp);
    } else {
      double temp = absRe/absIm;
      return absIm*Math.sqrt(1+temp*temp);
    }
  }

  /**

   * Returns the square of the absolute value (re*re+im*im).

   */
  public double abs2() {
    return re*re+im*im;
  }

  /**
   * Returns the magnitude of the complex number
   */
  public double mag() {
    return StrictMath.sqrt(re*re+im*im);
  }

  /**
   * Returns the argument of this complex number (Math.atan2(re,im))
   */
  public double arg() {
    return Math.atan2(im, re);
  }

  /** Convert text representation to a Complex.
   *  input format  (real_double,imaginary_double)
   */
  public static Complex parseComplex(String s) {
    int from = s.indexOf('(');
    if(from==-1) {
      return null;
    }
    int to = s.indexOf(',', from);
    double x = Double.parseDouble(s.substring(from+1, to));
    from = to;
    to = s.indexOf(')', from);
    double y = Double.parseDouble(s.substring(from+1, to));
    return new Complex(x, y);
  }

  /**
   * Returns the sum of two complex numbers
   */
  public Complex add(Complex z) {
    return new Complex(re+z.re, im+z.im);
  }

  /**
   * Returns the sum of a complex number and a double
   */
  public Complex add(double d) {
    return new Complex(re+d, im);
  }

  /**
   * Returns the subtraction of complex z from a complex number
   */
  public Complex subtract(Complex z) {
    return new Complex(re-z.re, im-z.im);
  }

  /**
   * Subtracts the double d from the complex number
   * Returns the subtraction of complex z from a complex number
   */
  public Complex subtract(double d) {
    return new Complex(re-d, im);
  }

  /**
   * Returns the negative value of this complex number.
   */
  public Complex neg() {
    return new Complex(-re, -im);
  }

  /**
   * Multiplies the complex number with a double value.
   * @return The result of the multiplication

   */
  public Complex mul(double b) {
    return new Complex(re*b, im*b);
  }

  /**
   * Multiplies the complex number with another complex value.
   * @return The result of the multiplication
   */
  public Complex mul(Complex b) {
    return new Complex(re*b.re-im*b.im, im*b.re+re*b.im);
  }

  /**
   * Returns the result of dividing this complex number by the parameter.
   */
  public Complex div(Complex b) {
    // Adapted from Numerical Recipes in C - The Art of Scientific Computing
    // ISBN 0-521-43108-5
    double resRe, resIm;
    double r, den;
    if(Math.abs(b.re)>=Math.abs(b.im)) {
      r = b.im/b.re;
      den = b.re+r*b.im;
      resRe = (re+r*im)/den;
      resIm = (im-r*re)/den;
    } else {
      r = b.re/b.im;
      den = b.im+r*b.re;
      resRe = (re*r+im)/den;
      resIm = (im*r-re)/den;
    }
    return new Complex(resRe, resIm);
  }

  /**

   * Devide the complex number by a double value.

   * @return The result of the division

   */
  public Complex div(double d) {
    return new Complex(re/d, im/d);
  }

  /**

   * Invert the complex number.

   * @return The result of the inversion

   */

  /**  */
  public Complex invert() {
    double r = re*re+im*im;
    return new Complex(re/r, -im/r);
  }

  /**
   * Conjugats the complex number
   */
  public Complex conjugate() {
    return new Complex(re, -im);
  }

  /**
   * Returns the value of this complex number raised to the power
   * of a real component (in double precision).<p>
   * This method considers special cases where a simpler algorithm
   * would return "ugly" results.<br>
   * For example when the expression (-1e40)^0.5 is evaluated without
   * considering the special case, the argument of the base is the
   * double number closest to pi. When sin and cos are used for the
   * final evaluation of the result, the slight difference of the
   * argument from pi causes a non-zero value for the real component
   * of the result. Because the value of the base is so high, the error
   * is magnified.Although the error is normal for floating
   * point calculations, the consideration of commonly occurring special
   * cases improves the accuracy and esthetics of the results.<p>
   * If you know a more elegant way to solve this problem, please let
   * me know at nathanfunk@hotmail.com .

   */
  public Complex power(double exponent) {
    // z^exp = abs(z)^exp * (cos(exp*arg(z)) + i*sin(exp*arg(z)))
    double scalar = Math.pow(abs(), exponent);
    boolean specialCase = false;
    int factor = 0;
    // consider special cases to avoid floating point errors
    // for power expressions such as (-1e20)^2
    if((im==0)&&(re<0)) {
      specialCase = true;
      factor = 2;
    }
    if((re==0)&&(im>0)) {
      specialCase = true;
      factor = 1;
    }
    if((re==0)&&(im<0)) {
      specialCase = true;
      factor = -1;
    }
    if(specialCase&&(factor*exponent==(int) (factor*exponent))) {
      short[] cSin = {0, 1, 0, -1}; // sin of 0, pi/2, pi, 3pi/2
      short[] cCos = {1, 0, -1, 0}; // cos of 0, pi/2, pi, 3pi/2
      int x = ((int) (factor*exponent))%4;
      if(x<0) {
        x = 4+x;
      }
      return new Complex(scalar*cCos[x], scalar*cSin[x]);
    }
    double temp = exponent*arg();
    return new Complex(scalar*Math.cos(temp), scalar*Math.sin(temp));
  }

  /**
   * Returns the value of this complex number raised to the power of
   * a complex exponent
   */
  public Complex power(Complex exponent) {
    if(exponent.im==0) {
      return power(exponent.re);
    }
    double temp1Re = Math.log(abs());
    double temp1Im = arg();
    double temp2Re = (temp1Re*exponent.re)-(temp1Im*exponent.im);
    double temp2Im = (temp1Re*exponent.im)+(temp1Im*exponent.re);
    double scalar = Math.exp(temp2Re);
    return new Complex(scalar*Math.cos(temp2Im), scalar*Math.sin(temp2Im));
  }

  /**
   * Returns e to the power of the complex number
   */
  public Complex exp() {
    double exp_x = StrictMath.exp(re);
    return new Complex(exp_x*StrictMath.cos(im), exp_x*StrictMath.sin(im));
  }

  /**
   * Returns the logarithm of this complex number.
   */
  public Complex log() {
    return new Complex(Math.log(abs()), arg());
  }

  /**
   * Calculates the square root of this object.
   */
  public Complex sqrt() {
    Complex c;
    double absRe, absIm, w, r;
    if((re==0)&&(im==0)) {
      c = new Complex(0, 0);
    } else {
      absRe = Math.abs(re);
      absIm = Math.abs(im);
      if(absRe>=absIm) {
        r = absIm/absRe;
        w = Math.sqrt(absRe)*Math.sqrt(0.5*(1.0+Math.sqrt(1.0+r*r)));
      } else {
        r = absRe/absIm;
        w = Math.sqrt(absIm)*Math.sqrt(0.5*(r+Math.sqrt(1.0+r*r)));
      }
      if(re>=0) {
        c = new Complex(w, im/(2.0*w));
      } else {
        if(im<0) {
          w = -w;
        }
        c = new Complex(im/(2.0*w), w);
      }
    }
    return c;
  }

  //------------------------------------------------------------------------
  // Trigonometric functions

  /**
   * Returns the sine of this complex number.
   */
  public Complex sin() {
    double izRe, izIm;
    double temp1Re, temp1Im;
    double temp2Re, temp2Im;
    double scalar;
    //  sin(z)  =  ( exp(i*z) - exp(-i*z) ) / (2*i)
    izRe = -im;
    izIm = re;
    // first exp
    scalar = Math.exp(izRe);
    temp1Re = scalar*Math.cos(izIm);
    temp1Im = scalar*Math.sin(izIm);
    // second exp
    scalar = Math.exp(-izRe);
    temp2Re = scalar*Math.cos(-izIm);
    temp2Im = scalar*Math.sin(-izIm);
    temp1Re -= temp2Re;
    temp1Im -= temp2Im;
    return new Complex(0.5*temp1Im, -0.5*temp1Re);
  }

  /**

   * Returns the cosine of this complex number.

   */
  public Complex cos() {
    double izRe, izIm;
    double temp1Re, temp1Im;
    double temp2Re, temp2Im;
    double scalar;
    //  cos(z)  =  ( exp(i*z) + exp(-i*z) ) / 2
    izRe = -im;
    izIm = re;
    // first exp
    scalar = Math.exp(izRe);
    temp1Re = scalar*Math.cos(izIm);
    temp1Im = scalar*Math.sin(izIm);
    // second exp
    scalar = Math.exp(-izRe);
    temp2Re = scalar*Math.cos(-izIm);
    temp2Im = scalar*Math.sin(-izIm);
    temp1Re += temp2Re;
    temp1Im += temp2Im;
    return new Complex(0.5*temp1Re, 0.5*temp1Im);
  }

  /**
   * Returns the tangent of this complex number.
   */
  public Complex tan() {
    // tan(z) = sin(z)/cos(z)
    double izRe, izIm;
    double temp1Re, temp1Im;
    double temp2Re, temp2Im;
    double scalar;
    Complex sinResult, cosResult;
    //  sin(z)  =  ( exp(i*z) - exp(-i*z) ) / (2*i)
    izRe = -im;
    izIm = re;
    // first exp
    scalar = Math.exp(izRe);
    temp1Re = scalar*Math.cos(izIm);
    temp1Im = scalar*Math.sin(izIm);
    // second exp
    scalar = Math.exp(-izRe);
    temp2Re = scalar*Math.cos(-izIm);
    temp2Im = scalar*Math.sin(-izIm);
    temp1Re -= temp2Re;
    temp1Im -= temp2Im;
    sinResult = new Complex(0.5*temp1Re, 0.5*temp1Im);
    //  cos(z)  =  ( exp(i*z) + exp(-i*z) ) / 2
    izRe = -im;
    izIm = re;
    // first exp
    scalar = Math.exp(izRe);
    temp1Re = scalar*Math.cos(izIm);
    temp1Im = scalar*Math.sin(izIm);
    // second exp
    scalar = Math.exp(-izRe);
    temp2Re = scalar*Math.cos(-izIm);
    temp2Im = scalar*Math.sin(-izIm);
    temp1Re += temp2Re;
    temp1Im += temp2Im;
    cosResult = new Complex(0.5*temp1Re, 0.5*temp1Im);
    return sinResult.div(cosResult);
  }

  //------------------------------------------------------------------------
  // Inverse trigonometric functions
  public Complex asin() {
    Complex result;
    double tempRe, tempIm;
    //  asin(z)  =  -i * log(i*z + sqrt(1 - z*z))
    tempRe = 1.0-((re*re)-(im*im));
    tempIm = 0.0-((re*im)+(im*re));
    result = new Complex(tempRe, tempIm);
    result = result.sqrt();
    result.re += -im;
    result.im += re;
    tempRe = Math.log(result.abs());
    tempIm = result.arg();
    result.re = tempIm;
    result.im = -tempRe;
    return result;
  }

  public Complex acos() {
    Complex result;
    double tempRe, tempIm;
    //  acos(z)  =  -i * log( z + i * sqrt(1 - z*z) )
    tempRe = 1.0-((re*re)-(im*im));
    tempIm = 0.0-((re*im)+(im*re));
    result = new Complex(tempRe, tempIm);
    result = result.sqrt();
    tempRe = -result.im;
    tempIm = result.re;
    result.re = re+tempRe;
    result.im = im+tempIm;
    tempRe = Math.log(result.abs());
    tempIm = result.arg();
    result.re = tempIm;
    result.im = -tempRe;
    return result;
  }

  public Complex atan() {
    // atan(z) = -i/2 * log((i-z)/(i+z))
    double tempRe, tempIm;
    Complex result = new Complex(-re, 1.0-im);
    tempRe = re;
    tempIm = 1.0+im;
    result = result.div(new Complex(tempRe, tempIm));
    tempRe = Math.log(result.abs());
    tempIm = result.arg();
    result.re = 0.5*tempIm;
    result.im = -0.5*tempRe;
    return result;
  }

  //------------------------------------------------------------------------
  // Hyperbolic trigonometric functions
  public Complex sinh() {
    double scalar;
    double temp1Re, temp1Im;
    double temp2Re, temp2Im;
    //  sinh(z)  =  ( exp(z) - exp(-z) ) / 2
    // first exp
    scalar = Math.exp(re);
    temp1Re = scalar*Math.cos(im);
    temp1Im = scalar*Math.sin(im);
    // second exp
    scalar = Math.exp(-re);
    temp2Re = scalar*Math.cos(-im);
    temp2Im = scalar*Math.sin(-im);
    temp1Re -= temp2Re;
    temp1Im -= temp2Im;
    return new Complex(0.5*temp1Re, 0.5*temp1Im);
  }

  public Complex cosh() {
    double scalar;
    double temp1Re, temp1Im;
    double temp2Re, temp2Im;
    //  cosh(z)  =  ( exp(z) + exp(-z) ) / 2
    // first exp
    scalar = Math.exp(re);
    temp1Re = scalar*Math.cos(im);
    temp1Im = scalar*Math.sin(im);
    // second exp
    scalar = Math.exp(-re);
    temp2Re = scalar*Math.cos(-im);
    temp2Im = scalar*Math.sin(-im);
    temp1Re += temp2Re;
    temp1Im += temp2Im;
    return new Complex(0.5*temp1Re, 0.5*temp1Im);
  }

  public Complex tanh() {
    double scalar;
    double temp1Re, temp1Im;
    double temp2Re, temp2Im;
    Complex sinRes, cosRes;
    //  tanh(z)  =  sinh(z) / cosh(z)
    scalar = Math.exp(re);
    temp1Re = scalar*Math.cos(im);
    temp1Im = scalar*Math.sin(im);
    scalar = Math.exp(-re);
    temp2Re = scalar*Math.cos(-im);
    temp2Im = scalar*Math.sin(-im);
    temp1Re -= temp2Re;
    temp1Im -= temp2Im;
    sinRes = new Complex(0.5*temp1Re, 0.5*temp1Im);
    scalar = Math.exp(re);
    temp1Re = scalar*Math.cos(im);
    temp1Im = scalar*Math.sin(im);
    scalar = Math.exp(-re);
    temp2Re = scalar*Math.cos(-im);
    temp2Im = scalar*Math.sin(-im);
    temp1Re += temp2Re;
    temp1Im += temp2Im;
    cosRes = new Complex(0.5*temp1Re, 0.5*temp1Im);
    return sinRes.div(cosRes);
  }

  //------------------------------------------------------------------------
  // Inverse hyperbolic trigonometric functions
  public Complex asinh() {
    Complex result;
    //  asinh(z)  =  log(z + sqrt(z*z + 1))
    result = new Complex(((re*re)-(im*im))+1, (re*im)+(im*re));
    result = result.sqrt();
    result.re += re;
    result.im += im;
    double temp = result.arg();
    result.re = Math.log(result.abs());
    result.im = temp;
    return result;
  }

  public Complex acosh() {
    Complex result;
    //  acosh(z)  =  log(z + sqrt(z*z - 1))
    result = new Complex(((re*re)-(im*im))-1, (re*im)+(im*re));
    result = result.sqrt();
    result.re += re;
    result.im += im;
    double temp = result.arg();
    result.re = Math.log(result.abs());
    result.im = temp;
    return result;
  }

  public Complex atanh() {
    //  atanh(z)  =  1/2 * log( (1+z)/(1-z) )
    double tempRe, tempIm;
    Complex result = new Complex(1.0+re, im);
    tempRe = 1.0-re;
    tempIm = -im;
    result = result.div(new Complex(tempRe, tempIm));
    tempRe = Math.log(result.abs());
    tempIm = result.arg();
    result.re = 0.5*tempRe;
    result.im = 0.5*tempIm;
    return result;
  }

  /**
   * Returns <tt>true</tt> if either the real or imaginary component of this
   * <tt>Complex</tt> is an infinite value.
   * <p>
   * @return                  <tt>true</tt> if either component of the <tt>Complex</tt> object is infinite; <tt>false</tt>, otherwise.
   * <p>
   */
  public boolean isInfinite() {
    return(Double.isInfinite(re)||Double.isInfinite(im));
  } // end isInfinite()

  /**
   * Returns <tt>true</tt> if either the real or imaginary component of this
   * <tt>Complex</tt> is a Not-a-Number (<tt>NaN</tt>) value.
   * <p>
   * @return                  <tt>true</tt> if either component of the <tt>Complex</tt> object is <tt>NaN</tt>; <tt>false</tt>, otherwise.
   * <p>
   */
  public boolean isNaN() {
    return(Double.isNaN(re)||Double.isNaN(im));
  } // end isNaN()

}

/*
 * Open Source Physics software is free software; you can redistribute
 * it and/or modify it under the terms of the GNU General Public License (GPL) as
 * published by the Free Software Foundation; either version 2 of the License,
 * or(at your option) any later version.

 * Code that uses any portion of the code in the org.opensourcephysics package
 * or any subpackage (subdirectory) of this package must must also be be released
 * under the GNU GPL license.
 *
 * This software is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
 * or view the license online at http://www.gnu.org/copyleft/gpl.html
 *
 * Copyright (c) 2007  The Open Source Physics project
 *                     http://www.opensourcephysics.org
 */
