#include <iostream>
#include <cstdlib>

#include <complex>
#include <vector>

#include "pSTO.h"

/* use the Nnm factor and work with normalized basis
*/
#define PSTO_NORMALIZED

using namespace std ;

#if 0
// copy constructor not needed because memberwise copy suffices
// Empty constructor without proper initialization. Do some minimum amount of work perhaps to recover from
// undetected hidden use by initialization by an invalid 'n'.
pSTO::pSTO() : Ylm()
{
	n= 999 ;
}
#endif

pSTO::pSTO(const int nn, const int ll, const int mm, const double expo, const complex<double> coef) 
		: Ylm(ll,mm,coef), n(nn), alpha(expo)
{
	if( nn >=0 && expo > 0. && abs(mm)<= ll)
		Nnm = normf(expo,n) ;
	else
	{
		cerr << "Faulty initialization pSTO n= " << nn << " m= " << mm << endl ;
		exit(EXIT_FAILURE) ;
	}
}

double pSTO::normf(const double alpha, const int n)
{
#ifdef PSTO_NORMALIZED
	if( n <= NMAX)
		return Nntab[n]*pow(2.*alpha,n+0.5) ;
	else
	{
		cerr << "Unknown Nnm for n= " << n << endl ;
		exit(EXIT_FAILURE) ;
	}
#else
	return 1.0 ;
#endif
}

pSTO & pSTO::operator*= (const double & fact)
{
	c *= fact ;
	return *this ;
}

pSTO & pSTO::operator*= (const complex<double> & fact)
{
	c *= fact ;
	return *this ;
}

void pSTO::shiftOrb(const int p)
{
	m += p ;
	while ( m > l)
		m -= 2*l+1 ;
}

// Nuclear attraction integral of a single primitive. No cc anywhere nor a factor Z of the nucleus.
// Integral d^3r (1/r)*in.c*Ylm*r^(n-1)*exp(-alpha*r) 
// = delta_{l,0}*delta_{m,0}*Int r^2 r^(n-2)exp(-alpha*r)dr int_0^pi sin(theta)dtheta int_0^(2pi} dphi 1/sqrt(4pi)
// [where the 1/sqrt(4pi) is Y00 and the integral over dOmega yields a factor 4*pi] 
// = delta(l,0)delta(m,0) Int r^n exp(-alpha*r)*2*sqrt(pi)
complex<double> pSTO::NpSTO() const
{
	if( l != 0 || m != 0 )	// normalization as if Y00=1/sqrt(4pi)
		return 0. ;
	else
	{
		return 4./M_2_SQRTPI*c* Nnm *radiInt(n+1,alpha) ;
	}
}

/** A table or sequence of 1/sqrt[(2n)!], n=0,1,2,...
* Maple code to establish this
> interface(prettyprint=0) ;
> Digits := 30 ;
> for n from 0 to 11 do
>    evalf(1./( (2*n)!)^(1/2)) ;
> end do ;
*/
const double pSTO::Nntab[NMAX+1] = { 
	1.,
	.707106781186547524400844362105,
	.204124145231931508183107006226,
	.372677996249964949401528944789e-1,
	.498011920555997349987007158205e-2,
	.524950656957260037797939633659e-3,
	.456910899277617353044394667089e-4,
	.338684891864543091458842524507e-5,
	.218620157633905874407878632395e-6,
	.124976825726157033253614529760e-7,
	.641117588536780424092550271881e-9
} ;
#undef NMAX

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

// allow comparison operator to handle cases with slight numerical jitters in the exponents; otherwise the
// strict member-by-member comparison would work.  No need to compare Nnm or c.
bool operator== (const pSTO & first, const pSTO & secnd)
{
	if( first.n != secnd.n || Ylm(first) != Ylm(secnd) )
			return false ;
	// assuming DBL_DIG=15 digits (as in limits.h) of precision of a "double"
	if( fabs(first.alpha-secnd.alpha) > 1.e-14 * first.alpha)
		return false ;
	else
		return true ;
}
ostream & operator<<(ostream &os, const pSTO & someSTO)
{
	os << " [" << someSTO.n << " " << someSTO.l << " " << someSTO.m << ", " << someSTO.alpha << "] " ;
	return os ;
}

