#include <iostream>
#include <cstdlib>

#include <gsl/gsl_sf.h>

#include "Ylm.h"

using namespace std ;

#if 0
Ylm::Ylm() : l(0), m(0), c(0.)
{
}
#endif

Ylm::Ylm(const int orbl, const int magm, const complex<double> coef) : l(orbl), m(magm), c(coef)
{
	if ( abs(magm) > orbl )
	{
		cerr << "Ylm ctor invalid m= " << magm << " l= "  << orbl << endl ;
		exit(EXIT_FAILURE) ;
	}
	if ( orbl < 0 )
	{
		cerr << "Ylm ctor invalid l= " << orbl << endl ;
		exit(EXIT_FAILURE) ;
	}
	// cout << __LINE__ << " ini "  << coef << endl ;
}

#if 0
/** complex conjugation, including the coefficient
* Edmonds, (2.5.6)
* Y(l,-m) = (-1)^m *cc [Y(l,m)]
* replaced by the public non-member function below.
*/
Ylm::Ylm conj() const
{
	if ( abs(m) %2 )	// (-1)^m=-1
		return Ylm(l,-m,-std::conj(c)) ;
	else		// (-1)^m=+1
		return Ylm(l,-m,std::conj(c)) ;
}
#endif

Ylm & Ylm::operator*= (const complex<double> & fact)
{
	// cout << __LINE__ << " mul "  << fact << endl ;
	c *= fact ;
	return *this ;
}

Ylm operator*(const complex<double> & fact, const Ylm & right)
{
	// return Ylm(right.l, right.m, right.c*fact) ;
	Ylm resul(right) ;
	return (resul *= fact) ;
}

/** complex conjugation, including the coefficient
* Edmonds, (2.5.6)
* Y(l,-m) = (-1)^m *cc [Y(l,m)]
*/
Ylm conj(const Ylm & arg)
{
	Ylm resul(arg.l, -arg.m, std::conj(arg.c)) ;
	if ( abs(arg.m) %2 )	// (-1)^m=-1
		resul.c *= -1 ;
	return resul ;
}

/** Wigner 3J coefficient in the Edmonds version, mapped on the corresponding
* version by Abramowitz-Stegun implemented in the GSL
*/
static double Wigner3j(const int j1, const int j2, const int j3, const int m1, const int m2, const int m3)
{
	const double resul = gsl_sf_coupling_3j(2*j1,2*j2,2*j3,2*m1,2*m2,-2*m3)/sqrt(2.*j3+1) ;
	return ((j1+j2+m3) % 2) ? -resul : resul ;
}

/**
* Product of the left Ylm, not taken complex conjugated, by the right Ylm.
* Edmonds (4.6.5):
* Y(l1,m1)*Y(l2,m2)=sum(l,m) sqrt[(2l1+1)(2l2+1)(2l+1)/4/pi] j3(l1 l2 l, m1 m2 m)*j3(l1 l2 l,0 0 0)*cc[Y(l,m)]
* j3(j1 j2 j3, m1 m2 m3) = (-1)^(j1-j2-m3)*CGordan<j1 m1 j2 m2|j3 -m3>/sqrt(2*j3+1)
*/
vector<Ylm> operator*(const Ylm & left, const Ylm & right)
{
	vector<Ylm> resul ;
	// selection rule of Wigner-3j coefficients: sum of the m must be zero to create nonzero first j3(...)
	const int m = -left.m-right.m ;
	// skip by 2 units since j1+j2+j3 odd means that the 2nd Wigner factor is zero
	for(int l = abs(left.l-right.l) ; l <= left.l+right.l ; l += 2)
	{
		if ( abs(m) <= l )
		{
			const double pref= sqrt((left.l+0.5)*(right.l+0.5)*(2*l+1)/M_PI)
				*Wigner3j(left.l,right.l,l,left.m,right.m,m) *Wigner3j(left.l,right.l,l,0,0,0) ;
			const complex<double> c = left.c*right.c*pref  ;
		
			Ylm term(l,m) ;
			// cout << __LINE__ << " pro " << c << " fro " << left.c << " " << right.c << endl ;
			resul.push_back( c*conj(term) ) ;
		}
	}
	return resul ;
}

/**
* equality does not check for the same coefficient in front, only for the same type.
*/
bool operator== (const Ylm & first, const Ylm & secnd)
{
	return ( first.l == secnd.l && first.m == secnd.m ) ;
}

bool operator!= (const Ylm & first, const Ylm & secnd)
{
	return ( first.l != secnd.l || first.m != secnd.m ) ;
}

ostream & operator<<(ostream &os, const Ylm & someylm)
{
	os << " (" << someylm.l << "," << someylm.m << ") " ;
	return os ;
}