#include "param.h"

/********************************************/
/*  Integration der Schroedinger-Gleichung  */
/*  Rueckgabewert: Energie-Eigenwert        */
/********************************************/

double schrod (double dr,     /* Schrittweite                      */
               int nr,        /* Anzahl Stuetzstellen              */
               double pot[],  /* Potential                         */
               double wf[],   /* Wellenfunktion                    */
               double ewert,  /* Startwert fuer Eigenwert-Suche    */
               double de,     /* Schrittweite auf Energieskala     */
               double egenau, /* Genauigkeit des Eigenwertes       */
               double emax)   /* Obergrennze: sinnvoller Eigenwert */
{
  double wfnorm, wf_o, v1, v2, v3, vfak, ewerto;
  int n;
  int anfang = 1;  /* Switch: Erste Iteration   */
  int sekant = 0;  /* Switch: Sekantenverfahren */
  int iter   = 0;  /* Zaehler fuer Iterationen  */
  ewerto = ewert;
  ewert  = ewert + de;

  while (ewert != ewerto && fabs (de) >= egenau) {
    if (ewert > emax) {
      fprintf (stderr, "ewert > emax:");
      exit (1);
    };
    iter ++;
    ewerto = ewert;

    /* Randbedingung am Rand */
    wf [nr]   = 0.0;
    wf [nr-1] = 1.0e-25;

    vfak = dr * dr / 6.0;
    v1 = vfak * (ewert - pot [nr]);
    v2 = vfak * (ewert - pot [nr-1]);

    /* Integration der Schroedingergleichung: Numerow-Algorithmus */
    for (n = nr-2; n>=0; n--) {
      v3 = vfak * (ewert - pot [n]);
      wf [n] = (  (2.0 - 10.0 * v2) * wf [n+1]
                - (1.0        + v2) * wf [n+2]) / (1.0 + v3);
      v1 = v2;
      v2 = v3;
    }

    /* Bei Vorzeichenwechsel umschalten auf Sekantenmethode */
    if (! (sekant || anfang)) sekant = wf [0] * wf_o < 0.0;
    anfang = 0;
    /* Sekantenmethode */
    if (sekant) de = - wf [0] * de / (wf [0] - wf_o);
    ewert = ewert + de;
    wf_o = wf [0];
  }

  printf ("Eigenwert = %f  de = %f,  %i Iterationen\n", ewert, de, iter);

  /* Normierung */
  wfnorm = 0.0;
  for (n = 0; n <= nr; n++) wfnorm += wf [n] * wf [n];
  wfnorm = 1.0 / sqrt (wfnorm * dr);
  for (n = 0; n <= nr; n++) wf [n] *= wfnorm;

  return ewert;
}