pSTO operator*(const complex<double> & fact, const pSTO & in)
{
	return pSTO(in.n, in.l, in.m, in.alpha, in.c*fact) ;
}

/** product of two STO's. The left one is *not* taken conj().
*/
vector<pSTO> operator*(const pSTO & left, const pSTO & right)
{
	const double gamma= left.alpha + right.alpha ;	// combined exponential
	const int n= left.n + right.n-1 ;		// combined new main r-power
	// cout << __LINE__ <<  " " << left << " " << right << endl;
	const vector<Ylm> ylm = Ylm(left)*Ylm(right) ;	// product expansion of the two spherical harmonics
	vector<pSTO> resul ;
		/* the original factors are left.Nnm and right.Nnm, the implicit factor in the product is
		*  normf(gamm,n), and the dangling portion is added to ylm[].c
		*/
	const double nresidNm = left.Nnm * right.Nnm/ pSTO::normf(gamma,n) ;
	for(int i=0 ; i < (int) ylm.size() ; i++)
		resul.push_back( pSTO(n, ylm[i].l, ylm[i].m, gamma, ylm[i].c * nresidNm ) ) ;
	return resul ;
}

/** A table of factorials, 0!, 1!, 2! etc.
*/
static const double facul[]={1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880.,3628800.,39916800.,479001600.,6227020800.,87178291200.};

#define FACUL (sizeof(facul)/sizeof(double))

// int(x=0..inf) x^(n-1) exp(-alpha*x) = Gamma(n)/alpha^n=(n-1)!/alpha^n, Gradsteyn-Ryzhik 3.381.4
double radiInt(const int n, const double alpha)
{
	if( n-1 >= FACUL || n-1 < 0)
	{
		cerr << "radiInt() facul[] n " << n << endl ;
		exit(EXIT_FAILURE) ;
	}
	return facul[n-1]/pow(alpha,n) ;
} 
#undef FACUL

// overlap integral of two primitives. The left one is taken cc implicitly.
complex<double> OpSTO(const pSTO & left, const pSTO &right)
{
	/* orthogonality of the Ylm
	*/
	if( left.l != right.l || left.m != right.m )
		return 0. ;
	else
	{
		// component r^(left.n-1)r^(right.n-1) and Jacobian r^2: need r^(left.n+right.n)=r^(n+1)
		return radiInt(right.n+left.n+1,right.alpha+left.alpha) *conj(left.c)*left.Nnm *right.c *right.Nnm;
	}
}

