PROGRAM d6r29(input,output);
(* driver for routine CEL *)

(*$I MODFILE.PAS *)

CONST
   pio2 = 1.5707963;
   n = 5;
TYPE
   RealArray55 = ARRAY [1..55] OF real;
   RealArrayNP = ARRAY [1..n] OF real;
VAR
   TrapzdIt: integer;
   Ran3Inext,Ran3Inextp: integer;
   Ran3Ma: RealArray55;
   FuncA,FuncB,FuncP,FuncAkc: real;
   ago,astop,s: real;
   i,idum: integer;

FUNCTION func(phi: real): real;
(* Programs using routine FUNC must declare the variables
VAR
   FuncA,FuncB,FuncP,FuncAkc: real;
in the main routine. *)
VAR
   cs,csq,ssq: real;
BEGIN
   cs := cos(phi);
   csq := cs*cs;
   ssq := 1.0-csq;
   func := (FuncA*csq+FuncB*ssq)/(csq+FuncP*ssq)
      /sqrt(csq+FuncAkc*FuncAkc*ssq)
END;

(*$I RAN3.PAS *)

(*$I TRAPZD.PAS *)

(*$I POLINT.PAS *)

(*$I QROMB.PAS *)

(*$I CEL.PAS *)

BEGIN
   writeln('complete elliptic integral');
   writeln('kc':7,'p':10,'a':10,'b':10,'cel':11,'integral':12);
   idum := -55;
   ago := 0.0;
   astop := pio2;
   FOR i := 1 TO 20 DO BEGIN
      FuncAkc := 0.1+ran3(idum);
      FuncA := 10.0*ran3(idum);
      FuncB := 10.0*ran3(idum);
      FuncP := 0.1+ran3(idum);
      qromb(ago,astop,s);
      writeln(FuncAkc:10:6,FuncP:10:6,FuncA:10:6,FuncB:10:6,
         cel(FuncAkc,FuncP,FuncA,FuncB):10:6,s:10:6)
   END
END.