``test_getPsi.c''

 

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

  /* Physics variables */
  double *V;
  double Z;
  int lmax;
  int nmaxmax,*nmax;
  double **E,*psi;

  int n,l;

  /* Working variables */ 
  double x,*integrand;
  int k;

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

  /* Specs for hydrogen */
  Z=1.;
  lmax=0;

  nmax=ivector(0,lmax);
  nmax[0]=0;

  nmaxmax=0;
  for (l=0; l<=lmax; l++)
    if (nmax[l]>nmaxmax) nmaxmax=nmax[l];

  /* Allocate remaining vectors */

  E=dmatrix(0,lmax,0,nmaxmax); /* E's */
  psi=dvector(0,Nmx); /* psi */

  r=dvector(0,Nmx); /* Grid vectors */
  dr=dvector(0,Nmx);
  V=dvector(0,Nmx);

  integrand=dvector(0,Nmx); /* Working grid for integrating things */

  /* Convenient loop over N */
  for (N=40; N<=Nmx; N*=10) {

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

    getallEs(E,lmax,nmax,Z,V,r,dr,N); /* Get E */
    getPsi(psi,E[0][0],0,V,r,dr,N); /* Get Psi for E[l=0][n=0] */

    for (k=0; k<=N-1; k++) /* Get integrand for root mean square (rms) error */
      integrand[k]=(psi[k]-2*r[k]*exp(-r[k]))*(psi[k]-2*r[k]*exp(-r[k]));

    printf("%6d %20.16f\n",N,sqrt(simpint(integrand,r,dr,N)));

  }
  free_dvector(r,0,Nmx);
  free_dvector(dr,0,Nmx);
  free_dvector(V,0,Nmx);

  free_dmatrix(E,0,lmax,0,nmaxmax);
  free_dvector(psi,0,Nmx);
  free_dvector(integrand,0,Nmx);
}