#include <iostream>
#include <fstream>
 
#include <cstring>
#include <cstdlib>

#include "linSTO.h"

using namespace std ;

#if 0
// constructor with number of elements given as argument
linSTO::linSTO(int nolin=0) : vector<pSTO>(nolin)
{
}
#else
linSTO::linSTO()
{
}
#endif

// constructor with initialization by single pSTO instance
linSTO::linSTO(const pSTO & origSTO)
{
	push_back(origSTO) ;
}

// constructor with initialization by pSTO instances
linSTO::linSTO(const vector<pSTO> & origSTO) : ::vector<pSTO>(origSTO)
{
	// *this = origSTO ;
}

// constructor by reading the file with orbitals
linSTO::linSTO (const char  *infile)
{
	double *expos= 0 ;
	int *nvect=0,
	    *lvect=0,
	    *mvect=0;
	int no=0 ;
	ifstream f(infile,ios::in) ;
	if( !f)
	{
		cerr << "cannot open " << infile << endl ;
		exit(EXIT_FAILURE) ;
	}
	while( !f.eof() )
	{
		string oneline ;
		getline(f,oneline) ;
		if ( f.eof() )
			break; 
		if( oneline[0] == '#')		// skip comment lines
			continue ;

		// cout << __FILE__ << " " << __LINE__ << endl ;
				// do not read beyond EOD (end of definition) at the start of line
		if ( oneline.find("EOD") == 0 )
			break;
		expos= (double*)realloc(expos,(no+1)*sizeof(double)) ;
		nvect= (int*)realloc(nvect,(no+1)*sizeof(int)) ;
		lvect= (int*)realloc(lvect,(no+1)*sizeof(int)) ;
		mvect= (int*)realloc(mvect,(no+1)*sizeof(int)) ;
		istringstream line(oneline) ;
		line >> nvect[no] >> lvect[no] >> mvect[no] >> expos[no] ;
				// case of m=2l+1 denotes implicit contraction. use all |m|<=l
		if( mvect[no] == 2*lvect[no]+1)
		{
			int m=mvect[no] ;
			const int l=lvect[no] ;
			expos= (double*)realloc(expos,(no+m)*sizeof(double)) ;
			nvect= (int*)realloc(nvect,(no+m)*sizeof(int)) ;
			lvect= (int*)realloc(lvect,(no+m)*sizeof(int)) ;
			mvect= (int*)realloc(mvect,(no+m)*sizeof(int)) ;
			mvect[no++]= -l ;
			/* have already installed m=-l above...*/
			for( m = -l+1; m <= l ; m++)
			{
				nvect[no]= nvect[no-1] ;
				lvect[no]= lvect[no-1] ;
				mvect[no]= m ;
				expos[no]=expos[no-1] ;
				no++ ;
			}
		}
		else
			no++ ;
	}
	f.close() ;
	for(int indx=0 ; indx < no ; indx++)
	{
		/* This yields a basis in the sense that all the c-factors are set to one.
		* The Nnm factors act according to the definition  of PSTO_NORMALIZED defined
		* in the pSTO class.
		*/
		pSTO tmp(nvect[indx],lvect[indx],mvect[indx],expos[indx]) ;
		push_back(tmp) ;
	}
	free(mvect) ;
	free(lvect) ;
	free(nvect) ;
	free(expos) ;
}

#if 0
// destructor, superfluous
linSTO::~linSTO()
{
}
#endif

// overloaded augmentation/superposition
// If toadd.psto equals the present comp[any].psto, add
// toadd.psto.coeff to comp[any].coeff. Otherwise append toadd to the list. This means expand the existing
// linSTO by 'toadd' if there is no compatible element (w.r.t. n,,m,s,alpha ) in the existing one,
// and add compatible elements algebraicly.
linSTO & linSTO::operator+= ( const vector<pSTO> & toadd)
{
	for( int indxnew=0 ; indxnew < (int) toadd.size() ; indxnew++)
	{
		int indxpres ;
		for( indxpres=0 ;indxpres< (int)size() ; indxpres++)
			if( at(indxpres) == toadd.at(indxnew) )	// this compares l,m,n and alpha, not c or Nlm.
			{
				// The present component has the form coe[indxpres]*c where c is the member of Ylm.
				// The new one has the form toadd.coe[indxnew]*c where c is the member of its Ylm.
				at(indxpres).c += toadd.at(indxnew).c ;
				break ;
			}
		if( indxpres == (int)size() )
		{
			// didn't find toadd in the list; append by creation of a new list
			push_back(toadd.at(indxnew)) ;			// add new one
		}
	}
	return * this ;		// call linSTO  copy constructor on return
}

// unary minus
linSTO linSTO::operator-() const
{
	linSTO resul = *this ;
	for( int indx=0 ; indx < (int)size() ; indx++)
		resul.at(indx).c *=  -1 ;
	return resul ;
}

linSTO & linSTO::operator-= ( const linSTO & toadd)
{
	const linSTO negLeft(-toadd) ;
	*this += negLeft ;
	return * this ;		// call linSTO  copy constructor on return
}

