/****************************************************************************
	This class contains the features of a H2+ molecule. 
	C.A. Bertulani - 05/20/2000   
/****************************************************************************/

typedef double Number;

class molecule{
public:

	molecule(Number);		/* Constructors */
	molecule();
	Number R();				/* Distance bewteen atoms */
	Number Eapp(int);		/* Approximate electronic energy */
	Number Eplus();			/* Exact electronic energy for gerade state */
	Number Eminus();		/* Exact electronic energy for ungerade state */
	Number Efit(int);	/* Fit to the exact energies */
	Number Aplus(Number);	/* Solution of transcendental equation for gerade state */ 
	Number Aminus(Number);	/* Solution of transcendental equation for ungerade state */ 
	Number Afit(Number, int); /* Fitting above solutions with a polynomial */
	void trans(int, Number[], Number[]);	/* Solves the transcendental eq. (gerade) */
	void trans2(int, Number[], Number[]);	/* Solves the transcendental eq. (ungerade) */
	Number un_wave(Number, Number, int );  /* Unnormalized wavefunction */	
	Number norm(int );
	Number wave(Number, Number, int );	/* Normalized wavefnction */
	Number wave_app(Number, Number, int );	/* Approximated wavefunction */
	Number a();			/* Optimized Bohr radius */
private:
	static const Number E2;
	static const Number HBARC;
	static const Number M_e;
	static const Number ME;
	static const Number a_B;
	static const Number E_0;
	Number the_R;
};


/* Constructors  */

molecule::molecule(Number R){
	the_R = fabs(R);
}

molecule::molecule(){
}
/*******************************************************************/
/* atomic constants  */
const Number molecule::E2 = 14.4;  /* elementary charge squared (eV.Angs) */
const Number molecule::HBARC = 1973 ;   /* Planck's const/(2.pi) times c */
const Number molecule::M_e = 511e3;  /* Electron mass */
const Number molecule::a_B = 0.529;  /* Bohr radius */
const Number molecule::E_0 = - 13.6;      /* Hydrogen g.s. energy */

/*******************************************************************/
/* Member functions */

/* Member function 'R': position of the molecule */

Number molecule::R(){return the_R;}

/*******************************************************************/
/* Member function 'energ': energy of the molecule */

Number molecule::Eapp(int pm){

	Number x, I, J, S, anew, Enew;
	anew = molecule::a();
	x = the_R/anew;
	I = E2/the_R*(1.-exp(-2.*x)*(1.+x));
	J = E2/anew*exp(-x)*(1.+x);
	S = exp(-x)*(1.+x+x*x/3.);
	Enew=E_0*a_B/anew*(2.*anew-a_B)/anew;
	if(pm == 1)return Enew-(I+J)/(1+S);
	else{
		Enew=Enew+(J-I)/(1-S);
		return Enew;
		}
}
/*******************************************************************/
/* Member function 'un_wave': unnormalized wavefunctions */

Number molecule::un_wave(Number mu, Number nu, int pm){

	Number  a2, a3, a4, a5, a6, a7, Ap, Am;
	Number fp, fm;
	Number gp = - M_e*molecule::Efit(+1)*the_R*the_R/2./HBARC/HBARC;
	Number gm = - M_e*molecule::Efit(-1)*the_R*the_R/2./HBARC/HBARC;
	Ap = molecule::Afit(gp,+1);
	Am = molecule::Afit(gm,-1);
	a2 = Ap/2.;
	a4 = (Ap*(Ap+6.)/2.-gp)/12.;
	a6 = ((Ap+20.)/12.*(Ap/2.*(Ap+6.)-gp)-gp*Ap/2.);
	a3 = (Am+2)/6.;
	a5 = ((Am+12.)*(Am+2.)/6.-gm)/20.;
	a7 = ((Am+30.)/20.*((Am+12.)*(Am+2.)/6.-gm)-gm*(Am+2.)/6.)/42.;
	/* wavefunctions */
	fp=(1.+a2*pow(nu,2)+a4*pow(nu,4)+a6*pow(nu,6))*exp(-sqrt(gp)*mu);
	fm=(nu+a3*pow(nu,3)+a5*pow(nu,5)+a7*pow(nu,7))*exp(-sqrt(gm)*mu);
	if(pm>0)return fp;
	else return fm;
}
/*******************************************************************/
/* Member function 'norm': norm of wavefunctions */

