/*
 * METRIC --- Mode expansion modeling in integrated optics / photonics
 * http://metric.computational-photonics.eu/ 
 */

/* 
 * complex.cpp
 * Handling complex numbers 
 */

#include<stdio.h>
#include<stdlib.h>
#include<cmath>
#include<complex>
#include"inout.h"
#include"complex.h"

/* error message */
void complexerror(const char *m)
{
	fprintf(stderr, "\ncomplex: %s.\n", m);
	exit(1);
}

/* argument, range [0, 2 PI[ */
double Complex::arg() const
{	
	if(im == 0.0)
	{
		if(re >= 0.0) return 0.0;
		return M_PI;
	}
	if(re == 0.0)
	{
		if(im >= 0.0) return M_PI_2;
		return 1.5*M_PI;
	}
	if(im > 0.0 && re > 0.0) return atan(im/re);
	if(im > 0.0 && re < 0.0) return M_PI-atan(fabs(im/re));
	if(im < 0.0 && re < 0.0) return M_PI+atan(fabs(im/re));
	if(im < 0.0 && re > 0.0) return 2.0*M_PI-atan(fabs(im/re));
	return 0.0;
}

/* argument, range ]-PI, PI] */
double Complex::args() const
{	
	if(im == 0.0)
	{
		if(re >= 0.0) return 0.0;
		return M_PI;
	}
	if(re == 0.0)
	{
		if(im >= 0.0) return M_PI_2;
		return -M_PI_2;
	}
	if(im > 0.0 && re > 0.0) return atan(im/re);
	if(im > 0.0 && re < 0.0) return M_PI-atan(fabs(im/re));
	if(im < 0.0 && re < 0.0) return -M_PI+atan(fabs(im/re));
	if(im < 0.0 && re > 0.0) return -atan(fabs(im/re));
	return 0.0;
}

/* output to FILE dat */
void Complex::write(FILE *dat) const
{
	if(dat == stderr || dat == stdout)
	{
		nice(dat);
		fprintf(dat, "\n");
	}
	else
	{
		fputdouble(re, dat);
		fputdouble(im, dat);
	}
	return;
}
void Complex::nice(FILE *dat) const
{
	if(isnotOK())
	{
		fprintf(dat, "NANorINF");
		return;
	}
	if(re == 0.0 && im == 0.0) 
	{
		fprintf(dat, "%g", 0.0);
		return;
	}
	if(fabs(re) > 0.0) fprintf(dat, "%g", re);
	if(im == 0.0) return;
	if(im > 0.0) fprintf(dat, "+i%g", im);
	else         fprintf(dat, "-i%g", -im);
	return;
}


/* input from FILE dat */
void Complex::read(FILE *dat)
{
	if(dat == stderr || dat == stdout)
	{
		complexerror("Complex: read");
	}
	else
	{
		re = fgetdouble(dat);
		im = fgetdouble(dat);
	}
	return;
}

/* ------------------------------------------------------------------------- */

/* trigonometric functions */
Complex sin(const Complex& z)
{
//	return Complex(sin(z.re)*cosh(z.im), cos(z.re)*sinh(z.im));
	std::complex<double> tz(z.re, z.im);
	std::complex<double> s = sin(tz);
	return Complex(real(s), imag(s));
}
Complex cos(const Complex& z)
{
//	return Complex(cos(z.re)*cosh(z.im), -sin(z.re)*sinh(z.im));
	std::complex<double> tz(z.re, z.im);
	std::complex<double> c = cos(tz);
	return Complex(real(c), imag(c));
}

/* squareroot 
   s = 0: sqrt(|c|)*exp(i arg(c)/2),      positive squareroot for real c 
   s = 1: sqrt(|c|)*exp(i arg(c)/2 + pi), negative squareroot for real c */
Complex sqrt(const Complex& c, int s)
{
	Complex z;
	double r, phi = 0.0;
	r = sqrt(c.abs());
	if( s <= 0) phi = c.arg()/2.0; 
	else        phi = c.arg()/2.0+M_PI;
	z = expi(phi)*r;
	return z;
}
Complex sqrt_pos(const Complex& c) 
{ 
	return sqrt(c, 0);
}
Complex sqrt_neg(const Complex& c) 
{ 
	return sqrt(c, 1);
}
/* principal value; +/-sqrt(|c|)*exp(i arg(c)/2),
   sign such that real part of the result is positive;
   sign of pure imaginary result z.im equal to the sign of c.im */
Complex sqrt(const Complex& c)
{
/*
	Complex z;
	z = expi(c.arg()/2.0 )*sqrt(c.abs());
	if(z.re == 0.0)
	{
		if(c.im > 0.0)
		{
			if(z.im >= 0.0) return z;
			else            return -z;
		}
		if(z.im >= 0.0) return -z;
		else            return z;
	}
	if(z.re >= 0.0) return z;
	return -z;
*/
	std::complex<double> tc(c.re, c.im);
	std::complex<double> z = sqrt(tc);
	return Complex(real(z), imag(z));
}

/* natural logarithm, imaginary part in ]-PI, PI] */
Complex log(const Complex& z)
{
//	return Complex(log(z.abs()), z.args());
	std::complex<double> tz(z.re, z.im);
	std::complex<double> l = log(tz);
	return Complex(real(l), imag(l));
}

/* ------------------------------------------------------------------------- */
/* powers of complex numbers */

Complex pow(const Complex& z, const int& n)
{
	if(z.re == 0.0 && z.im == 0.0) return CC0;
	Complex r;
	int i = n;
	if(n == 0) return CC1;
	if(n>0)
	{
		r = z; --i;
		while(i>0) { r = r*z; --i; }
		return r;
	}
	r = z; ++i;
	while(i<0) { r = r*z; ++i; }
	return CC1/r;
}
Complex pow(const double& d, const Complex& c)
{
//	if(d == 0.0) return CC0;
//	if(d >  0.0) return exp(c*log(d));
//	return exp(c*(log(-d)+CCI*M_PI));
	std::complex<double> tc(c.re, c.im);
	std::complex<double> p = pow(d, tc);
	return Complex(real(p), imag(p));
}
Complex pow(const Complex&z, const double& d)
{
//	if(z.re == 0.0 && z.im == 0.0) return CC0;
//	return exp(d*Complex(log(z.abs()), z.args()));
	std::complex<double> tz(z.re, z.im);
	std::complex<double> p = pow(tz, d);
	return Complex(real(p), imag(p));
}
Complex pow(const Complex& z, const Complex& c)
{
//	if(z.re == 0.0 && z.im == 0.0) return CC0;
//	return exp(c*Complex(log(z.abs()), z.args()));
	std::complex<double> tz(z.re, z.im);
	std::complex<double> tc(c.re, c.im);
	std::complex<double> p = pow(tz, tc);
	return Complex(real(p), imag(p));
}

/* ckeck validity: re, im in {nan, -nan, inf, -inf} is not OK */
bool Complex::isOK() const
{
	if(std::isfinite(re) && std::isfinite(im)) return true;
	return false;
}
bool Complex::isnotOK() const
{
	if(!std::isfinite(re) || !std::isfinite(im)) return true;
	return false;
}