linSTO & linSTO::operator*= (const complex<double> & fact)
{
	for( int indx=0 ; indx < (int)size() ; indx++)
		at(indx).c *= fact ;
	return *this ;
}

#if 0
// not needed
linSTO & linSTO::operator/= (const complex<double> & fact)
{
	for( int indx=0 ; indx < (int)size() ; indx++)
		at(indx).c /= fact ;
	return *this ;
}
#endif

// multiply two linear STO's: the 'left' factor is used without complex conjugation
linSTO & linSTO::operator*= (const pSTO & right)
{
	// linSTO resul(0) ;
	linSTO resul ;
	for(int i=0; i< (int)size()  ; i++)
	{
		vector<pSTO> tmp = at(i)*right ;
		resul += tmp ;
	}
	return *this = resul ;
}

// unary minus
void linSTO::shiftOrb(const int p)
{
	for( int indx=0 ; indx < (int)size() ; indx++)
		at(indx).shiftOrb(p) ;
}

#if 0
// Duplicate all elements to allow for a implicit spin-labelling at each second index
// altern=true means spins alternate with up,down,up,down, otherwise there are two blocks
linSTO::linSTO dup(bool altern)
{
	linSTO resul ;
	if (altern )
		for(int i=0; i< (int)size()  ; i++)
		{
			resul.push_back( at(i) ) ;
			resul.push_back( at(i) ) ;
		}
	else
	{
		for(int i=0; i< (int)size()  ; i++)
			resul.push_back( at(i) ) ;
		for(int i=0; i< (int)size()  ; i++)
			resul.push_back( at(i) ) ;
	}
	return resul ;
}
#endif


/** complex conjugation.
*/
linSTO conj(const linSTO & arg)
{
	vector<pSTO> resul ;
	for(int i=0; i < (int)arg.size() ; i++)
		resul.push_back( conj(arg[i]) ) ;
	//return linSTO(resul) ;
	return resul ;
}

ostream & operator<<(ostream &os, const linSTO & someSTO)
{
	for( int indx=0 ; indx< (int)someSTO.size() ; indx++)
		os << someSTO.at(indx).c << " : " << someSTO.at(indx) << " | " ;
	return os ;
}

linSTO operator+(const linSTO & left, const linSTO & right)
{
	linSTO resul(left) ;
	resul += right ;
	return resul ;
}
linSTO operator-(const linSTO & left, const linSTO & right)
{
	linSTO resul(left) ;
	resul -= right ;
	return resul ;
}

linSTO operator*(const complex<double> & fact, const linSTO & in)
{
	linSTO resul(in) ;
	resul *= fact ;
	return resul ;
}

linSTO operator*(const double & fact, const linSTO & in)
{
	linSTO resul(in) ;
	resul *= fact ;
	return resul ;
}

linSTO operator*(const pSTO & left, const linSTO & right)
{
	// linSTO resul(0) ;
	linSTO resul ;
	for( int j=0 ; j< (int) right.size(); j++)
		resul += left * right[j] ;
	return resul ;
}

#if 0
// not used
linSTO operator/(const linSTO &in, const complex<double> & fact)
{
	linSTO resul(in) ;
	resul /= fact ;
	return resul ;
}

linSTO operator/(const linSTO &in, const double & fact)
{
	linSTO resul(in) ;
	resul /= fact ;
	return resul ;
}
#endif

// multiply two linear STO's: the 'left' factor is not used with its complex conjugate
linSTO operator*(const linSTO & left, const linSTO & right)
{
	// linSTO resul(0) ;
	linSTO resul ;
	for(int i=0; i< (int)left.size() ;i++)
		resul += left[i] * right ;
	return resul ;
}

# if 0
// not needed
linSTO & linSTO::operator*= (const linSTO & right)
{
	linSTO resul = *this * right ;
	return *this = resul ;
}
#endif

// kinetic energy integral of two linear combinations. The left factor is taken complex conjugated explicitly.
// The factor (-1/2) in front of the Laplacian is already included.
complex<double> TlinSTO(const linSTO & left, const linSTO & right)
{
	complex<double> resul=0. ;
	for(int i=0 ; i< (int)left.size(); i++)
		for(int j=0; j< (int)right.size(); j++)
			resul += TpSTO(left[i],right[j]) ;
	return resul ;
}

// Nuclear attraction integral. The left is taken cc implicitly.
// Integral d^3r cc[left] (1/r) right .
complex <double> NlinSTO(const linSTO &left, const linSTO &right)
{
	// conj(left)*right returns a vector<pSTO> which matches the signature in pSTO.h 
	return NlinSTO( conj(left) * right) ;
}

// one-particle Coulomb potential. The left operator is used complex conjugated on-the-fly within the implementation.
// No factor 1/2 but bare 1/|r-r'| d^3r d^3r'
complex<double> ClinSTO(const linSTO & left, const linSTO & right)
{
	complex<double> resul=0. ;
	for(int leindx=0 ; leindx<left.size() ; leindx++)
		for(int reindx=0; reindx<right.size() ; reindx++)
			resul += CpSTO(left[leindx],right[reindx]) ;
	return resul ;
}