//---------------------------------------------------------------------------
#include <vcl\vcl.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream.h>
#include <fstream.h>
#include <iomanip.h>
#pragma hdrstop
#include "schroedinger.h"
//---------------------------------------------------------------------------
#pragma resource "*.dfm"
	TForm1 *Form1;
	double maxx,minx,maxy,miny;
	double xmod,ymod;

//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)

// declaration of constants
        {
		#define MHBARSQ 13.12508203
		#define ESQ 1.439964392
        }
    ofstream outFile    ;
//---------------------------------------------------------------------------
void __fastcall TForm1::pixw(double x,double y)
	{
		if(x>maxx || x<minx) return;
		if(y>maxy || y<miny) return;
		int gx=(int)((x-minx)*xmod+0.05);
		int gy=(int)((maxy-y)*ymod+0.05);
		Image1->Canvas->MoveTo(gx,gy);
		Image1->Canvas->LineTo(gx+1,gy+1);
	}
//---------------------------------------------------------------------------

void TForm1::initgraph(double xscal,double yscal)
	{
		double x,y;
		minx=atof(min_x->Text.c_str());
		maxx=atof(max_x->Text.c_str());
		miny=atof(min_y->Text.c_str());
		maxy=atof(max_y->Text.c_str());
		xmod=(Image1->Width-1)/(maxx-minx);
		ymod=(Image1->Height-1)/(maxy-miny);
		y=miny;                           // draw scale
		if((minx<=0) && (maxx>=0))
		{
			while(y <= maxy)
			{
				pixw(0,y);
				y+=yscal;
			}
		}
		x=minx;
		if((miny<=0) && (maxy>=0))
			{
			while(x <= maxx)
				{
					pixw(x,0);
 					x+=xscal;
				}
			}

}

double sqr (double InVal)
// calculates and returns the square of the input value ( InVal*InVal)
{
	return InVal*InVal;
}
//---------------------------------------------------------------------------
//****MODULE FOR DATA FILES***********************************************
int TForm1::Files()
// This module will open and prepare a file for the output of data.
// The file will in the end contain the results of the calculation.
{

	outFile.open(outfile1->Text.c_str(),(ios::ate));
	if(outFile.fail())
	// If unable to open file .... report error and quit.
	{
		Application->MessageBox("The file was not successfully opened","Error",MB_OK);

	}
	// ....Otherwise set output file stream formats
	outFile << setiosflags(ios::fixed)
		    << setiosflags(ios::showpoint)
		    << setprecision(5);
	return 0;
}

//************MODULE FOR COULOMB POTENTIAL V[i]*****************************
float* TForm1::POTENTIAL(double k,int j)
// Module will calculate the coulomb potential
{
	int i, // an index counter
	    m; // counter to keep track of maximum iterations.
	double R=0;

	m=j+2;  // m is initialised before the pointers to the array are used.
		// This is to avoid warning errors.
	float *V= new float[m];
	V[0] = 0;
	for(i=1;i<j;i++)
	{
		R=R+k;
        V[i] = -ESQ/R ;
		if(fabs(V[i])>1e10)
			break;
		else if(graphtype==1) pixw(R,V[i]);

	}
	return (V);

}
//*******************ITERATION VALUES***************************************
double TForm1::ITERATE(double sign,double E1, double DELTA_E,int* m)
{
 double Et = E1;

		if((sign) > 0)
        {
        	Et=Et+(DELTA_E);
            (*m)++;

        }
        else if((sign) < 0)
        {
           Et=Et-(DELTA_E/2);
           (*m)++;

		}


	return (Et);
}

