PROGRAM d4r8(input,output);
(* driver for routine GAULEG *)
CONST
   npoint=10;
   x1=0.0;
   x2=1.0;
   x3=10.0;
   nval=10;
TYPE
   double = real;
   darray=ARRAY [1..npoint] OF double;
VAR
   i,j : integer;
   xx : real;
   x,w : darray;

FUNCTION func(x: real): real;
BEGIN
   func := x*exp(-x)
END;

FUNCTION sngl(x : real): real;
BEGIN
   sngl := x
END;

(*$I GAULEG.PAS *)

BEGIN
   gauleg(x1,x2,x,w,npoint);
   writeln;
   writeln('#':2,'x[i]':11,'w[i]':12);
   FOR i := 1 to npoint DO BEGIN
      writeln(i:2,x[i]:12:6,w[i]:12:6)
   END;
(* demonstrate the use of gauleg for an integral *)
   gauleg(x1,x3,x,w,npoint);
   xx := 0.0;
   FOR i := 1 to npoint DO BEGIN
      xx := xx+sngl(w[i])*func(sngl(x[i]))
   END;
   writeln;
   writeln('Integral from GAULEG: ',xx:12:6);
   writeln('Actual value: ',((1.0+x1)*exp(-x1)-(1.0+x3)*exp(-x3)):12:6)
END.