PROGRAM d4r7(input,output);
(* driver for routine QGAUS *)
CONST
   x1=0.0;
   x2=5.0;
   nval=10;
VAR
   dx,ss,x : real;
   i : integer;

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

(*$I MODFILE.PAS *)
(*$I QGAUS.PAS *)

BEGIN
   dx := (x2-x1)/nval;
   writeln;
   writeln('0.0 to','qgaus':10,'expected':13);
   writeln;
   FOR i := 1 to nval DO BEGIN
      x := x1+i*dx;
      qgaus(x1,x,ss);
      writeln(x:5:2,ss:12:6,(-(1.0+x)*exp(-x)
         +(1.0+x1)*exp(-x1)):12:6)
   END
END.

PROGRAM d4r4(input,output);
(* driver for routine QROMB *)
CONST
   pio2=1.5707963;
   n=5;
TYPE
   glnarray = ARRAY [1..n] OF real;
VAR
   glit : integer;
   a,b,s : real;

FUNCTION func(x: real): real;
(* Test function *)
BEGIN
   func := sqr(x)*(sqr(x)-2.0)*sin(x)
END;

FUNCTION fint(x: real): real;
(* Integral of test function func *)
BEGIN
   fint := 4.0*x*(sqr(x)-7.0)*sin(x)
      -(sqr(sqr(x))-14.0*sqr(x)+28.0)*cos(x)
END;

(*$I MODFILE.PAS *)
(*$I TRAPZD.PAS *)

(*$I POLINT.PAS *)

(*$I QROMB.PAS *)

BEGIN
   a := 0.0;
   b := pio2;
   writeln ('Integral of func computed with QROMB');
   writeln;
   writeln ('Actual value of integral is',fint(b)-fint(a):12:6);
   qromb(a,b,s);
   writeln ('Result from routine QROMB is',s:12:6);
END.

PROGRAM d4r6(input,output);
(* Driver for routine QROMO *)
CONST
   x1=0.0;
   x2=1.5707963;
   n=5;
TYPE
   glnarray = ARRAY [1..n] OF real;
VAR
   glit : integer;
   result : real;

FUNCTION func(x: real): real;
BEGIN
   func :=sqrt(x)/sin(x)
END;

(*$I MODFILE.PAS *)
(*$I MIDSQL.PAS *)

(*$I POLINT.PAS *)

(*$I QROMO.PAS *)

BEGIN
   writeln('Improper integral:');
   writeln;
   qromo(x1,x2,result);
   writeln('Function: sqrt(x)/sin(x)       Interval: (0,pi/2)');
   writeln('Using: MIDSQL                  Result:',result:8:4);
END.

PROGRAM d4r3(input,output);
(* driver for routine QSIMP *)
CONST
   pio2=1.5707963;
VAR
   glit : integer;
   a,b,s : real;

FUNCTION func(x : real) : real;
(* Test function *)
BEGIN
   func := sqr(x)*(sqr(x)-2.0)*sin(x)
END;

FUNCTION fint(x : real) : real;
(* Integral of test function *)
BEGIN
   fint := 4.0*x*(sqr(x)-7.0)*sin(x)-
         (sqr(sqr(x))-14.0*sqr(x)+28.0)*cos(x);
END;

(*$I MODFILE.PAS *)
(*$I TRAPZD.PAS *)

(*$I QSIMP.PAS *)

BEGIN
   a := 0.0;
   b := pio2;
   writeln ('Integral of func computed with QSIMP');
   writeln;
   writeln ('Actual value of integral is',fint(b)-fint(a):12:6);
   qsimp(a,b,s);
   writeln ('Result from routine QSIMP is',s:12:6)
END.

PROGRAM d4r2(input,output);
(* driver for routine QTRAP *)
CONST
   pio2=1.5707963;
VAR
   glit : integer;
   a,b,s : real;

FUNCTION func(x : real) : real;
(* Test function *)
BEGIN
   func := sqr(x)*(sqr(x)-2.0)*sin(x)
