PROGRAM d6r29(input,output);
(* driver for routine CEL *)
CONST
   pio2=1.5707963;
VAR
   glit : integer;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gla,glb,glp,glakc : real;
   ago,astop,s : real;
   i,idum : integer;

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

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

(*$I TRAPZD.PAS *)

(*$I QSIMP.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
      glakc := 0.1+ran3(idum);
      gla := 10.0*ran3(idum);
      glb := 10.0*ran3(idum);
      glp := 0.1+ran3(idum);
      qsimp(ago,astop,s);
      writeln(glakc:10:6,glp:10:6,gla:10:6,glb:10:6,
         cel(glakc,glp,gla,glb):10:6,s:10:6)
   END
END.