/***************************************************************************
                          complex.h  -  description
                             -------------------
    begin                : Fri Jan 17 19:41:58 CST 2003
    copyright            : (C) 2003 by Daniel Eric Smith
    email                : 
 ***************************************************************************/

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

/*
fix asin and acos
*/

#include <math.h>
#include <iomanip.h>


class complex{
  double real,imag;
 public:
  complex() {real=0; imag=0;}
  complex(double x) {real=x; imag=0;}
  complex(double r, double i) {real=r; imag=i;}  //parameterized construtor
  void set(double r, double i);                  //initializes z
  double Re();                                //returns the real part of z
  double Im();                                //returns the imagaginary part of z
  double mod();                                  //returns the modulus of z
  double arg();                                  //returns the argument of z
/*
CONJUGATE AND ARITMETIC
*/
  complex operator~();                           //complex conjugate operator ~(a+bi)=>a-bi
  friend complex operator+(complex z, complex w); //complex addition z+w
  friend complex operator+(complex z, double r);
  friend complex operator+(double r, complex z);
  friend complex operator-(complex z, complex w);                  //complex subtraction z-w
  friend complex operator-(complex z, double r);
  friend complex operator-(double r, complex z);
  friend complex operator*(complex z, complex w);                  //complex multiplication
  friend complex operator*(complex z, double r);
  friend complex operator*(double r, complex z);
  friend complex operator/(complex z, complex w);                  //complex division
  friend complex operator/(complex z, double r);
  friend complex operator/(double r, complex z);
  complex operator=(complex w);                  //complex equality
/*
OUTPUT
*/
  friend istream& operator>>(istream& s, complex& z);
  friend ostream& operator<<(ostream& s, complex z);
/*
EXPONENTIAL
*/  
  friend complex exp(complex z);
/*
TRIG
*/
  friend complex sin(complex z);
  friend complex cos(complex z);
  friend complex tan(complex z);
  friend complex cot(complex z);
  friend complex sec(complex z);
  friend complex csc(complex z);
  friend complex asin(complex z);
  friend complex acos(complex z);
  friend complex atan(complex z);
/*
HYPERBOLIC TRIG
*/
  friend complex sinh(complex z);
  friend complex cosh(complex z);
  friend complex tanh(complex z);
  friend complex coth(complex z);
  friend complex sech(complex z);
  friend complex csch(complex z);
  friend complex asinh(complex z);
  friend complex acosh(complex z);
  friend complex atanh(complex z);
/*
LOG
*/
  friend complex log(complex z);
/*
EXPONENT
*/
 friend complex pow(complex z, complex w);
 friend complex pow(complex z, double x);
 friend complex pow(double x, complex w);
};

/*****************************
*                            *
* I defines imagaginary unit *
*                            *
*****************************/
const complex I(0,1);





void complex::set(double r, double i)
{
  real=r;
  imag=i;
}

double complex::Re()
{
  return real;
}

double complex::Im()
{
  return imag;
}

double complex::mod()
{
  return sqrt(pow(real,2)+pow(imag,2));
}

double complex::arg()
{
  return atan2(imag,real);
}




/**************************
* CONJUGATE AND ARITMETIC *
**************************/

/****************
*               *
* ~z  CONJUGATE *
*               *
****************/
complex complex::operator~()
{
  complex temp;

  temp=*this;
  temp.imag=(-1)*imag;

  return temp;
}
/****************
*               *
* z+w  ADDITION * 
*               *
****************/
complex operator+(complex z, complex w)
{
	complex zsum;
	
	zsum.real=z.real+w.real;
	zsum.imag=z.imag+w.imag;

	return zsum;
}
complex operator+(complex z, double r)
{
  complex zsum,t(r);

  zsum.real=z.real+t.real;
  zsum.imag=z.imag+t.imag;

  return zsum;
}
complex operator+(double r, complex z)
{
  return z+r;
}
/******************
*                 *
* z-w SUBTRACTION *
*                 *
******************/
complex operator-(complex z, complex w)
{
  complex zdiff;

  zdiff.real=z.real-w.real;
  zdiff.imag=z.imag-w.imag;

  return zdiff;
}
complex operator-(complex z, double r)
{
  complex zdiff,t(r);

  zdiff.real=z.real-t.real;
  zdiff.imag=z.imag-t.imag;

  return zdiff;
}
complex operator-(double r, complex z)
{
  complex zdiff,t(r);

  zdiff.real=t.real-z.real;
  zdiff.imag=t.imag-z.imag;

  return zdiff;  
}
/*********************
*                    *
* z*w MULTIPLICATION *
*                    *
*********************/
complex operator*(complex z, complex w)
{
  complex zprod;

  zprod.real=z.real*w.real-z.imag*w.imag;
  zprod.imag=z.real*w.imag+z.imag*w.real;

  return zprod;
}
complex operator*(complex z, double r)
{
	complex zprod;

	zprod.real=z.real*r;
	zprod.imag=z.imag*r;

	return zprod;
}
complex operator*(double r, complex z)
{
	return z*r;
}
/***************
*              *
* z/w DIVISION *
*              *
***************/
complex operator/(complex z, complex w)
{
  complex zquot,num,denom,wc;
  wc=~w;
  
  num = z * wc;
  denom = w * wc;
  
  zquot.real=num.real/denom.real;
  zquot.imag=num.imag/denom.real;

  return zquot;
}
complex operator/(complex z, double r)
{
	complex zquot;

	zquot.real=z.real/r;
	zquot.imag=z.imag/r;

	return zquot;
}
complex operator/(double r, complex z)
{
	complex zquot, zc, denom,num;

	zc=~z;
	num=r*zc;
	denom=z*zc;

	zquot.real=num.real/denom.real;
	zquot.imag=num.imag/denom.real;

  return zquot;
}
/***************
*              *
* z=w EQUALITY *
*              *
***************/
complex complex::operator=(complex w)
{
  real=w.real;
  imag=w.imag;

  return *this;
}




