/******************************************************************

			COMPLEX

*******************************************************************

File contains routines to manipulate complex numbers
*/

#include <math.h>
#include "complex.h"
#include "pi.h"

/********************************************************************

			cmpy()

*********************************************************************

Function multiplies two complex numbers to give complex result
*/

COMPLEX cmpy(a,b)
COMPLEX a,b;
{
	static COMPLEX c;

	c.real = a.real * b.real - a.imag * b.imag;
	c.imag = a.real * b.imag + a.imag * b.real;

	return(c);
}

/********************************************************************

			cadd()

**********************************************************************

Function adds two complex numbers to give a complex result
*/

COMPLEX cadd(a,b)
COMPLEX a,b;
{
	static COMPLEX c;

	c.real = a.real + b.real;
	c.imag = a.imag + b.imag;

	return(c);
}

/**********************************************************************

			csub()

************************************************************************

Function subtracts two complex numbers to give a complex result
*/

COMPLEX csub(a,b)
COMPLEX a,b;
{
	static COMPLEX c;

	c.real = a.real - b.real;
	c.imag = a.imag - b.imag;

	return(c);
}

/***********************************************************************

			conj()

***********************************************************************

Function returns the complex conjugate of a complex argument
*/

COMPLEX conj(a)
COMPLEX a;
{
	static COMPLEX b;

	b.real = a.real;
	b.imag = - a.imag;

	return(b);
}

/******************************************************************

			cdiv()

*******************************************************************

Function divides two complex numbers
*/

COMPLEX cdiv(a,b)
COMPLEX a,b;

{
	static COMPLEX c;
	COMPLEX an,bn;	/* normalized versions of a and b */

	double magb; /* magnitude of b */

	magb = cmag(b);

/* now normalize a and b by magb */

	an.real = a.real/magb;
	an.imag = a.imag/magb;
	bn.real = b.real/magb;
	bn.imag = b.imag/magb;

/* finally, a/b is an * conj(bn) */

	c = cmpy(an,conj(bn));

	return(c);
}

/***********************************************************************

			csqrt()

************************************************************************

Function returns the square root of a complex number
*/

COMPLEX csqrt(a)
COMPLEX a;

{
	static COMPLEX b;
	double mag,phase,x,y;

/* convert to polar coordinates */

/* first the magnitude */

	mag = cmag(a);

/* calculate the phase */

	phase = atan2(a.imag,a.real);

/* take the square root */

	mag = sqrt(mag);
	phase /= 2.0;

/* convert back to rectangular coordinates */

	b.real = mag * cos(phase);
	b.imag = mag * sin(phase);

	return(b);
}

/****************************************************************

			cmag()

*****************************************************************

Function calculates the magnitude of a complex number in such a
way that it doesnt overflow if large */

double cmag(a)
COMPLEX a;

{
	double maga,m1,m2; /* magnitude of a */


/* if a.real and a.imag both small do it the easy way */

if ( fabs(a.real) < 1.e10 && fabs(a.imag) < 1.e10 )
	maga = sqrt( a.real*a.real + a.imag*a.imag );

else if ( fabs(a.real) < fabs(a.imag) ) {
	maga = a.real/a.imag;
	maga = fabs(a.imag) * sqrt( 1.0 + maga*maga );
	}

else {
	maga = a.imag/a.real;
	maga = fabs(a.real) * sqrt( 1.0 + maga*maga );
	}

	return(maga);
}

/****************************************************************

			cexp()

*****************************************************************

Function returns a complex exponential

Single argument is the angle in radians  */

COMPLEX cexp(omega)
double omega;

{
	static COMPLEX c;

	c.real = cos(omega);
	c.imag = sin(omega);

	return(c);
}

/**********************************************************************

			cpl()

**********************************************************************

Function turns two real numbers into a complex number

Two arguments are the real and imaginary parts */

COMPLEX cpl(real,imag)
double real, imag;

{
	static COMPLEX c;

	c.real = real;
	c.imag = imag;

	return(c);
}
