// solve time independent Schrodinger equation
#include <iostream>
#include <cmath>

using namespace std;
/*     Difinition of function g(x,y,E)    */
double g(double r, double A, double E, double L)
{
        const double hbar=197.327 ;
        const double m=510.999e3;
        const double alpha=1./137.036;
        const double c1=-2.*m/(hbar*hbar);
        const double c2=alpha*hbar;
        return c1*(E+c2/r)*A+L*(L+1.)*A/(r*r);
}
// x,y,z were substituted by r,A,B. B wa A no r ni yoru bibun.
void runge( double *r, double *A, double *B, double E, double h, double L)
{
        double k1, k2, k3, k4, l1, l2, l3, l4;

        k1 = h*(*B);        l1= h*g(*r, *A, E, L);
        k2 = h*(*B+l1/2.0); l2= h*g(*r+h/2.0, *A+k1/2.0, E, L);
        k3 = h*(*B+l2/2.0); l3= h*g(*r+h/2.0, *A+k2/2.0, E, L);
        k4 = h*(*B+l3);     l4= h*g(*r+h, *A+k3, E, L);
        *A += k1/6.0 + k2/3.0 + k3/3.0 + k4/6.0;
        *B += l1/6.0 + l2/3.0 + l3/3.0 + l4/6.0;
        *r += h;
        //cout <<"k's: "<< k1 <<" " << k2 << " " << k3 << " " << k4 << endl;
        //cout <<"l's: "<< l1 <<" " << l2 << " " << l3 << " " << l4 << endl;
}

int main(int argc, char *argv[])
{
	// d2A/dr2=-2*m/hbar^2[E+alpha*hbar*c/r]A, A=rR ,R:radial function
  // dA/dr=B
  // time independent solution using Runge Kutta method
  //const double hbar=197. ; 
  //const double m=511e3;
  //const double alpha=1./137.;
  if(argc!=3) {
    cout << "Usage: ./Hydrogen n L \n"<<"./Hydorogen 2 1 (for n l \n";
    return 1;
  }
  double n =atof(argv[1]);
  double L =atof(argv[2]);
  //cout << "input Energy: " << E << endl;
  const double h=0.0001;
  double A=0.01;
  double B=10.;
  double r=0.001;
  double E=-13.6057/(n*n);
  //cout << "input Energy: ";
  //cin >> E;
  for ( int i=0; i < 25000; i++){
    runge( &r, &A, &B, E, h , L);
    if( i%100 == 0)  cout << r << " " << A << endl;
  }
  
  return 0;
}