/****************************************************************************
	Calculates features of an H2+  diatomic molecule. 
	Use method of E. Teller, Zeit. Phys. 61, 440 (1930)
	C.A. Bertulani - 05/19/2000    
*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include <iostream.h>
#include <iomanip.h>
#include <fstream.h>

typedef double Number;

#include "molecule.h"

void main(){
    Number R, nu, Np, Nm, dR, dg, mu;
	int val;
	molecule h2plus;
	Number gamma;
	cout << "\nEnter:\n";
	cout << "1 - to calculate electronic energy\n";
	cout << "2 - to solve transcendental equations for separation constant\n";
	cout << "3 - to get wavefunctions\n";
	cin >> val;
	dR = 0.1;
	R=dR;
    ofstream out1("h2plus1.dat");
    ofstream out2("h2plus2.dat");
    ofstream out3("h2plus3.dat");
	switch(val){
	  case 1:
		/* finds the exact energy distance R and compare with approximation*/
		  out1.setf(ios::fixed);
          out1.precision(2);
		  for(R=0.1;R<=7.;R+=dR){
			h2plus = molecule(R);
			out1 << R  << setw(10) << h2plus.Eplus();
			out1 << setw(10) << h2plus.Eapp(+1); 
			out1 << setw(10) << h2plus.Eminus() << setw(10) << h2plus.Eapp(-1) << endl;
		  }
		cout << "\nOutput of R, E+, EApprox+, E-, and EApprox- are in file h2plus1.dat\n";
		break;
	  case 2:
		/* This test other parts of the program */
		/* Solves transcendental eqs. (23) and (25) of  Teller */
		cout << "\n\nTest solution of the transcendental eqs. (23) and (25) of\n";
		cout <<	"E. Teller, Zeit. Phys. 61, 440 (1930)\n";
		h2plus = molecule(dR);
		dg = 0.1;
		/* this result is well fitted by the 4th-order polynomial*/
		out2.setf(ios::fixed);
		out2.precision(2);
		for(gamma=0.1;gamma<=10.;gamma+=dg){
			out2 << gamma << setw(10) << h2plus.Aplus(gamma) << setw(10) 
				<< h2plus.Afit(gamma,+1) << 
				setw(10) << h2plus.Aminus(gamma) << setw(10) << 
				h2plus.Afit(gamma,-1) << endl;
		}
		cout <<"\nOutput of gamma, A+, Approx+, A- and Approx- are in file h2plus2.dat\n";
		break;
	  default:
		cout << "Wavefunctions along line joining the two protons. Origin at the middle.\n"; 
		R=2.;
		h2plus = molecule(R);
		Np = h2plus.norm(+1);
		Nm = h2plus.norm(-1);
		out3.setf(ios::fixed);
		out3.precision(2);
		for(mu=5.;mu>1.;mu-=0.1){
			out3 << -R*mu/2. << setw(10) << h2plus.wave(mu,-1.,1)/Np << setw(10)  
				<< h2plus.wave_app(mu,-1.,1) << setw(10)  
				<< h2plus.wave(mu,-1.,-1)/Nm   << setw(10)  << h2plus.wave_app(mu,-1.,-1)<<endl;
		}
		for(nu=-1.;nu<1.;nu+=0.05){
			out3 << R*nu/2. << setw(10)  << h2plus.wave(1.,nu,+1)/Np << setw(10)  
				<< h2plus.wave_app(1.,nu,+1) << setw(10)  
				<< h2plus.wave(1.,nu,-1)/Nm   << setw(10)  << h2plus.wave_app(1.,nu,-1)<<endl; 
		}
		for(mu=1.;mu<=5.;mu+=0.1){
			out3 << R*mu/2. << setw(10)  << h2plus.wave(mu,1.,1)/Np << setw(10)  
				<< h2plus.wave_app(mu,1.,1) << setw(10)  
				<< h2plus.wave(mu,1.,-1)/Nm   << setw(10)  << h2plus.wave_app(mu,1.,-1)<<endl;
		}
		cout <<"\nOutput of r, Psi+, Psi_Approx+, Psi- and Psi_Approx- are in file h2plus3.dat\n";
		cout <<"\nNOTE: \n(1) R=2 Angst. \n(2) Wavefunction is simmetric around 0.\n";
	}
/* end tests */
}