//  file:shooting_method_normalization.cpp
//
//  This program solves the schroedinger equation for an infinite square well using the shooting method
//
//  Programmer:  XXXXX
//
//
//  Notes:
//   * compile with: the corresponding make file, make_shooting_method_normalization
//

//*********************************************************************//

// include files
#include <iostream> // this has the cout, cin definitions
#include <cmath> // useful math functions
#include <fstream> // allows me to stream things into data files
#include <iomanip> // use things such as setprecision and setw
using namespace std; // if omitted, then need std::cout, std::cin

//*********************************************************************//


int
main ()
{
//first define the variables needed for finding the wave function and normalizing it.
  double L = 1.0;                  //length of well
  double hc = 2.0;
  double psi = 0.0;                //guess wavefunction
  double dpsi = 1.0;               //first derivative of the wavefunction
  double ddpsi = 0.0;              //second derviative of the wavefunction
  double E = 2.40;                 //Energy found through shooting method
  double x = 0;
  double dx = .01;
  double V = 0;                     //potential energy
  double xt = 0;
  double b = 4.37;                  //test normalization factor for numerical method. Change around to find correct factor.
  double wavefunct = 0;             //known wavefunction
  double normalize = 0;
  double A = 0;
  double pi=atan(1.0)*4.0;
  double psinormal=0;

//will first integrate the analytical solution to find a normalized probability density
  ofstream out2;
out2.open ("analytical_solution.dat");
out2 << "x*t       wavefunction^2"<<endl;
for (double j=0.0; j<=L; j +=.01)
{
xt=0;
wavefunct=sqrt(2/L)*sin(pi*j/L);
normalize=normalize+(wavefunct*wavefunct*dx);
xt=xt+j;

out2<<fixed<<setprecision(8)<<setw(5)<<xt<< " "
    <<fixed<<setprecision(8)<<setw(5)<<wavefunct*wavefunct<< " "
     <<endl;
}
out2.close();

cout<<" The integrated probability density of the analytical function is:"<<normalize<<endl;




//now will integrate my numerical algorithm to find to normalize it.
ofstream out;
out.open ("numerical_solution.dat");
out<< " x          psinormal^2 "<<endl;
for (double i=0.0; i<=L; i +=.01)
{
 x=0;                      //for reinitialization
 ddpsi=-2*hc*(E-V)*psi;    //second derivative of the wave function
 dpsi=dpsi+ddpsi*dx;        //first derivative of the wave function
 psi=psi+dpsi*dx;           //wave function
 x=x+i;                    //calculates the x value fo the particle from 0 to L

 psinormal=b*psi;

 A=A+psinormal*psinormal*dx;


out<<fixed<<setprecision(8)<<setw(5)<<x<< " "
   <<fixed<<setprecision(8)<<setw(5)<<psinormal*psinormal<< " "
   <<endl;
}