``getphidebug.c''

#include <stdio.h>
#include <math.h>
#include "p480.h"
#include "nrutil.h"

#define Nmx 40000
main()
{
  /* Change of variables info */
  double *r,*dr;
  int N;

  /* Physics variables */
  double *V,*Rho,*phi,*integrand,Z;
  int lmax,*nmax,nmaxmax;
  double **E,***Psi,**F;

  /* Working variables */
  int n,l,k;
  double x;

  /* Value of pi */
  const double pi=4.*atan(1.);

  /* Specs for hydrogen */
  Z=1.; /* Charge on nucleus */
  lmax=0; /* Only up to 's' states (l=0) are filled in H */

  nmax=ivector(0,lmax); /* Declare space for max nodes for each l */
  nmax[0]=0; /* Specify that only up to the first (zero node)
                's' states are filled in H */

  /* General loop for finding nmaxmax (needed to define arrays)...
     Note: This is a little silly for H, but is handy when you have
     some random element. */
  nmaxmax=0;
  for (l=0; l<=lmax; l++)
    if (nmax[l]>nmaxmax) nmaxmax=nmax[l];

  F=dmatrix(0,lmax,0,nmaxmax); /* Allocate space for filling data */
  F[0][0]=1.; /* For H, there is 1 electron in the '1s' (l=0,n=0) state */

  /* Allocate space for table for E's */
  E=dmatrix(0,lmax,0,nmaxmax);

  for (N=40; N<=Nmx; N*=10) {
    /* Allocate space for N dependent physics arrays */
    Psi=d3tensor(0,lmax,0,nmaxmax,0,N);
    Rho=dvector(0,N);
    phi=dvector(0,N);
    
    integrand=dvector(0,N);
    
    /* Space for grid vectors */
    r=dvector(0,N); 
    dr=dvector(0,N);
    V=dvector(0,N);

    /* Set up grid */
    for (k=0; k<=N; k++) {
      x=((double) k)/((double) N);
      r[k]=1./(1.-x)-1.-x-x*x;
      dr[k]=1./(1.-x)/(1.-x)-1.-2*x;
      V[k]=-Z/r[k];
    }
    V[0]=0.; dr[N]=0.;

    /* Test section */  
    getallEs(E,lmax,nmax,Z,V,r,dr,N);
    getallPsis(Psi,E,lmax,nmax,V,r,dr,N);
    getRho(Rho,Psi,F,lmax,nmax,N);
    getphi(phi,Rho,r,dr,N);

    for (k=0; k<=N; k++)
      integrand[k]=phi[k]*Rho[k];
    printf("%6d %20.12f (Hartree Energy)\n",N,simpint(integrand,r,dr,N)/2);

    /* Clean up N-dependent vectors BEFORE next value of N! */
    free_dvector(r,0,Nmx);
    free_dvector(dr,0,Nmx);
    free_dvector(V,0,Nmx);
    
    free_d3tensor(Psi,0,lmax,0,nmaxmax,0,Nmx);
    free_dvector(Rho,0,Nmx);
    free_dvector(phi,0,Nmx);
    free_dvector(integrand,0,N);
  }

  /* Free up all remaining objects */
  free_dmatrix(E,0,lmax,0,nmaxmax);
  free_ivector(nmax,0,lmax); 
}