home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / pascal / library / dos / nrpas / gauleg.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1994-04-11  |  959 b   |  35 lines
  1. PROCEDURE gauleg(x1,x2: double; VAR x,w: darray; n: integer);
  2. (* Programs using routine GAULEG must define the type
  3. TYPE
  4.    darray=ARRAY [1..n] OF double;
  5. in the calling program *)
  6. CONST
  7.    eps = 3.0e-11; (* adjust to your floating precision *)
  8. VAR
  9.        m,j,i: integer;
  10.        z1,z,xm,xl,pp,p3,p2,p1: double;
  11. BEGIN
  12.    m := (n+1) DIV 2;
  13.    xm := 0.5*(x2+x1);
  14.    xl := 0.5*(x2-x1);
  15.    FOR i := 1 TO m DO BEGIN
  16.       z := cos(3.141592654*(i-0.25)/(n+0.5));
  17.       REPEAT
  18.              p1 := 1.0;
  19.              p2 := 0.0;
  20.              FOR j := 1 TO n DO BEGIN
  21.                p3 := p2;
  22.                p2 := p1;
  23.                p1 := ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
  24.              END;
  25.              pp := n*(z*p1-p2)/(z*z-1.0);
  26.              z1 := z;
  27.              z := z1-p1/pp;
  28.       UNTIL (abs(z-z1) <= eps);
  29.       x[i] := xm-xl*z;
  30.       x[n+1-i] := xm+xl*z;
  31.       w[i] := 2.0*xl/((1.0-z*z)*pp*pp);
  32.       w[n+1-i] := w[i]
  33.    END
  34. END;
  35.