#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); }