END;

FUNCTION fint(x : real) : real;
(* Integral of test function *)
BEGIN
   fint := 4.0*x*(sqr(x)-7.0)*sin(x)-
         (sqr(sqr(x))-14.0*sqr(x)+28.0)*cos(x);
END;

(*$I MODFILE.PAS *)
(*$I TRAPZD.PAS *)

(*$I QTRAP.PAS *)

BEGIN
   a := 0.0;
   b := pio2;
   writeln ('Integral of func computed with QTRAP');
   writeln;
   writeln ('Actual value of integral is',fint(b)-fint(a):12:6);
   qtrap(a,b,s);
   writeln ('Result from routine QTRAP is',s:12:6);
END.
PROGRAM d9r12(input,output);
(* driver for routine QROOT *)
CONST
   n=7;
   nv=3;
   eps=1.0e-6;
   ntry=10;
   tiny=1.0e-5;
TYPE
   glnarray = ARRAY [1..n] OF real;
   glnvarray = ARRAY [1..nv] OF real;
VAR
   i,j,nflag,nroot : integer;
   p : glnarray;
   b,c : ARRAY [1..ntry] OF real;

(*$I MODFILE.PAS *)
(*$I POLDIV.PAS *)

(*$I QROOT.PAS *)

BEGIN
   p[1] := 10.0; p[2] := -18.0; p[3] := 25.0; p[4] := -24.0;
   p[5] := 16.0; p[6] := -6.0; p[7] := 1.0;
   writeln;
   writeln('P(x) := x^6-6x^5+16x^4-24x^3+25x^2-18x+10');
   writeln('Quadratic factors x^2+bx+c');
   writeln;
   writeln('factor':6,'b':10,'c':12);
   writeln;
   nroot := 0;
   FOR i := 1 to ntry DO BEGIN
      c[i] := 0.5*i;
      b[i] := -0.5*i;
      qroot(p,n,b[i],c[i],eps);
      IF  (nroot = 0)  THEN BEGIN
         writeln(nroot:4,'   ',b[i]:12:6,c[i]:12:6);
         nroot := 1
      END ELSE BEGIN
         nflag := 0;
         FOR j := 1 to nroot DO BEGIN
            IF ((abs(b[i]-b[j]) < tiny)
               AND (abs(c[i]-c[j]) < tiny)) 
               THEN nflag := 1
         END;
         IF  (nflag = 0)  THEN BEGIN
            writeln(nroot:4,'   ',b[i]:12:6,c[i]:12:6);
            nroot := nroot+1
         END
      END
   END
END.
PROGRAM d4r9(input,output);
(* driver for routine QUAD3D *)
CONST
   pi=3.1415926;
   nval=10;
VAR
   i : integer;
   s,xmax,xmin,xmax5 : real;
   glx,gly : real;

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(sqr(xmax)-sqr(x)-sqr(y))
END;

FUNCTION z2(x,y: real): real;
BEGIN
   z2 := sqrt(sqr(xmax)-sqr(x)-sqr(y))
END;

FUNCTION y1(x: real): real;
BEGIN
   y1 := -sqrt(sqr(xmax)-sqr(x))
END;

FUNCTION y2(x: real): real;
BEGIN
   y2 := sqrt(sqr(xmax)-sqr(x))
END;

(*$I MODFILE.PAS *)

(*$I QUAD3D.PAS *)

BEGIN
   writeln('Integral of r^2 over a spherical volume');
   writeln;
   writeln('radius':13,'QUAD3D':9,'Actual':10);
   FOR i := 1 to nval DO BEGIN
      xmax := 0.1*i;
      xmin := -xmax;
      quad3d(xmin,xmax,s);
      xmax5 := sqr(sqr(xmax))*xmax;
      writeln(xmax:12:2,s:10:4,4.0*pi*(xmax5)/5.0:10:4)
   END
END.