//*******************TEST EIGENVALUE****************************************
float TForm1::EIGENVALUE(double E,double DELTA_E,double k,int j)
{
//A trial energy may be chosen which lies below ground state
 //u(r) = R(r)/r
	float *V = POTENTIAL(k,j);
 	int i;
 	double sign,TEST0,R,E1,Tol;

    double *u = new double [j+2];
	#define KSQ sqr(k)
 	#define B (KSQ/12)

	if(E < 0.000)
// if E is less than 0.000 (Zero)

// A, C & D can only be defined after all values have been received and
//  initial value set.
 		#define A ((2*MHBARSQ)*(V[i]-E))
    	#define C ((2*MHBARSQ)*(V[i+1]-E))
    	#define D ((2*MHBARSQ)*(V[i-1]-E))
// initalise values for u[0] and u[1]... prior to calculating u[i+1]

        u[0]=0;
		u[1]=1; R=k;

		for(i=1;i<j;i++)
        {
        	R=R+k;
            u[i+1]=(((2*(1+(5*B*A)))*u[i]-(1-(B*D))*u[i-1])/(1-C*B));

        }


		TEST0 = u[i];
        E1 = E + DELTA_E;

        #define F ((2*MHBARSQ)*(V[i]-E1))
        #define G ((2*MHBARSQ)*(V[i+1]-E1))
        #define H ((2*MHBARSQ)*(V[i-1]-E1))

        cout << "Please enter the tolerance level for the eigenvalue iteration:\n";
        cin >> Tol;
        Tol=_atold(tolerance->Text.c_str());
        outFile << "This is the tolerance level of the iteration:\t"<<Tol<<endl;
        int m=0;
        double e[3];
        double u_prev,u_curr=TEST0;
        double simpson_1,simpson_2=0,simpson_3;
        int ei=0; e[0]=E1; e[1]=0; e[2]=0;
        u[0]=0;
        u[1]=1; R=k;

        do {
        u_prev=u_curr;
        simpson_1=u[0]*u[0];
          for(i=1;i<j;i++)
        {
        	R=R+k;
            u[i+1]= (((2*(1+(5*B*F)))*u[i] -(1-(B*H))*u[i-1])/(1-B*G));
		}

        u_curr=u[i];
		sign=u_curr*u_prev;
        ei++; if(ei==3) { e[0]=e[1]; e[1]=e[2]; ei=2;}
		e[ei]=ITERATE(sign,e[ei-1],DELTA_E,&m);
        if(ei==2)
        {
        if(e[1]>e[0] && e[1]>e[2])
        DELTA_E=-(DELTA_E/2);
        if(e[1]<e[0] && e[1]<e[2])
        DELTA_E=-(DELTA_E/2);
        }
        E1=e[ei];
        } while(fabs(e[2]-e[1])>Tol);

        #define M ((2*MHBARSQ)*(V[i+1]-E1))
        outFile<<"C(r):\t"<<M<<endl;


		for(i=1;i<j-1;i++)
        {
        	if(i%2==0) simpson_2=simpson_2+2*u[i];
            else simpson_2=simpson_2+4*u[i];
        }
        simpson_3=(k/3)*(simpson_1+simpson_2+u[i]*u[i]);
        simpson_3=sqrtl(fabs(simpson_3));
        outFile	<<setiosflags(ios::left)<<"r/nm:"
                <<setiosflags(ios::right)<<setw(60)<<"C(r)u[r+1]_normalized:"
				<<setiosflags(ios::right)<<setw(150)<<"u[r+1]_normalized:"<<endl;
        R=0;
        for(i=0;i<j;i++)
        {
           R=R+k;
           u[i]=u[i]/simpson_3;
           if(fabs(u[i])>1e10)
				break;
				if(graphtype==0)
				pixw(R,u[i]);
           outFile<<setiosflags(ios::left)<<R;
           outFile<<setiosflags(ios::right)<<setw(70)<<M*u[i];
           outFile<<setiosflags(ios::right)<<setw(150)<<u[i]<<endl;
        }

    // write out the column headers in to the file.

    outFile <<"Final DELTA_E:"<<DELTA_E<<endl;
	outFile <<"This is the counter index of the iteration of the step size of the energy:\t"
            <<m<<endl;
    outFile <<"The E eigenvalue:\n"<<ITERATE(sign,E1,DELTA_E,&m)<<endl;


	//otherwise report error
     if (E>0)
	{
       	Application->MessageBox("Value that you have entered is invalid! Please try again.\n","Error",MB_OK);
       	exit (1);
    }
    else
    {
    Application->MessageBox("Well, your turn for fun!\n","The End",MB_YESNO);
    }
    //close the file
	outFile.close();
    delete u;
	return 0;
 }

//*********MODULE FOR WAVEFUNCTION GENERATION*******************************
TForm1::main()
// To Generate a wavefunction for any atom given any Eigenvalue.
{
 	int   j;              //angular momentum's variable and radius variable

   double 	k,            // number of step size to be generated
    		E,            // Trial value of Energy
            DELTA_E;
//-----------------------------------------------------------------
 cout << setiosflags (ios::fixed);
 cout << setiosflags (ios::scientific);
    // Start by preparing the Output File using the Files function.
	Files();
	k=_atold(Step_size->Text.c_str());
	j=atoi(num_step_size->Text.c_str());
    E=atof(tr_energy->Text.c_str());
    DELTA_E=atof(Step_energy->Text.c_str());

// Write out the information just read from the user in to the file.
   outFile << "This is the initial trial energy entered:\t"<<E<<endl;
    outFile <<"Step size chosen:\t"<<k<<endl;
    outFile << "Number of step size generated :\t"<<j<<endl;
    outFile << "The delta E chosen:\t"<<DELTA_E<<endl;

// Do actual processing
    EIGENVALUE(E,DELTA_E,k,j);

    return 0;
}
//********End Of File***********************************************************

//---------------------------------------------------------------------------
void __fastcall TForm1::onclick(TObject *Sender)
{
	Image1->Canvas->FillRect(Image1->Canvas->ClipRect);
	initgraph(atof(x_scal->Text.c_str()),atof(y_scal->Text.c_str()));
	graphtype=0;
	y_label->Caption="u[i+1]";
	main();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::v_coul(TObject *Sender)
{
	Image1->Canvas->FillRect(Image1->Canvas->ClipRect);
	initgraph(atof(x_scal->Text.c_str()),atof(y_scal->Text.c_str()));
	graphtype=1;
	y_label->Caption="V[i+1]/V";
	main();
}