getallEs()

The routine getEs() from the previous problem set finds all wave functions for a given one-dimensional potential. As discussed in class, we can generate all of the solutions (spherically symmetric or not) to a given radial potential by solving for the one-dimensional states in the effective potential

$$V_{eff}(r) \equiv V(r)+\frac{\ell(\ell+1)}{2r^2},$$

where $\ell$ ($\ell$=0, 1, 2, ...) is the angular momentum of the given state. For atoms, we know how many states we need for each angular momentum. Thus to find all the states, all we need do is to call getEs() once for each value of $\ell$. For convenience, we package this capability into a new routine, which you should place in ``physics.c'' along with our other physics subroutines:

void getallEs(double **E,
              int lmax, int nmax[],
              double Z, double V[],
              double r[], double dr[], int N)

Output:

• E[l][n]: matrix of energies of states of angular momentum l and n=0 nodes through n=nmax[l] nodes (inclusive!)

Input:

• lmax: maximum angular momentum to consider (l=0 through l=lmax)
• nmax[l]: maximum number of nodes sought for angular momentum l
• Z: charge on the nucleus
• V[]: ``self-consistent potential'' (V[k]=-Z/r[k], for now ...)
• r[], dr[], N: standard grid information

Write a routine of the above prototype which returns the energies of all requested electronic states in an atom with radial potential V[]. The routine should return the results in the form of a double precision two dimensional matrix of values E[l][n]. It should determine the following table of values,

l=0	E[0][0]	E[0][1]	...	E[0][nmax[0]]
l=1	E[1][0]	E[1][1]	...	E[1][nmax[1]]
l=lmax	E[lmax][0]	E[lmax][1]	...	E[0][nmax[lmax]]

If you are new to C, you should note here how two- (and higher-) dimensional arrays are accessed.

A suitable, but not the most efficient, code fragment is

  double *Veff;
                       .
                       .
  Veff=dvector(0,N); /* Must declare space for Veff */
                       .
                       .
    for (k=0; k<=N; k++) { /* Compute Veff */
      Veff[k]=V[k]+l*(l+1)/(2*r[k]*r[k]);
    }
    Veff[0]=0.; /* Mathematically proper for origin */
                       .
                       .
    getEs(E[l],nmax[l],-Z*Z,Veff,r,dr,N); /* Solve for states of ang mom l */
    /* Note: Above, passing in E[l] sends the pointer to the appropriate row
             of E[][], which getEs() just fills in appropriately as a
             one-dimensional array */
                       .
                       .
  free_dvector(Veff,0,N); /* Must free Veff at end of routine! */

Note that, here, as discussed in lecture, we use -Z*Z as a guaranteed lower bound. Note the convenient use of the Numerical Recipes dvector() routine to declare space for Veff.

Hints on working with two-dimensional matrices:

Note that the above code fragment will only work if E is properly passed in as an array of pointers, which can get quite tricky unless you have a lot of experience with dealing with pointers in C. Fortunately, the Numerical Recipes dmatrix() function can take care of all of the needed set-up for you. A typical use of dmatrix() in a program would be

{
  double **E; /* In the declarations section */
                   .                
                   .                
  E=dmatrix(0,lmax,0,nmaxmax); /* Among the first executables -
                   .              lmax is the maximum l value we expect;
                   .              nmaxmax is the maximum expected of all nmax[l]'s */
    getallEs(E,lmax,nmax,Z,V,r,dr,N); /* Example call */
    printf("E[l=0][n=0]=%f, E[l=1][n=0]=%f\n",E[0][0],E[1][0]); /* Use E[][] like this */
                   .                
                   .                
  free_dmatrix(E,0,lmax,0,nmaxmax); /* At end of routine or main, must call free
                                       for each dvector and dmatrix before exiting! */
}