#ifndef _SS_COMPLEX_H_
#define _SS_COMPLEX_H_

//
//   File : ss_complex.h
//   Creation date : 28.03.2006 by Szymon Stefanek
//
//   Copyright (C) 2006 Szymon Stefanek (pragma at kvirc dot net)
//
//   This program is FREE software. You can redistribute it and/or
//   modify it under the terms of the GNU General Public License
//   as published by the Free Software Foundation; either version 2
//   of the License, or (at your opinion) any later version.
//
//   This program is distributed in the HOPE that it will be USEFUL,
//   but WITHOUT ANY WARRANTY; without even the implied warranty of
//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
//   See the GNU General Public License for more details.
//
//   You should have received a copy of the GNU General Public License
//   along with this program. If not, write to the Free Software Foundation,
//   Inc. ,59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//

#include <math.h>


template<class T> class SSComplex
{
protected:
	T m_tRe;
	T m_tIm;
public:
	SSComplex()
	{
	}
	SSComplex(const T &tRe,const T &tIm)
	{
		m_tRe = tRe;
		m_tIm = tIm;
	}
	SSComplex(const SSComplex<T> &d)
	{
		m_tRe = d.m_tRe;
		m_tIm = d.m_tIm;
	}
	~SSComplex()
	{
	}
public:
	const T & re() const
	{
		return m_tRe;
	}
	const T & im() const
	{
		return m_tIm;
	}
	void setRe(const T &r)
	{
		m_tRe = r;
	}
	void setIm(const T &i)
	{
		m_tIm = i;
	}
	SSComplex<T> conjugate()
	{
		return SSComplex<T>(m_tRe,-m_tIm);
	}
	void setConjugate()
	{
		m_tIm = -m_tIm;
	}
	void setOpposite()
	{
		m_tRe = -m_tRe;
		m_tIm = -m_tIm;
	}
	T modulus()
	{
		return (T)sqrt((double)((m_tRe * m_tRe) + (m_tIm * m_tIm)));
	}
	T phase()
	{
		return (T)atan((double)(m_tIm / m_tRe));
	}
	void operator += (const SSComplex<T> &y)
	{
		m_tRe += y.re();
		m_tIm += y.im();
	}
	void operator -= (const SSComplex<T> &y)
	{
		m_tRe -= y.re();
		m_tIm -= y.im();
	}
	void operator *= (const SSComplex<T> &y)
	{
		// (a+ib)(c+id) = ac + iad + ibc - bd
		// -> re = ac - bd
		// -> im = ad + bc
		m_tRe = (m_tRe * y.re()) - (m_tIm * y.im());
		m_tIm = (m_tRe * y.im()) + (m_tIm * y.re());
	}
	void operator /= (const SSComplex<T> &y)
	{
		// (a+ib)/(c+id) = (a+ib)(c-id)/(c+id)(c-id) = (a+ib)(c-id)/c^2+d^2
		// -> re = (ac + bd) / (c^2+d^2)
		// -> im = (bc - ad) / (c^2+d^2)
		T den = (y.re() * y.re()) + (y.im() * y.im());
		m_tRe = ((m_tRe * y.re()) + (m_tIm * y.im())) / den;
		m_tIm = ((m_tIm * y.re()) - (m_tIm * y.im())) / den;
	}

	void operator = (const SSComplex<T> &d)
	{
		m_tRe = d.m_tRe;
		m_tIm = d.m_tIm;
	}
	
	void operator = (const T &d)
	{
		m_tRe = d;
		m_tIm = 0;
	}
};

template<class T> SSComplex<T> operator + (const SSComplex<T> &x,const SSComplex<T> &y)
{
	return SSComplex<T>(x.re() + y.re(),x.im() + y.im());
}

template<class T> SSComplex<T> operator - (const SSComplex<T> &x,const SSComplex<T> &y)
{
	return SSComplex<T>(x.re() - y.re(),x.im() - y.im());
}

template<class T> SSComplex<T> operator * (const SSComplex<T> &x,const SSComplex<T> &y)
{
	// (a+ib)(c+id) = ac + iad + ibc - bd
	// -> re = ac - bd
	// -> im = ad + bc
	return SSComplex<T>(
			(x.re() * y.re()) - (x.im() * y.im()),
			(x.re() * y.im()) + (x.im() * y.re())
	);
}

template<class T> SSComplex<T> operator / (const SSComplex<T> &x,const SSComplex<T> &y)
{
	// (a+ib)/(c+id) = (a+ib)(c-id)/(c+id)(c-id) = (a+ib)(c-id)/c^2+d^2
	// -> re = (ac + bd) / (c^2+d^2)
	// -> im = (bc - ad) / (c^2+d^2)
	T den = (y.re() * y.re()) + (y.im() * y.im());
	return SSComplex<T>(
		((x.re() * y.re()) + (x.im() * y.im())) / den,
		((x.im() * y.re()) - (x.re() * y.im())) / den
	);
}

template<class T> bool operator == (const SSComplex<T> &x,const SSComplex<T> &y)
{
	return (x.m_tRe == y.m_tRe) && (x.m_tIm == y.m_tIm);
}

template<class T> bool operator == (const SSComplex<T> &x,const T &y)
{
	return (x.m_tRe == y) && (x.m_tIm == (T)0);
}


template<class T> bool operator != (const SSComplex<T> &x,const SSComplex<T> &y)
{
	return (x.m_tRe != y.m_tRe) || (x.m_tIm != y.m_tIm);
}

template<class T> bool operator != (const SSComplex<T> &x,const T &y)
{
	return (x.m_tRe != y) || (x.m_tIm != (T)0);
}


// helpers for the matrix class

template<class T> SSComplex<T> matrix_helper_complexConjugate(const SSComplex<T> &t){ return t.conjugate(); };


template<class T> bool matrix_helper_valueNearlyZero(SSComplex<T> &t)
{
	if(t.re() > 0)
	{
		if(t.re() > 0.0000001)return false;
	} else {
		if(t.re() < -0.0000001)return false;
	}
	if(t.im() > 0)
	{
		if(t.im() > 0.0000001)return false;
	} else {
		if(t.im() < -0.0000001)return false;
	}
	return true;
}

#endif //!_SS_COMPLEX_H_