(*-------------------------------------------------------------------------*)
(*                GammaIn -- Incomplete Gamma integral                     *)
(*-------------------------------------------------------------------------*)

FUNCTION GammaIn(     Y, P    : REAL;
                      Dprec   : INTEGER;
                      MaxIter : INTEGER;
                  VAR Cprec   : REAL;
                  VAR Iter    : INTEGER;
                  VAR Ifault  : INTEGER ) : REAL;

(*-------------------------------------------------------------------------*)
(*                                                                         *)
(*       Function:  GammaIn                                                *)
(*                                                                         *)
(*       Purpose:   Evaluates Incomplete Gamma integral                    *)
(*                                                                         *)
(*       Calling Sequence:                                                 *)
(*                                                                         *)
(*            P     := GammaIn(       Y, P    : REAL;                      *)
(*                                    Dprec   : INTEGER;                   *)
(*                                    MaxIter : INTEGER;                   *)
(*                               VAR  Cprec   : REAL;                      *)
(*                               VAR  Iter    : INTEGER;                   *)
(*                               VAR  Ifault  : INTEGER ) : REAL;          *)
(*                                                                         *)
(*                 Y      --- Gamma distrib. value                         *)
(*                 P      --- Degrees of freedom                           *)
(*                 Ifault --- error indicator                              *)
(*                                                                         *)
(*                 P      --- Resultant probability                        *)
(*                                                                         *)
(*       Calls:                                                            *)
(*                                                                         *)
(*            None                                                         *)
(*                                                                         *)
(*       Remarks:                                                          *)
(*                                                                         *)
(*            Either an infinite series summation or a continued fraction  *)
(*            approximation is used, depending upon the argument range.    *)
(*                                                                         *)
(*-------------------------------------------------------------------------*)

CONST
   Oflo    = 1.0E+37;
   MinExp  = -87.0;

VAR
   F:        REAL;
   C:        REAL;
   A:        REAL;
   B:        REAL;
   Term:     REAL;
   Pn:       ARRAY[1..6] OF REAL;
   Gin:      REAL;
   An:       REAL;
   Rn:       REAL;
   Dif:      REAL;
   Eps:      REAL;
   Done:     BOOLEAN;

LABEL 9000;

BEGIN (* GammaIn *)
                                   (* Check arguments *)
   Ifault  := 1;
   GammaIn := 1.0;

   IF( ( Y <= 0.0 ) OR ( P <= 0.0 ) ) THEN GOTO 9000;

                                   (* Check value of F *)
   Ifault := 0;
   F      := P * LN( Y ) - ALGama( P + 1.0 ) - Y;
   IF ( F < MinExp ) THEN GOTO 9000;

   F := EXP( F );
   IF( F = 0.0 ) THEN GOTO 9000;

                                  (* Set precision *)
   IF Dprec > MaxPrec THEN
      Dprec := MaxPrec
   ELSE IF Dprec <= 0 THEN
      Dprec := 1;

   Cprec  := Dprec;

   Eps    := PowTen( -Dprec );

                                   (* Choose infinite series or *)
                                   (* continued fraction.       *)

   IF( ( Y > 1.0 ) AND ( Y >= P ) ) THEN

      BEGIN (* Continued Fraction *)

         A       := 1.0 - P;
         B       := A + Y + 1.0;
         Term    := 0.0;
         Pn[ 1 ] := 1.0;
         Pn[ 2 ] := Y;
         Pn[ 3 ] := Y + 1.0;
         Pn[ 4 ] := Y * B;
         Gin     := Pn[ 3 ] / Pn[ 4 ];
         Done    := FALSE;
         Iter    := 0;

         REPEAT

            Iter   := Iter + 1;
            A      := A + 1.0;
            B      := B + 2.0;
            Term   := Term + 1.0;
            An     := A * Term;

            Pn[ 5 ] := B * Pn[ 3 ] - An * Pn[ 1 ];
            Pn[ 6 ] := B * Pn[ 4 ] - An * Pn[ 2 ];

            IF( Pn[ 6 ] <> 0.0 ) THEN
               BEGIN

                  Rn     := Pn[ 5 ] / Pn[ 6 ];
                  Dif    := ABS( Gin - Rn );

                  IF( Dif <= Eps ) THEN
                     IF( Dif <= ( Eps * Rn ) ) THEN Done := TRUE;

                  Gin    := Rn;

               END;

            Pn[ 1 ] := Pn[ 3 ];
            Pn[ 2 ] := Pn[ 4 ];
            Pn[ 3 ] := Pn[ 5 ];
            Pn[ 4 ] := Pn[ 6 ];

            IF( ABS( Pn[ 5 ] ) >= Oflo ) THEN
               BEGIN
                  Pn[ 1 ] := Pn[ 1 ] / Oflo;
                  Pn[ 2 ] := Pn[ 2 ] / Oflo;
                  Pn[ 3 ] := Pn[ 3 ] / Oflo;
                  Pn[ 4 ] := Pn[ 4 ] / Oflo;
               END;

         UNTIL ( Iter > MaxIter ) OR Done;

         Gin    := 1.0 - ( F * Gin * P );

         GammaIn := Gin;
                                   (* Calculate precision of result *)

         IF Dif <> 0.0 THEN
            Cprec := -LogTen( Dif )
         ELSE
            Cprec := MaxPrec;

      END   (* Continued Fraction *)

   ELSE

      BEGIN (* Infinite series *)

         Iter    := 0;

         Term    := 1.0;
         C       := 1.0;
         A       := P;
         Done    := FALSE;

         REPEAT

            A       := A + 1.0;
            Term    := Term * Y / A;
            C       := C + Term;

            Iter    := Iter + 1;

         UNTIL ( Iter > MaxIter ) OR ( ( Term / C ) <= Eps );

         GammaIn := C * F;
                                   (* Calculate precision of result *)

         Cprec := -LogTen( Term / C );

      END   (* Infinite series *);

9000: (* Error exit *)

END    (* GammaIn *);