#include <math.h>
#include <vector>
#include <stdio>
#include <iomanip>

double c0 = 1;
double e = 1;
double pi = 3.14;
double ep_0 = 1;
double hbar = 1;
double m = 1;

double a = 4*pi*ep_0*hbar*hbar/m/e/e;

// hydrogen atom model

// kappa = sqrt(-2*m*E/hbar)
// kappa = m * e^2 / (4 * pi * ep_0 * hbar^2)

int factorial(int n)
{
    return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
}

double kappa(double n) 
{
    double result = m * e^2 / (4 * pi * ep_0 * hbar^2);
    result /= n;
    return result;
}

double c_j (double j, double n, double l) 
{
  double result = (2 *(j+l+1-n)/(j+1)/(j+2*l+2));
  result *= c_j(j);
} 

double v(double n, double l, double rho)
{
    double jmax = n-l-1;
    std::vector<double> c;
    for (double power=0; power < jmax + 1; power++) {
        c[power] = c0
    }
}

double R(double n, double l,
            double r)
{
    double rho = r*a/n;
    double result = 1/r 
        * pow(rho,l+1) 
        * exp(-rho)
        * v(n,l,rho);
}

double Y(double l, double m, 
            double theta, double phi);

double psi(double n, double l, double m,
            double r, double theta, double phi)
{
    norm(n,l) * R(n,l,r) * Y(l,m,theta,phi);
}
double norm(double n, double l) 
{
    double result_squared = pow((2*a/n),3) * factorial(n-l-1) / (2*n) / pow(factorial(n+1),3);
    return sqrt(result_squared);
}


int main(int argc, char const *argv[])
{
    // nanometers
    for (double r=0; r<100; r += 0.005)
    std::cout
        << std::scientific << std::setprecision(6)
        << r << '\t'
        << R(1,0,r) << '\t'
        << R(2,0,r) << '\t'
        << R(3,0,r) << '\t'
        << R(4,0,r) << '\t'
        << R(5,0,r) << '\t'
        << R(6,0,r) << '\t'
        << R(7,0,r) << '\t'
        << R(8,0,r) << '\t'
        << R(9,0,r) << '\t'
        << R(10,0,r) << '\t'
        << R(11,0,r) << '\t'
        << R(12,0,r) << '\t'
        << R(13,0,r) << '\t'
        << R(14,0,r) << '\t'
        << R(15,0,r);


    //for (double r=0; r<100; r += 0.005)
    //std::cout
    //    << std::scientific << std::setprecision(6)
    //    << r << '\t' 
    //    << psi(1,0,0,r,0,0) << '\t'
    //    << psi(2,0,0,r,0,0) << '\t'
    //    << psi(3,0,0,r,0,0) << '\t'
    //    << psi(4,0,0,r,0,0) << '\t'
    //    << psi(5,0,0,r,0,0) << '\t'
    //    << psi(6,0,0,r,0,0) << '\t'
    //    << psi(7,0,0,r,0,0) << '\t'
    //    << psi(8,0,0,r,0,0) << '\t'
    //    << psi(9,0,0,r,0,0) << '\t'
    //    << psi(10,0,0,r,0,0) << '\t'
    //    << psi(11,0,0,r,0,0) << '\t'
    //    << psi(12,0,0,r,0,0) << '\t'
    //    << psi(13,0,0,r,0,0) << '\t'
    //    << psi(14,0,0,r,0,0) << '\t'
    //    << psi(15,0,0,r,0,0);
    return 0;
}