#include "param.h"

/************************************************************/
/* Loesung der Poisson-Gleichung fuer die Ladungsdichte Rho */
/************************************************************/

double pois3d (double dr,     /* Schrittweite                      */
               int nr,        /* Anzahl Stuetzstellen              */
               double rf[],   /* Stuetzstellengitter               */
               double rho[],  /* Ladungsdichte                     */
               double pot[])  /* Potential                         */
{
  int n;
  double s3, m, vh;

  /* Anfangswerte fuer die Vorwaertsintegration */
  double con = dr * dr / 12.0;
  double s1 = 0.0;
  double s2 = - 4.0 * PI * rho[1] * rf[1];
  pot[0] = 0.0;
  pot[1] = dr;

  /* Numerow-Algorithmus */
  for (n = 1; n <= nr-1; n++) {
    s3 = - 4.0 * PI * rho[n+1] * rf[n+1];
    pot[n+1] = 2.0 * pot[n] - pot[n-1] + con*(10.0 * s2 + s3 + s1);
    s1 = s2;
    s2 = s3;
  }

  /* Ein lineares Verhalten der Loesung wird subtrahiert */
  m = (pot [nr] - pot [nr-10]) / (10.0 * dr);
  vh = 0.0;
  for (n = 1; n <= nr; n++) {
    pot [n] = pot [n] / rf[n] - m;
    vh += pot[n] * rho [n] * rf [n] * rf [n];
  }
  return (4.0 * PI * dr * vh);
}