// Kinetic energy integral of two primitives. The left is implicitly converted to its cc.
// The factor (-1/2) is already included.
// (-1/2)Laplace right = [p_r^2/2 + L^2/(2r^2)] right  (Messiah 9.11)
// p_r^2/2 = -(1/2r) d^2/dr^2(r * right) after (9.18) where right = c*r^(n-1)exp(-alpha*r)*Ylm(phi,theta)
// r*right = c*r^n exp(-alpha*r)*Ylm(phi,theta)
// d^2(r*right)/dr^2 = c*Ylm(phi,theta)*[n(n-1)r^(n-2) -2*alpha*n*r^(n-1)+alpha^2*r^n]exp(-alpha*r)
// (-1/2r)*d^2(r*right)/dr^2 = (-c/2)*Ylm(phi,theta)*[n(n-1)r^(n-3) -2*alpha*n*r^(n-2)+alpha^2*r^(n-1)]exp(-alpha*r)
// L^2/(2r^2) right = (c/2)*l(l+1)*Ylm(phi,theta)*r^(n-3)exp(-alpha*r)
// The phi-theta integral over the product of the Ylm is the delta-condidition on the l and m pairs.
// Remains with Jacobian r^2: int r^2 cc[left.c] r^(left.n-1)exp(-left.apha*r)exp(-right.alpha*r)
//            *[-right.c/2{n(n-1)r^(n-3)-2*alpha*n*r^(n-2)+alpha^2*r^(n-1)}+right.c/2*right.l(right.l+1)*r^(n-3)]
// =
// int cc[left.c]*right.c r^(left.n-1)exp(-{left.apha+right.apha}*r)
//            *[-1/2{n(n-1)r^(right.n-1)-2*alpha*n*r^right.n+alpha^2*r^(right.n+1)}+1/2*right.l(right.l+1)*r^(n-1)]
// =
// int cc[left.c]*right.c/2 r^(left.n-1)exp(-{left.apha+right.apha}*r)
//            *[-{n(n-1)r^(right.n-1)-2*alpha*n*r^right.n+alpha^2*r^(right.n+1)}+right.l(right.l+1)*r^(n-1)]
complex<double> TpSTO(const pSTO &left, const pSTO &right)
{
	// using cc(left) we have the standard orthogonality relation between the Ylm components
	// which survive the Laplace operator as seen above.
	if( left.l != right.l || left.m != right.m)
		return 0. ;
	else
	{
		const double gam = right.alpha + left.alpha ;
		const int n = right.n + left.n-1 ;

		// component -n(n-1)r^(left.n-1)r^(right.n-1)
		// and component l*(l+1)*r^(left.n-1)r^(right.n-1)
		double resul = (right.l*(right.l+1)-right.n*(right.n-1))*radiInt(n,gam) ;

		// component +2*alpha*right.n*r^(left.n-1)r^(right.n)
		resul += 2.*right.alpha*right.n*radiInt(n+1,gam) ;

		// component -alpha^2*r^(left.n-1)r^(right.n+1)
		resul -= pow(right.alpha,2)*radiInt(n+2,gam) ;

		return 0.5*resul *conj(left.c)*left.Nnm * right.c*right.Nnm ;
	}
}

// Nuclear attraction integral. No cc or nuclear charge factor included.
complex <double> NlinSTO(const vector<pSTO> & in)
{
	complex<double> resul=0. ;
	for(int indx=0 ; indx< (int)in.size(); indx++)
		resul += in[indx].NpSTO() ;
	return resul ;
}

// Nuclear attraction integral. The left is taken cc implicitly.
// Integral d^3r cc[left] (1/r) right .
complex<double> NpSTO(const pSTO & left, const pSTO & right)
{
	if( left.l != right.l || left.m != right.m ) // orthogonality
		return 0. ;
	else
	{
		/* factors Nnm are evaluated in this overloaded multiplication operator 
		*/
		vector<pSTO> prod = conj(left)*right ;
		return NlinSTO(prod) ;
	}
}

#define MATH_HERON1(c0,c1,x) ((c1)*(x)+(c0))
#define MATH_HERON2(c0,c1,c2,x) (((c2)*(x)+(c1))*(x)+(c0))
#define MATH_HERON3(c0,c1,c2,c3,x) ((((c3)*(x)+(c2))*(x)+(c1))*(x)+(c0))
#define MATH_HERON4(c0,c1,c2,c3,c4,x) (((((c4)*(x)+(c3))*(x)+(c2))*(x)+(c1))*(x)+(c0))
#define MATH_HERON5(c0,c1,c2,c3,c4,c5,x) (MATH_HERON4(c1,c2,c3,c4,c5,x)*(x)+(c0))
#define MATH_HERON6(c0,c1,c2,c3,c4,c5,c6,x) (MATH_HERON5(c1,c2,c3,c4,c5,c6,x)*(x)+(c0))
#define MATH_HERON7(c0,c1,c2,c3,c4,c5,c6,c7,x) (MATH_HERON6(c1,c2,c3,c4,c5,c6,c7,x)*(x)+(c0))
#define MATH_HERON8(c0,c1,c2,c3,c4,c5,c6,c7,c8,x) (MATH_HERON7(c1,c2,c3,c4,c5,c6,c7,c8,x)*(x)+(c0))
#define MATH_HERON9(c0,c1,c2,c3,c4,c5,c6,c7,c8,c9,x) (MATH_HERON8(c1,c2,c3,c4,c5,c6,c7,c8,c9,x)*(x)+(c0))
#define MATH_HERON10(c0,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,x) (MATH_HERON9(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,x)*(x)+(c0))
#define MATH_HERON11(c0,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,x) (MATH_HERON10(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,x)*(x)+(c0))
#define MATH_HERON12(c0,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,x) (MATH_HERON11(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c10,x)*(x)+(c0))