/***************
* INPUT/OUTPUT *
***************/
istream& operator>>(istream& s, complex& z)
{
  s>>z.real>>z.imag;
  return s;
}

ostream& operator<<(ostream& s, complex z)
{
 double imag2;

  if(z.imag>=0)
    {  
      s<<z.real<<"+"<<z.imag<<"i";
    }
  else
    {
      imag2=(-1)*z.imag;
      s <<z.real<<"-"<<imag2<<"i";
    }
  return s;
}

/*********************
*                    *
* exp(z) EXPONENTIAL *
*                    *
*********************/
complex exp(complex z)
{
	complex e;

	e.real=exp(z.real)*cos(z.imag);
	e.imag=exp(z.real)*sin(z.imag);

	return e;
}


/*************************
* TRIGONMETRIC FUNCTIONS *
*************************/

/*********
*        *
* sin(x) *
*        *
*********/
complex sin(complex z)
{
	return ((exp(z*I)-exp((-1)*z*I))/(2*I));
}

/*********
*        *
* cos(x) *
*        *
*********/
complex cos(complex z)
{
	return ((exp(z*I)+exp((-1)*z*I))/2);
}	

/*********
*        *
* tan(x) *
*        *
*********/
complex tan(complex z)
{
	return sin(z)/cos(z);
}

/*********
*        *
* cot(x) *
*        *
*********/
complex cot(complex z)
{
	return cos(z)/sin(z);
}	

/*********
*        *
* sec(x) *
*        *
*********/
complex sec(complex z)
{
	return 1/cos(z);
}

/*********
*        *
* csc(x) *
*        *
*********/
complex csc(complex z)
{
	return 1/sin(z);
}
/**********
*         *
* asin(x) *
*         *
**********/
complex asin(complex z)
{
  complex t;

  return (-1)*I*log(I*z+pow(1-z*z,.5));
}

/**********
*         *
* acos(x) *
*         *
**********/
complex acos(complex z)
{
	return (-1)*I*log(z+I*pow(1-z*z,.5));
}	

/**********
*         *
* atan(x) *
*         *
**********/
complex atan(complex z)
{
	return I/2*log((I+z)/(I-z));
}



/************************************
* HYPERBOLIC TRIGONMETRIC FUNCTIONS *
************************************/

/**********
*         *
* sinh(x) *
*         *
**********/
complex sinh(complex z)
{
	return ((exp(z)-exp((-1)*z))/2);
}

/**********
*         *
* cosh(x) *
*         *
**********/
complex cosh(complex z)
{
	return ((exp(z)+exp((-1)*z))/2);
}	

/**********
*         *
* tanh(x) *
*         *
**********/
complex tanh(complex z)
{
	return sinh(z)/cosh(z);
}

/**********
*         *
* coth(x) *
*         *
**********/
complex coth(complex z)
{
	return cosh(z)/sinh(z);
}	

/**********
*         *
* sech(x) *
*         *
**********/
complex sech(complex z)
{
	return 1/cosh(z);
}

/**********
*         *
* csch(x) *
*         *
**********/
complex csch(complex z)
{
	return 1/sinh(z);
}
/**********
*         *
* asinh(x) *
*         *
**********/
complex asinh(complex z)
{
	return log(z+pow(z*z+1,.5));
}

/**********
*         *
* acosh(x) *
*         *
**********/
complex acosh(complex z)
{
	return log(z+pow(z*z-1,.5));
}	

/**********
*         *
* atanh(x) *
*         *
**********/
complex atanh(complex z)
{
/*  complex t,t1,t2;

  t1=1+z;
  t2=1-z;
  t=t1/t2;
  t=log(t);
  t=.5*t;

  return t;*/
  return .5*log((1+z)/(1-z));
}

/***********************
*                      *
* LOGARITHMIC FUNCTION *
*                      *
***********************/
complex log(complex z)
{
  double m,a;
  complex l;

  m=z.mod();
  a=z.arg();
  
  l.real=log(m);
  l.imag=a;

  return l;
}

/***********
*          *
* EXPONENT *
*          *
***********/
complex pow(complex z, complex w)
{
  return exp(w*log(z));
}
complex pow(complex z, double x)
{
  return exp(x*log(z));
}
complex pow(double x, complex w)
{
  return exp(w*log(x));
}