Number molecule::norm(int pm){

	Number  dnu, dmu, sump, summ, tempp, tempm, fp, fm, mu, nu;
	int i, j, int1;
	int n=201;
	int int2=4;
	dnu=2./(n-1);
	dmu=20./(n-1);
	sump=0.;
	summ=0.;
	for(i=1;i<=n-1;i++){
		mu=1.+i*dmu;
		/* Simpson's integration - contribution of end points */
		fp=molecule::un_wave(mu,-1.,+1);
		fm=molecule::un_wave(mu,-1.,-1);
		tempp=2.*fp*fp*(1.-mu*mu)*dnu/3.;
		tempm=2.*fm*fm*(1.-mu*mu)*dnu/3.;
		tempp=0.;
		tempm=0.;
		int1=4;
		for(j=2;j<=n-2;j++){
			nu=-1.+j*dnu;
			fp=molecule::un_wave(mu,nu,+1);
			fm=molecule::un_wave(mu,nu,-1);
			tempp=tempp+int1*fp*fp*(mu*mu-nu*nu)*dnu/3.;
			tempm=tempm+int1*fm*fm*(mu*mu-nu*nu)*dnu/3.;
			int1 = 6 - int1;		/* Alternate weights: 2 or 4 */
		}
		/* In this integral end points were neglected. 
		        This is OK for large nmax. */
		summ = summ + int2*tempm*dmu/3.;
		sump = sump + int2*tempp*dmu/3.;
		int2 = 6 - int2;
	}
	if(pm>0)return sqrt(3.14159/4.*pow(the_R,3)*fabs(sump));
	else return sqrt(3.14159/4.*pow(the_R,3)*fabs(summ));
}
/*******************************************************************/
/* Member function 'wave': normalized wavefunctions */

Number molecule::wave(Number mu1, Number nu1, int pm){

	return molecule::un_wave(mu1, nu1, pm);
}
/*******************************************************************/
/* Member function 'wave_app': approximated wavefunctions */

Number molecule::wave_app(Number mu, Number nu, int pm){

	Number x, S, anew, c, temp;
	anew = molecule::a();
	x = the_R/anew;
	S = exp(-x)*(1.+x+x*x/3.);
	c = 1./sqrt(2.*3.14159*pow(anew,3)*(1+pm*S));
	temp = exp(-x/2*mu)*(pm*exp(-nu*x/2.)+exp(x/2.*nu));
	return c*temp;
}
/*******************************************************************/
/* Member function 'a': solves transcendental eq. for optimal  Bohr radius */

