Avoiding overflow in ``schint''

 

As we discussed in lab, we can exploit the linearity of the Schrödinger equation to avoid overflow of Psi and Psip in ``schint''. In particular, you should add the ``if'' statement in the code fragment below to your ``schint'' routine.

#define SMALL (1.e-300)
#define SQRTSMALL (1.e-150)
#define BIG (1.e300)
#define SQRTBIG (1.e150)
#define TOL (1.e-12)

int schint(double *Psi, double *Psip, double Psiout[],
	   int k1, int k2, double Veff[], double E,
	   double r[], double dr[], int N)
{
                             .
                             .
    rk4p480(y,dydx,2,x,h,y,derivs_Schrodinger,k,dk,T,r,dr,N);
    if (fabs(*Psi)>SQRTBIG) {
      y[1]=SMALL*y[1];
      y[2]=SMALL*y[2];
    }
                             .
                             .

Note that above, we ``#define'' the various constants (``parameter'' in Fortran). To avoid debugging problems, such constants should never appear explicitly in your code. Rather, as above, they should be defined in a way which is easily seen and easily modified.