implementation
//PROGRAM D3R10
//Driver for routine QUA3D
uses
    unit2;
  {$R *.DFM}
Function FUNC(X, Y, Z:real):real;
begin
    FUNC:= Sqr(X) + Sqr(Y) + Sqr(Z);
end;

Function Z1(X, Y:real):real;
begin
    Z1:= -Sqrt(Abs(Sqr(XMAX) - Sqr(X) - Sqr(Y)));
end;

Function Z2(X, Y:real):real;
begin
    Z2:= Sqrt(Abs(Sqr(XMAX ) - Sqr(X ) - Sqr(Y)));
end;

Function Y1(X:real):real;
begin
    Y1:= -Sqrt(Abs(Sqr(XMAX) - Sqr(X)));
end;

Function Y2(X:real):real;
begin
    Y2:= Sqrt(Abs(Sqr(XMAX) - Sqr(X)));
end;

procedure QGAUSX(A, B:real; var SS:real);
var
    X1,W:array[0..5] of real;
    XM,XR,DX:real;
    J:INTEGER;
begin
    X1[1]:=0.1488743389;  X1[2]:=0.4333953941; X1[3]:=0.6794095682;
    X1[4]:=0.8650633666;  X1[5]:=0.9739065285;
    W[1]:=0.2955242247;   W[2]:=0.2692667193;  W[3]:=0.2190863625;
    W[4]:=0.1494513491;   W[5]:=0.0666713443;
    XM:=0.5 * (B + A);
    XR:=0.5 * (B - A);
    SS:=0;
    For J:=1 To 5 do
    begin
        DX:=XR * X1[J];
        SS:=SS + W[J] * (H(XM + DX) + H(XM - DX));
    end;
    SS:=XR * SS;
end;

procedure QGAUSY(A, B:real;var SS:real);
var
    X1,W:array[0..5] of real;
    XM,XR,DX:real;
    J:INTEGER;
begin
    X1[1]:= 0.1488743389; X1[2]:=0.4333953941; X1[3]:=0.6794095682;
    X1[4]:=0.8650633666;  X1[5]:=0.9739065285;
    W[1]:=0.2955242247;   W[2]:=0.2692667193;  W[3]:=0.2190863625;
    W[4]:=0.1494513491;   W[5]:=0.0666713443;
    XM:=0.5 * (B + A);
    XR:=0.5 * (B - A);
    SS:=0;
    For J:=1 To 5 do
    begin
        DX:=XR * X1[J];
        SS:=SS + W[J] * (G(XM + DX) + G(XM - DX));
    end;
    SS:=XR * SS;
end;

procedure QGAUSZ(A, B:real;var SS:real);
var
    X1,W:array[0..5] of real;
    XM,XR,DX:real;
    J:INTEGER;
begin
    X1[1]:= 0.1488743389; X1[2]:=0.4333953941; X1[3]:=0.6794095682;
    X1[4]:=0.8650633666;  X1[5]:=0.9739065285;
    W[1]:=0.2955242247;   W[2]:=0.2692667193;  W[3]:=0.2190863625;
    W[4]:=0.1494513491;   W[5]:=0.0666713443;
    XM:=0.5 * (B + A);
    XR:=0.5 * (B - A);
    SS:=0;
    For J:=1 To 5 do
    begin
        DX:=XR * X1[J];
        SS:=SS + W[J] * (F(XM + DX) + F(XM - DX));
    end;
    SS:=XR * SS;
end;

procedure TForm1.Button1Click(Sender: TObject);
var
  F:TextFile;   I:Integer;
  XMIN,S:real;
const
  s1='%12.4f'; s2='0.#0';
  PI=3.1415926;  NVAL=10;
begin
  //輸出計算結果到文件
  AssignFile(F, 'd:\delphi_shu\p3\d3r10.dat');
  Rewrite(F);
  Writeln(F, 'Integral of r^2 over a spherical volume');
  Writeln(F, '   Radius     QUAD3D      Actual');
  For I:= 1 To NVAL do
  begin
      XMAX:=0.1 * I;
      XMIN:=-XMAX;
      QUAD3D(XMIN, XMAX, S);
      Writeln(F,'    ',FormatFloat(s2,XMAX),Format(s1,[S]),
           Format(s1,[4*PI*EXP(5*Ln(XMAX))/5]));
  end;
  CloseFile(F);
  //屏幕顯示計算結果
  memo1.Lines.LoadFromFile('d:\delphi_shu\p3\d3r10.dat');
end;