/* computed by residDel.mp:
* 8*(n-l)!*(n'-l)!*(2*alpha)^n/alpha^l (2*alpha')^n'/alpah'^l*int_0^infty sum(s,s') k^(2l)*omega(s,n,l)*omega(s',n',l)/()^n+1-s)...
* That is the Coulomb integral without the factors of the coefficients.
*/
static double CpSTOaux( const double al, const double alpr, const int n, const int npr, const int l)
{
	/* tabulated below only for n'<=n. Exchange by symmetry for the others.
	*/
	if ( npr > n )
		return CpSTOaux(alpr,al,npr,n,l) ;

	const double x = alpr/al ;
	double res = M_PI/(al*alpr*pow(al+alpr,n+npr+1)) ;
	switch(n)
	{
#if 0
	// correct but never needed
	case 0 :
		switch(npr)
		{
		case 0 :
			switch(l)
			{
			case 0 :
				res *= 4. ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		default:
			cerr << __FILE__ << " " << __LINE__ << " n= " << n << " not implemented\n" ;
			exit(EXIT_FAILURE) ;
		}
		break ;
#endif
	case 1 :
		switch(npr)
		{
#if 0
		// correct but never needed
		case 0 :
			switch(l)
			{
			case 0 :
				res *= 4.*(alpr+2.*al)/al ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
#endif
		case 1 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON2(8,24,8,x)/x ;
				break ;
			case 1 :
				res *= 8. ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		default:
			cerr << __FILE__ << " " << __LINE__ << " n= " << n << " not implemented\n" ;
			exit(EXIT_FAILURE) ;
		}
		break ;
	case 2 :
		switch(npr)
		{
#if 0
		// correct but never needed
		case 0 :
			switch(l)
			{
			case 0 :
				res *= 8.*(alpr*alpr+3.*al*alphS)/(al*al) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
#endif
		case 1 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON3(24,96,64,16,x)/x ;
				break ;
			case 1 :
				res *= MATH_HERON1(32,8,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 2 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON4(48,240,480,240,48,x)/(x*x) ;
				break ;
			case 1 :
				res *= MATH_HERON2(32,160,32,x)/x ;
				break ;
			case 2 :
				res *= 96.;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		default:
			cerr << __FILE__ << " " << __LINE__ << " n= " << n << " not implemented\n" ;
			exit(EXIT_FAILURE) ;
		}
		break ;
	case 3 :
		switch(npr)
		{
#if 0
		// correct but never needed
		case 0 :
			switch(l)
			{
			case 0 :
				res *= 24.*(alpr+2.*al)*(alpr*alpr+2.*al*alphS)/pow(al,3) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
#endif
		case 1 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON4(96,480,480,240,48,x)/x ;
				break ;
			case 1 :
				res *= MATH_HERON2(160,80,16,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 2 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON5(192,1152,2880,2160,864,144,x)/(x*x) ;
				break ;
			case 1 :
				res *= MATH_HERON3(160,960,384,64,x)/x ;
				break ;
			case 2 :
				res *= MATH_HERON1(576,96,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 3 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON6(576,4032,12096,20160,12096,4032,576,x)/pow(x,3) ;
				break ;
			case 1 :
				res *= MATH_HERON4(320,2240,6720,2240,320,x)/(x*x) ;
				break ;
			case 2 :
				res *= MATH_HERON2(576,4032,576,x)/x ;
				break ;
			case 3 :
				res *= 2880. ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		default:
			cerr << __FILE__ << " " << __LINE__ << " n= " << n << " not implemented\n" ;
			exit(EXIT_FAILURE) ;
		}
		break ;
	case 4 :
		switch(npr)
		{
		case 1 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON5(480,2880,3840,2880,1152,192,x)/x ;
				break ;
			case 1 :
				res *= MATH_HERON3(960,720,288,48,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 2 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON6(960,6720,20160,20160,12096,4032,576,x)/(x*x) ;
				break ;
			case 1 :
				res *= MATH_HERON4(960,6720,4032,1344,192,x)/x ;
				break ;
			case 2 :
				res *= MATH_HERON2(4032,1344,192,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 3 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON7(2880,23040,80640,161280,129024,64512,18432,2304,x)/(x*x*x) ;
				break ;
			case 1 :
				res *= MATH_HERON5(1920,15360,53760,26880,7680,960,x)/(x*x) ;
				break ;
			case 2 :
				res *= MATH_HERON3(4032,32256,9216,1152,x)/x ;
				break ;
			case 3 :
				res *= MATH_HERON1(23040,2880,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 4 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON8(11520,103680,414720,967680,1451520,967680,414720,103680,11520,x)/pow(x,4) ;
				break ;
			case 1 :
				res *= MATH_HERON6(5760,51840,207360,483840,207360,51840,5760,x)/pow(x,3) ;
				break ;
			case 2 :
				res *= MATH_HERON4(8064,72576,290304,72576,8064,x)/(x*x) ;
				break ;
			case 3 :
				res *= MATH_HERON2(23040,207360,23040,x)/x ;
				break ;
			case 4 :
				res *= 161280. ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		default:
			cerr << __FILE__ << " " << __LINE__ << " n= " << n << " not implemented\n" ;
			exit(EXIT_FAILURE) ;
		}
		break ;
	case 5 :
		switch(npr)
		{
		case 1 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON6(2880,20160,33600,33600,20160,6720,960,x)/x ;
				break ;
			case 1 :
				res *= MATH_HERON4(6720,6720,4032,1344,192,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 2 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON7(5760,46080,161280,201600,161280,80640,23040,2880,x)/(x*x) ;
				break ;
			case 1 :
				res *= MATH_HERON5(6720,53760,43008,21504,6144,768,x)/x ;
				break ;
			case 2 :
				res *= MATH_HERON3(32256,16128,4608,576,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 3 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON8(17280,155520,622080,1451520,1451520,967680,414720,103680,11520,x)/pow(x,3) ;
				break ;
			case 1 :
				res *= MATH_HERON6(13440,120960,483840,322560,138240,34560,3840,x)/(x*x) ;
				break ;
			case 2 :
				res *= MATH_HERON4(32256,290304,124416,31104,3456,x)/x ;
				break ;
			case 3 :
				res *= MATH_HERON2(207360,51840,5760,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 4 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON9(69120,691200,3110400,8294400,14515200,12096000,6912000,2592000,576000,57600,x)
						/pow(x,4) ;
				break ;
			case 1 :
				res *= MATH_HERON7(40320,403200,1814400,4838400,2764800,1036800,230400,23040,x)/pow(x,3) ;
				break ;
			case 2 :
				res *= MATH_HERON5(64512,645120,2903040,1088640,241920,24192,x)/(x*x) ;
				break ;
			case 3 :
				res *= MATH_HERON3(207360,2073600,460800,46080,x)/x ;
				break ;
			case 4 :
				res *= MATH_HERON1(1612800,161280,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 5 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON10(345600,3801600,19008000,57024000,114048000,159667200,114048000,57024000,19008000,3801600,345600,x)
						/pow(x,5) ;
				break ;
			case 1 :
				res *= MATH_HERON8(161280,1774080,8870400,26611200,53222400,26611200,8870400,1774080,161280,x)
						/pow(x,4) ;
				break ;
			case 2 :
				res *= MATH_HERON6(193536,2128896,10644480,31933440,10644480,2128896,193536,x)/pow(x,3) ;
				break ;
			case 3 :
				res *= MATH_HERON4(414720,4561920,22809600,4561920,414720,x)/(x*x) ;
				break ;
			case 4 :
				res *= MATH_HERON2(1612800,17740800,1612800,x)/x ;
				break ;
			case 5 :
				res *= 14515200 ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		default:
			cerr << __FILE__ << " " << __LINE__ << " n= " << n << " not implemented\n" ;
			exit(EXIT_FAILURE) ;
		}
		break ;
	case 6 :
		switch(npr)
		{
		case 1 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON7(20160,161280,322560,403200,322560,161280,46080,5760,x)/x ;
				break ;
			case 1 :
				res *= MATH_HERON5(53760,67200,53760,26880,7680,960,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 2 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON8(40320,362880,1451520,2177280,2177280,1451520,622080,155520,17280,x)/(x*x) ;
				break ;
			case 1 :
				res *= MATH_HERON6(53760,483840,483840,322560,138240,34560,3840,x)/x ;
				break ;
			case 2 :
				res *= MATH_HERON4(290304,193536,82944,20736,2304,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 3 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON9(120960,1209600,5443200,14515200,17418240,14515200,8294400,3110400,691200,
					69120,x)/(x*x*x) ;
				break ;
			case 1 :
				res *= MATH_HERON7(107520,1075200,4838400,4032000,2304000,864000,192000,19200,x)/(x*x) ;
				break ;
			case 2 :
				res *= MATH_HERON5(290304,2903040,1658880,622080,138240,13824,x)/x ;
				break ;
			case 3 :
				res *= MATH_HERON3(2073600,777600,172800,17280,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 4 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON10(483840,5322240,26611200,79833600,159667200,159667200,114048000,
					57024000,19008000,3801600,345600,x)/pow(x,4) ;
				break ;
			case 1 :
				res *= MATH_HERON8(322560,3548160,17740800,53222400,38016000,19008000,6336000,
					1267200,115200,x)/(x*x*x) ;
				break ;
			case 2 :
				res *= MATH_HERON6(580608,6386688,31933440,15966720,5322240,1064448,96768,x)/(x*x) ;
				break ;
			case 3 :
				res *= MATH_HERON4(2073600,22809600,7603200,1520640,138240,x)/x ;
				break ;
			case 4 :
				res *= MATH_HERON2(17740800,3548160,322560,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 5 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON11(2419200,29030400,159667200,532224000,1197504000,1916006400,
					1642291200,1026432000,456192000,136857600,24883200,2073600,x)/pow(x,5) ;
				break ;
			case 1 :
				res *= MATH_HERON9(1290240,15482880,85155840,283852800,638668800,399168000,
					177408000,53222400,9676800,806400,x)/pow(x,4) ;
				break ;
			case 2 :
				res *= MATH_HERON7(1741824,20901888,114960384,383201280,170311680,51093504,
					9289728,774144,x)/(x*x*x) ;
				break ;
			case 3 :
				res *= MATH_HERON5(4147200,49766400,273715200,82114560,14929920,1244160,x)/(x*x) ;
				break ;
			case 4 :
				res *= MATH_HERON3(17740800,212889600,38707200,3225600,x)/x ;
				break ;
			case 5 :
				res *= MATH_HERON1(174182400,14515200,x) ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		case 6 :
			switch(l)
			{
			case 0 :
				res *= MATH_HERON12(14515200,188697600.,1132185600.,4151347200.,10378368000.,18681062400.,
					24908083200.,18681062400.,10378368000.,4151347200.,1132185600.,188697600.,14515200.,x)/pow(x,6) ;
				break ;
			case 1 :
				res *= MATH_HERON10(6451200.,83865600.,503193600.,1845043200.,4612608000.,8302694400.,
					4612608000.,1845043200.,503193600.,83865600.,6451200.,x)/pow(x,5) ;
				break ;
			case 2 :
				res *= MATH_HERON8(6967296.,90574848.,543449088.,1992646656.,4981616640.,1992646656.,
					543449088.,90574848.,6967296.,x)/pow(x,4) ;
				break ;
			case 3 :
				res *= MATH_HERON6(12441600.,161740800.,970444800.,3558297600.,970444800.,161740800.,
					12441600.,x)/(x*x*x) ;
				break ;
			case 4 :
				res *= MATH_HERON4(35481600.,461260800.,2767564800.,461260800.,35481600.,x)/(x*x) ;
				break ;
			case 5 :
				res *= MATH_HERON2(174182400.,2264371200.,174182400.,x)/x ;
				break ;
			case 6 :
				res *= 1916006400. ;
				break ;
			default:
				cerr << __FILE__ << " " << __LINE__ << " n= " << n << " l=" << l << " not implemented\n" ;
				exit(EXIT_FAILURE) ;
			}
			break ;
		default:
			cerr << __FILE__ << " " << __LINE__ << " n= " << n << " not implemented\n" ;
			exit(EXIT_FAILURE) ;
		}
	default:
		cerr << __FILE__ << " " << __LINE__ << " n= " << n << " not implemented\n" ;
		exit(EXIT_FAILURE) ;
	}
	return res ;
}
#undef MATH_HERON10
#undef MATH_HERON9
#undef MATH_HERON8
#undef MATH_HERON7
#undef MATH_HERON6
#undef MATH_HERON5
#undef MATH_HERON4
#undef MATH_HERON3
#undef MATH_HERON2
#undef MATH_HERON1

// 1-particle  Coulomb integral. No factor 1/2 included, but only the bare operator d^3r d^3r'/|r-r'| .
// Introduces an implicit cc of the left factor.
// use left(k)=int exp(i k r) left(r) d^3 -> left(r)=int d^3k/(2pi)^3 exp(-i k r) left(k)and
//   the cc[left(r)]= int d^3k/(2pi)^3 exp(i k r) cc[left(k)]
// Also use the Fourer Rep of the Coulomb potential, 1/|r-r'|=4*pi*int exp[-ik(r-r')]/k^2 d^3k
// int cc[left(r)] right(r')/|r-r'|d^3r d^3r' = 4pi int d^3k/(2pi)^3 cc[left(k)] right(k)/k^2
// One may switch from k to -k in the integrand to get cc[left(-k)] right(-k) where the
// parity of Ylm means cc[left(-k)]=cc[(-)^l left(k)]=(-)^l cc[left(k)] etc...
// This is essentially the same as (5.14) in my IJQC article, and (10) in Safouhi J Comp Phys 2000.
// See also Filter PRA 18.

complex<double> CpSTO(const pSTO &left, const pSTO &right)
{
	if( left.l != right.l || left.m != right.m)
		return 0. ;
	else
	{
		const double res= CpSTOaux( left.alpha, right.alpha, left.n, right.n, left.l) ;

		if ( res > 0.)
			return conj(left.c)*left.Nnm *right.c*right.Nnm *res ;
		else
		{
			// cout << __LINE__ << " n " << left.n << " n' " << right.n << " l " << left.l << " not impl" << endl ;
			exit(EXIT_FAILURE) ;
		}
	}
}