/*        ATK IV - Numerical Programming (2006)       */
/*    Excercise 11.1 -  Spectrum of a hydrogen atom   */

#include<stdio.h>
#include<math.h>
#define NRANSI
#include<nr.h>
#include<nrutil.h>

#define maxn 60
#define Pi 3.141926

int l;    /*angular momentum*/
float energy;
float r0,rinf,rm;

float ivzero[2+1], ivinf[2+1];   /*initial values*/
float **u, *r;
int imax;

/*These are for odeint, not used otherwise*/
int kmax = 0, kount;
float *xp,**yp,dxsav;

float delta(float);

float normsqr(void);

void integrate(float *, float *, float, float);

void schroedinger(float,float[],float[]);

int main(void){
  float ul[2+1], ur[2+1];  /* Used to get the final solution*/
  float energies[maxn];
  float step,dold,dnew,newe,cnorm,h,p;

  int ne,ifit;
  int i,j,n;

  FILE *file;

  file=fopen("print.txt","w");

  l = 1;
  r0 = 1.0e-4;
  rm = 4.0;
  rinf = 100.0;

  ivzero[1] = 0.0;
  ivzero[2] = 1e6;
  ivinf[1] = 0.0;
  ivinf[2] = 1e-10;

  imax = 1001;
  ifit = (imax-1)*((rm-r0)/(rinf-r0))+1;
  u = matrix(1,2,1,imax);
  r = vector(1,imax);

  n = maxn;
  ne = 0; /* number of eigenvalues found */

  step = 1.0/n;
  energy = -1.0; /* All eigenvalues are found from the interval -1 to 0 */

  dold = delta(energy);

  for(i=1;i<=n;i++){
    energy = i*step -1.0;

    dnew = delta(energy);

    if (dnew*dold < 0.0){
      newe = rtbis(delta,energy-step,energy,1e-6);

      /* To obtain true roots we must check...*/
      if(fabs(delta(newe))<1e-3){
	ne++;
	energies[ne]=newe;
      }
    }

    dold = dnew;
  }

  for(i=1;i<=ne;i++){
    printf("%f %f\n",energies[i], -1/SQR(i+1));
  }
  for(i=1;i<=ne;i++){
    energy = energies[i];

    integrate(ul,ivzero,r0,rm);
    integrate(ur,ivinf,rinf,rm);

    ivinf[2]*=ul[1]/ur[1];  /* u(r) and its derivative continuous at rm */
    
    n = ifit;
    p = r0;
    h = (rm-r0)/(n);

    ul[1] = ivzero[1];
    ul[2] = ivzero[2];

    for(j=1;j<=n;j++){
      r[j] = p;
      u[1][j] = ul[1];
      u[2][j] = ul[2];

      integrate(ul,ul,p,p+h);
      p = h*j + r0;
    }

    n = imax - ifit;
    p = rinf;
    h = (rm - rinf)/(n-1);
    ur[1] = ivinf[1];
    ur[2] = ivinf[2];

    for(j=1;j<=n;j++){
      r[imax-j+1] = p;
      u[1][imax-j+1] = ur[1];
      u[2][imax-j+1] = ur[2];
      integrate(ur,ur,p,p+h);
      p = h*j + rinf;
    }

    cnorm = 1.0/sqrt(4.0*Pi*normsqr()/3.0);

    /*Normalization*/

    for(j=1;j<=imax;j++){
      u[1][j]*=cnorm;
      u[2][j]*=cnorm;
      }
   
    for(j=1;j<=imax;j++){
      fprintf(file,"%12.7f %12.7f %12.7f\n",r[j],SQR(u[1][j]/r[j]),u[2][j]);
    }
    fprintf(file,"\n\n");
  }

  fclose(file);
  free_vector(r,1,imax);
  free_matrix(u,1,2,1,imax);

  return 0;
}

float delta(float e){
  float ul[2+1];
  float ur[2+1];
  float logdl,logdr;

  energy = e;

  integrate(ul,ivzero,r0,rm);
  integrate(ur,ivinf,rinf,rm);

  logdl = ul[2]/ul[1];
  logdr = ur[2]/ur[1];

  return logdl - logdr;
}

void integrate(float *u, float *ivs, float a, float b){
  int ok,bad;
  float h;

  h = fabs(b-a);

  u[1] = ivs[1];
  u[2] = ivs[2];

  odeint(u,2,a,b,(float)1e-6,h,0.,&ok,&bad,schroedinger,rkqs);
}

float normsqr(void){
  int i;

  float s=0;

  float y0[2+1],y1[2+1],h;

  y0[1] = SQR(u[1][1]);
  y0[2] = 2*u[1][1]*u[2][1];

  for(i=2;i<=imax;i++){
    h = -r[i-1]+r[i];
    y1[1] = SQR(u[1][i]);
    y1[2] = 2*u[1][i]*u[2][i];

    s+=h*h*(y0[2]-y1[2])/12.0 + h*(y0[1]+y1[1])/2.0;

    y0[1] = y1[1];
    y0[2] = y1[2];

  }
  return s;
}

void schroedinger(float r, float u[], float du[]){
  du[1] = u[2];
  du[2] = (-energy -2.0/r + l*(l+1)/(r*r))*u[1];
}