#ifndef COMPLEX_H
#define COMPLEX_H

/*
 * METRIC --- Mode expansion modeling in integrated optics / photonics
 * Manfred Hammer
 * University of Twente, Department of Applied Mathematics 
 * P.O. Box 217, 7500AE Enschede, The Netherlands
 * (2004)
 */

/* 
 * complex.c
 * Handling complex numbers 
 */

#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; }

	// inversion 
	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;

	// 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; } 

	// 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;
};

/* 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));
}

/* 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);

#endif   // COMPLEX_H