Number molecule::a(){

	Number anew, aold;
    anew = a_B;
    aold = 0.;
	while(fabs(anew-aold)>1.e-6){
		 aold = anew;
		 anew = a_B/(1.+exp(- the_R / aold));
	}

	return anew;
}
/*******************************************************************************
Solves transcendental eq. (23) of E. Teller, Zeit. Phys. 61, 440 (1930) 
********************************************************************************/
Number molecule::Aplus(Number gamma){
    Number *g, *A, res;
	int N=501;
	Number dg=50./N;
	g = new Number[N];
	A = new Number[N];
	trans(N, g, A);
	int l = int(gamma/dg);   // linear interpolation
	Number al = gamma - l*dg;
	if(l>=0&&l<N-1) res = (1.-al)*A[l]+al*A[l+1];
	return res;
}
/**/
void molecule::trans(int N, Number gamma[], Number Ac[]){

	Number *a, *b, f, dAc, fold, Ac0;
	int i, n, j, k;
	n = 10;    /* number of highest b[n-1] coefficient in continued fraction */
	b = new Number[n];
	a = new Number[n];
	/* find Ac for a given gamma - simple search */
	dAc=0.1;
	Ac[0]=0.;
	Ac0=0.;
	gamma[0]=0.;
	for(k=1;k<=N-1;k++){
		gamma[k]=k*50./N;
		Ac[k]=Ac0-dAc;
		fold=0.;
		for(j=1;j<=50;j++){
			Ac[k]=Ac[k]+dAc;
			/* first build up a's and b's coefficients */
			a[0]=0.;
			a[1]=gamma[k]/6.;
			b[0]=-Ac[k]/2.;
			b[1]=1.+Ac[k]/6.;
			for(i=2;i<=n-1;i++){
				a[i]=-gamma[k]/(2*(i-1))/(2*(i-1)+3);
				b[i]=1+Ac[k]/(2*i)/(2*i+1);
			}
			/* construct fraction backwards */
			f=b[n-1];
			for(i=n-2;i>=0;i--)f = b[i]+a[i+1]/f;
			if(f*fold<0.){
				Ac[k]=Ac[k]-dAc;
				dAc=dAc/2.;
			}
			else fold=f;
			if(fabs(dAc)<1.e-6){ Ac0=Ac[k]; dAc=0.1;break;}
		}
	}
}
/***********************************************************************
Solves transcendental eq. (25) of E. Teller, Zeit. Phys. 61, 440 (1930) 
***********************************************************************/
Number molecule::Aminus(Number gamma){
    Number *g, *A, res;
	int N=501;
	Number dg=50./N;
	g = new Number[N];
	A = new Number[N];
	trans2(N, g, A);
	int l = int(gamma/dg);	// linear interpolation
	Number al = gamma - l*dg;
	if(l>=0&&l<N-1) res = (1.-al)*A[l]+al*A[l+1];
	return res;
}
/**/
void molecule::trans2(int N, Number gamma[], Number Ac[]){

	Number *a, *b, f, dAc, fold, Ac0;
	int i, n, j, k;
	n = 10;    /* number of highest b[n-1] coefficient in continued fraction */
	b = new Number[n];
	a = new Number[n];
	/* find Ac for a given gamma - simple search */
	dAc=0.3;
	Ac[0]=-2.;
	Ac0=-2.;
	gamma[0]=0.;
	for(k=1;k<=N-1;k++){
		gamma[k]=k*50./N;
		Ac[k]=Ac0-dAc;
		fold=0.;
		for(j=1;j<=50;j++){
			Ac[k]=Ac[k]+dAc;
			/* first build up a's and b's coefficients */
			a[0]=0.;
			a[1]=gamma[k]/12.;
			b[0]=-(Ac[k]+2.)/6.;
			b[1]=1.+Ac[k]/12.;
			for(i=2;i<=n-1;i++){
				a[i]=-gamma[k]/(2*i-1)/(2*i+2);
				b[i]=1+Ac[k]/(2*i+1)/(2*i+2);
			}
			/* construct fraction backwards */
			f=b[n-1];
			for(i=n-2;i>=0;i--)f = b[i]+a[i+1]/f;
			if(f*fold<0.){
				Ac[k]=Ac[k]-dAc;
				dAc=dAc/2.;
			}
			else fold=f;
			if(fabs(dAc)<1.e-6){ Ac0=Ac[k]; dAc=0.1;break;}
		}
	}
}
/*******************************************************************************
	Fits to the solution of transcendental eqs. (23) and (25) of E. Teller. 
********************************************************************************/
Number molecule::Afit(Number gamma, int pm){
    Number a0, a1, a2, a3, a4, b0, b1, b2, b3, b4, fit;
	a0= 0. ;			// symmetric state
	a1=   0.3398071573;
	a2=  0.01986316544;
	a3=   - 0.0003751132225;
	a4=   2.758904302E-006;
	b0= -2. ;			// antisymmetric state
	b1=    0.5975772439;
	b2= 0.007941739587;
	b3=   - 0.000181237624;
	b4=  1.246589754E-006;
	if(pm>0)return fit=a0+a1*gamma+a2*pow(gamma,2)+a3*pow(gamma,3)+a4*pow(gamma,4);
	else
	return fit=b0+b1*gamma+b2*pow(gamma,2)+b3*pow(gamma,3)+b4*pow(gamma,4);
}
/*******************************************************************************
	Finds the exact energy for a given R - symmetric (gerade) state
********************************************************************************/
Number molecule::Eplus(){
	Number rho=the_R/a_B;
	/* do a simple search. */
	Number dg=0.001;
	Number gam=0.;
	Number fold=0.;
	for(int j=1;j<=1000;j++){
		gam=gam+dg;
//		Number A=molecule::Afit(gam,+1);
		Number A=molecule::Aplus(gam);
		Number min =(gam+2.*sqrt(gam)+1.-A)/2./(1./2./sqrt(gam)+1.);
		Number f = rho - min;
		if(f*fold<0){
			gam=gam-dg;
			dg=dg/4.;
		}
		else {fold=f; dg = 2.*dg;}
		if(dg<1.e-6)break;
	}
	return 4.*gam/rho/rho*E_0;
}
/*******************************************************************************
	Finds the exact energy for a given R - antisymmetric (ungerade) state
********************************************************************************/
Number molecule::Eminus(){
	Number rho=the_R/a_B;
	/* do a simple search. */
	Number dg=0.001;
	Number gam=0.;
	Number fold=0.;
	for(int j=1;j<=1000;j++){
		gam=gam+dg;
//		Number A=molecule::Afit(gam,-1);
		Number A=molecule::Aminus(gam);
		Number min=1/4.*(7.*gam+8.*sqrt(gam)+4.+6.*pow(gam,1.5)+2.*gam*gam
			-2.*A*gam-2.*A*sqrt(gam)-A);
		min = min/(gam+1.5*sqrt(gam)+1.5+3./4./sqrt(gam));
		Number f = rho - min;
		if(f*fold<0){
			gam=gam-dg;
			dg=dg/4.;
		}
		else {fold=f; dg = 2.*dg;}
		if(gam>62.){cerr<<"\nSorry. Method doesn't work for this R =" <<
			the_R <<"\n" ; break;}
		if(dg<1.e-6)break;
	}
	return 4.*gam/rho/rho*E_0;
}
/*******************************************************************************
	Fits to the exact energy for a given R 
********************************************************************************/
Number molecule::Efit(int pm){
	Number a0, a1, a2, a3, a4;
	if(pm>0){
		a0=-53.83615459; 	/*     Eplus      */
		a1=31.97172445;
		a2=-10.80648891;
		a3= 1.629868337;
		a4= - 0.08990237871;
	}
	else{
		a0= -12.58636789; 	/*     Eminus      */
		a1= - 6.793532919;
		a2=2.680435046;
		a3=-0.3970004434;
		a4=0.02059612589;
	}
	return a0+a1*the_R+a2*pow(the_R,2)+a3*pow(the_R,3)+a4*pow(the_R,4);
}
/****************    End class molecule   ****************************************/