/*
 *
 * 	Solution of the Schrodinger equation for harmonic oscillator
 * 	Numerov algorithm; forward & backward integration; 
 * 	eigenvalue search using the shooting method.
 *
 *  dimensionless units: x = (m k/hbar^2)^(1/4) X
 *                       e = E/(hbar omega)
 *
 *  where omega = sqrt(k/m) - frequency of the oscillator,
 *  m - its mass, k - ellastic constant (V(X) = 1/2 k X^2)
 *
 *  For the explanation of the code see:
 *  http://www.fisica.uniud.it/~giannozz/Corsi/MQ/LectureNotes/mq-cap1.pdf
 */

#include <stdio.h>
#include <math.h>

#define MSHX 2000		/* Max number of grid points */
#define EERR 1.e-10		/* Error in finding eigenvalues */

int main ()
{
  int mesh, i, icl = 0;
  int nodes, hnodes = 0, ncross, parity, iter;
  double xmax, dx, ddx12, norm;
  double elw, eup, e = 0.;
  double x[MSHX+1], y[MSHX+1], vpot[MSHX+1], f[MSHX+1];

  /*  Read input data */

  fprintf (stderr, "Max value of x (typical value: 10) ? ");
  scanf ("%lf", &xmax);
  fprintf (stderr, "Number of grid points (max:%5d) ? ", MSHX);
  scanf ("%d", &mesh);
  if (mesh > MSHX) {
      mesh = MSHX;
  }

  /*  Initialize grid  */

  dx = xmax/mesh;
  ddx12 = dx*dx/12.;

  /*  set up the potential  */

  for (i = 0; i <= mesh; ++i) {
    x[i] = i*dx;
    vpot[i] = x[i]*x[i]/2;
  }

  /* eigenvalue search */
  while (1 == 1) {

    /* Read number of nodes (stop if < 0) */
    fprintf (stderr, "Number of nodes (-1=exit) ? ");
    scanf ("%d", &nodes);
    if (nodes < 0) {
	  return (0);
	}

    hnodes = nodes/2;

    /*  
     *  if nodes is even, there are 2*hnodes nodes
     *  if nodes is odd,  there are 2*hnodes+1 nodes (one is in x=0)
     *  hnodes is thus the number of nodes in the x>0 semi-axis (x=0 excepted) 
     */

    if (2*hnodes == nodes)
  	  parity = +1;
    else
  	  parity = -1;

    /* 
     * set initial lower and upper bounds of the eigenvalue 
     * to the smallest and the largest valuest of the potential respectively
     */

    eup = vpot[mesh];
    elw = eup;
    for (i = 0; i <= mesh; i++) {
	  if (elw > vpot[i])
	    elw = vpot[i];
	  if (eup < vpot[i])
	    eup = vpot[i];
	}

    iter = 0;
	while (eup - elw > EERR) {

      /* search eigenvalues using bisection */
	  e = 0.5 * (elw + eup);

      iter += 1;

      /*
       *   set up the f-function used by the Numerov algorithm
       *   and determine the position of its last crossing, i.e. change of sign
       *   f < 0 means classically allowed   region
       *   f > 0 means classically forbidden region
       */
      f[0] = 2*ddx12*(vpot[0] - e);
      icl = -1;

      for (i = 1; i <= mesh; ++i) {
	    f[i] = 2*ddx12*(vpot[i] - e);
        /*
         *   if f[i] is exactly zero the change of sign is not detected.
         *   the following line is a trick to prevent missing a change of sign 
         *   in this unlikely but not impossible case:
         */
	    if (f[i] == 0.) {
	      f[i] = 1e-20;
	    }
        /* store the index 'icl' where the last change of sign has been found */
	    if (f[i] != copysign(f[i], f[i-1])) {
	      icl = i;
	    }
	  }

      /* check the validity of xmax choice */

      if (icl >= mesh - 2) {
	    fprintf (stderr, "last change of sign too far.");
	    return (1);
	  }
      if (icl < 1) {
	    fprintf (stderr, "no classical turning point?");
	    return (1);
	  }

      /* f(x) as required by the Numerov algorithm  */
      for (i = 0; i <= mesh; ++i) {
	    f[i] = 1. - f[i];
	  }

      /* Determination of the wave-function in the first two points  */

      if (parity == 1) {
	    /*  even number of nodes: wavefunction is even  */
	    y[0] = 1.;
	    /*  assume f(-1) = f(1)  */
	    y[1] = 0.5*(12. - 10*f[0])/f[1];
	  } else {
	    /*  odd number of nodes: wavefunction is odd  */
	    y[0] = 0.;
	    y[1] = 1.;
	  }

      /*  Outward integration */
      ncross = 0;
      for (i = 1; i < mesh; i++) {
	    y[i+1] = ((12. - 10*f[i])*y[i] - f[i-1]*y[i-1])/f[i+1];
	    if (y[i] != copysign (y[i], y[i+1]))
	      ++ncross;
	  }

      /*  Check number of crossings  */

      fprintf (stderr, "%4d %4d %14.8f %14.8f %14.8f\n", iter, ncross, elw, eup, e);
	  if (ncross > hnodes) {
        /* Too many crossings: current energy is too high -> lower the upper bound */
	    eup = e;
      } else {
	    /* Too few or correct number of crossings: current energy is too low -> raise the lower bound */
	      elw = e;
      }
    }
    
    /* Inward (stable) integration */

    /*
     *   Determination of the wave-function in the last two points 
     *   assuming y(mesh+1) = 0
     */      
    y[mesh] = dx; 
    y[mesh-1] = (12. - 10*f[mesh])*y[mesh]/f[mesh-1];
      
    for (i = mesh-1; i >= 0; i--) {
        y[i-1] = ((12. - 10*f[i])*y[i] - f[i+1]*y[i+1])/f[i-1];
    }   

    /* Wavefunction normalization */

    norm = 0.;
    for (i = 1; i <= mesh; ++i) {
        norm += y[i]*y[i];
    }
    norm = dx * (2*norm + y[0]*y[0]);
    norm = sqrt(norm);
    for (i = 0; i <= mesh; ++i) {
        y[i] /= norm;
    }

    /*
     * output calculated wavefunction
     */

    /* x<0 region: */
    for (i = mesh; i >= 1; --i) {
	  fprintf (stdout, "%7.3f %16.8e %16.8e %12.6f\n",
		   -x[i], parity * y[i], y[i] * y[i], vpot[i]);
	}

    /* x>0 region: */
    for (i = 0; i <= mesh; ++i) {
	  fprintf (stdout, "%7.3f %16.8e %16.8e %12.6f\n",
		   x[i], y[i], y[i] * y[i], vpot[i]);
	}
    fprintf (stdout, "\n\n");
  }

  return(0);
}