/* The functions in this file allow for arithmetic operations on 
   complex numbers */

/* a complex number is defined in "complex.h" as a structure 
   of two double values : the real (.re) and imaginary (.im) parts */


#include <complex.h>


/* the function Complex() returns a complex number with specified real and 
   imaginary parts */

dcomplex Complex(double real, double imag)
{
     dcomplex c;
     c.re = real;
     c.im = imag;
     return c;
}


/* the function Cecho() echoes a complex number to the screen */

void Cecho(dcomplex z)
{
     printf("(%.4e,%.4e )", z.re, z.im );
}


/* the function Cre() returns the real part of a complex number */

double Cre(dcomplex z)
{
     double a;
     a = z.re;
     return a;
}


/* the function Cim() returns the imaginary part of a complex number */

double Cim(dcomplex z)
{
     double a;
     a = z.im;
     return a;
}


/* the function Conjg() returns the complex conjugate of a complex number */

dcomplex Conjg(dcomplex z)
{
     dcomplex c;
     c.re = z.re;
     c.im = - z.im;
     return c;
}


/* the function Cabs() returns the modulus of a complex number*/

double Cabs(dcomplex z)
{
     double x,y,ans;
     x = fabs(z.re);
     y = fabs(z.im);
     if (x == 0.0)
	ans = y;
     else
	if (y == 0.0)
	   ans = x;
	else
	   ans = sqrt(x * x + y * y);
     return ans;
}


/* the function Cinv() returns the inverse of a complex number */

dcomplex Cinv(dcomplex z)
{
     dcomplex c;
     c.re = z.re / (z.re * z.re + z.im * z.im);
     c.im = - z.im / (z.re * z.re + z.im * z.im);
     return c;
}


/* the function Carg() returns the argument of a complex number in radians 
   in the range 0 <= arg < 2 PI */

double Carg(dcomplex z)
{
     double arg;
     if (z.im >= 0)
	arg =  acos(z.re / Cabs(z) );
     else
	arg = 2 * PI - acos(z.re / Cabs(z) );
     return arg;
}


/* the function Cadd() adds two complex numbers */

dcomplex Cadd(dcomplex a, dcomplex b)
{
     dcomplex c;
     c.re = a.re + b.re;
     c.im = a.im + b.im;
     return c;
}

/* the function Csub() subtracts two complex numbers */


dcomplex Csub(dcomplex a, dcomplex b)
{
     dcomplex c;
     c.re = a.re - b.re;
     c.im = a.im - b.im;
     return c;
}


/* the function Cmul() multiplies two complex numbers */

dcomplex Cmul(dcomplex a, dcomplex b)
{
     dcomplex c;
     c.re = a.re * b.re - a.im * b.im;
     c.im = a.re * b.im + a.im * b.re;
     return c;
}


/* the function RCmul() returns the complex product of a real number
   and a complex number */

dcomplex CRmul(double x, dcomplex a)
{
     dcomplex c;
     c.re = x * a.re;
     c.im = x * a.im;
     return c;
}


/* the function Cpow() returns the complex number z raised to a real power */

dcomplex Cpow(dcomplex z, double p)
{
     dcomplex c;
     c.re = pow(Cabs(z), p) * cos( p * Carg(z));
     c.im = pow(Cabs(z), p) * sin( p * Carg(z));
     
     return c;
}


/* the function Cexp() returns the exponential function 
   with a complex argument */
   
dcomplex Cexp(dcomplex z)
{
     double r, i;
     
     r = z.re;
     i = z.im;
     
     return CRmul( exp(r), Complex(cos(i), sin(i)) );
}

/* the function Clog() returns the logarithm function 
   with a complex argument */
   
dcomplex Clog(dcomplex z)
{
     double mod, phase;
     
     mod = Cabs(z);
     phase = Carg2(z);
     
     return Cadd( Complex(log(mod),0.), Complex(0., phase) );
}

/* the function Carg2() returns the argument of a complex number in radians 
   in the range -PI <= arg < PI */

double Carg2(dcomplex z)
{
     double arg;
	
	if (z.im >= 0)
	arg =  acos(z.re / Cabs(z) );
      else
	arg = - acos(z.re / Cabs(z) );

     return arg;
}

dcomplex Csin(dcomplex z)
{
        return Cmul(Complex(0.,-.5),Csub(Cexp(Cmul(Complex(0.,1.),z)),
                Cexp(Cmul(Complex(0.,-1.),z))));
}