#ifndef COMPLEX_H
#define COMPLEX_H

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

/* 
 * complex.h
 * Handling complex numbers 
 */

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

/* Complex: complex number */
class Complex 
{
public:
	// real & imaginary part
	double re;
	double im;

	// initialize
	inline Complex() { re=0.0; im=0.0; }
	inline Complex(const double& x) { re=x; im=0.0; }
	inline Complex(const double& x, const double& y) { re=x; im=y; }

	// copy & assignment
	inline Complex(const Complex& z) { re = z.re; im = z.im; }
	inline Complex& operator=(const Complex& z)
		{ if(this != &z) { re = z.re; im = z.im; } return *this; }

	// the negative 
	inline Complex operator-() const { return Complex(-re, -im); }

	// complex conjugation 
	inline Complex conj() const { return Complex(re, -im); }

	// modulus
	inline double abs() const { return sqrt(re*re+im*im); }

	// modulus squared
	inline double sqabs() const { return re*re+im*im; }

	// argument;
	double arg() const;  // range [0, 2 PI[
	double args() const; // range ]-PI, PI]

	// multiply by (-1) & assignment  
	inline void neg() { re = -re; im = -im; } 
	// multiply by (i) & assignment  
	inline void muli() { double t=re; re = -im; im = t; } 
	// multiply by (-i) & assignment  
	inline void mulmi() { double t=re; re = im; im = -t; } 
	// 1/(z) & assignment  
	inline void inv() { double t = re*re+im*im; re = re/t; im = -im/t; } 

	// addition 
	inline Complex operator+(const Complex& z) const
		{ return Complex(re+z.re, im+z.im); }
	inline Complex operator+(const double& x) const
		{ return Complex(re+x, im); }
	friend inline Complex operator+(const double& x, const Complex& z)
		{ return Complex(x+z.re, z.im); }

	// addition & assignment
	inline void operator+=(const Complex& z) { re+=z.re; im+=z.im; return;}
	inline void operator+=(const double& x)  { re+=x; return; }

	// subtraction 
	inline Complex operator-(const Complex& z) const
		{ return Complex(re-z.re, im-z.im); }
	inline Complex operator-(const double& x) const
		{ return Complex(re-x, im); }
	friend inline Complex operator-(const double& x, const Complex& z)
		{ return Complex(x-z.re, -z.im); }

	// subtraction & assignment 
	inline void operator-=(const Complex& z)
		{ re-=z.re; im-=z.im; return; }
	inline void operator-=(const double& x)
		{ re-=x; return; }

	// multiplication 
	inline Complex operator*(const Complex& z) const
		{ return Complex(re*z.re-im*z.im, re*z.im+im*z.re); }
	inline Complex operator*(const double& x) const
		{ return Complex(re*x, im*x); }
	friend inline Complex operator*(const double& x, const Complex& z)
		{ return Complex(x*z.re, x*z.im); }

	// multiplication & assignment
	inline void operator*=(const Complex& z)
	{ 	
		double r = re, i = im; 
		re = r*z.re-i*z.im;
		im = r*z.im+i*z.re; 
		return;
	}
	inline void operator*=(const double& x) { re*=x; im*=x; return; }

	// division 
	inline Complex operator/(const Complex& z) const
	{
		double den=z.re*z.re+z.im*z.im;
		return Complex((re*z.re+im*z.im)/den, (im*z.re-re*z.im)/den);
	}
	inline Complex operator/(const double& x) const
	{
		return Complex(re/x, im/x);
	}
	friend inline Complex operator/(const double& x, const Complex& z)
	{
		double den=z.re*z.re+z.im*z.im;
		return Complex(x*z.re/den, -x*z.im/den);
	}

	// division & assignment
	inline void operator/=(const Complex& z)
	{ 	
		double r = re, i = im; 
		double den=z.re*z.re+z.im*z.im;
		re = (r*z.re+i*z.im)/den; 
		im = (i*z.re-r*z.im)/den;
		return;
	}
	inline void operator/=(const double& x)
		{ re/=x; im/=x; return; }

	// complex conjugate & multiply
	inline Complex conjmult(const Complex& z)
		{ return Complex(re*z.re+im*z.im, re*z.im-im*z.re); }
	inline Complex conjmult(const double& x)
		{ return Complex(re*x, -im*x); }

	// input 
	void read(FILE *dat);
	// output 
	void write(FILE *dat) const;
	void nice(FILE *dat) const;

	// ckeck validity: re, im in {nan, -nan, inf, -inf} is not OK
	bool isOK() const;
	bool isnotOK() const;
};

/* some important constants ... */
#define CC0 (Complex(0.0, 0.0))
#define CC1 (Complex(1.0, 0.0))
#define CCm1 (Complex(-1.0, 0.0))
#define CCi (Complex(0.0, 1.0))
#define CCI (Complex(0.0, 1.0))
#define CCmi (Complex(0.0, -1.0))
#define CCmI (Complex(0.0, -1.0))

/* exponential function */
/* exp(c.re)*(cos(c.im) + i sin(c.im)) */
inline Complex exp(const Complex& c)
{
	double e = exp(c.re);
	return Complex(e*cos(c.im), e*sin(c.im));
}
/* cos(phi) + i sin(phi) */
inline Complex expi(double phi)
{
	return Complex(cos(phi), sin(phi));
}
/* exp(-c.im)*(cos(c.re) + i sin(c.re)) */
inline Complex expi(const Complex& c)
{
	double e = exp(-c.im);
	return Complex(e*cos(c.re), e*sin(c.re));
}
/* exp(c.im)*(cos(c.re) - i sin(c.re)) */
inline Complex expmi(const Complex& c)
{
	double e = exp(c.im);
	return Complex(e*cos(c.re), -e*sin(c.re));
}

/* trigonometric functions */
Complex sin(const Complex& z);
Complex cos(const Complex& z);

/* 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 sqrt_pos(const Complex& c);
Complex sqrt_neg(const Complex& c);
/* 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);

/* natural logarithm, imaginary part in ]-PI, PI] */
Complex log(const Complex& z);

/* powers of complex numbers */
Complex pow(const Complex& z, const int& n);
Complex pow(const double& d, const Complex& c);
Complex pow(const Complex&z, const double& d);
Complex pow(const Complex& z, const Complex& c);

#endif   // COMPLEX_H