PROGRAM d6r28(input,output);
(* driver for routine EL2 *)
VAR
   glit : integer;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gla,glb,glakc : real;
   ago,astop,s,x : real;
   i,idum : integer;

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

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

(*$I TRAPZD.PAS *)

(*$I QSIMP.PAS *)

(*$I EL2.PAS *)

BEGIN
   writeln('general elliptic integral of second kind');
   writeln('x':7,'kc':10,'a':10,'b':10,'el2':11,'integral':12);
   idum := -55;
   ago := 0.0;
   FOR i := 1 to 20 DO BEGIN
      glakc := 5.0*ran3(idum);
      gla := 10.0*ran3(idum);
      glb := 10.0*ran3(idum);
      x := 10.0*ran3(idum);
      astop := arctan(x);
      qsimp(ago,astop,s);
      writeln(x:10:6,glakc:10:6,gla:10:6,glb:10:6,
         el2(x,glakc,gla,glb):10:6,s:10:6)
   END
END.