Time sequence and spectral analysis - Earth Dynamics: Deformations and Oscillations of the Rotating Earth - page 191

positive A and non-negative X, using up to 100 terms in the power series expansion.
The expansion is continued until the last term is negligible compared to the series
sum to double precision accuracy.

SUBROUTINE SEGAMI(SEGI,A,X)

C
C Subroutine to calculate the incomplete gamma function by power series.
C

IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C Check if X is negative or A is not positive.
IF(X.LT.0.D0.OR.A.LE.0.D0) GO TO 10
C Find natural logarithm of gamma(x).
ALG=ALOGAM(A)
C Set maximum number of terms.
ITM=100
C Test for zero argument.
IF(X.EQ.0.D0) GO TO 11
C Set value of first term of series sum.
SUM=1.D0/A
C Set value of first term.
TERM=SUM
C Set initial value of inverse of coefficient of X N .
COEF=A
C Sum series.
DO 12 N=1,ITM
C Update inverse of coefficient of X N .
COEF=COEF+1.D0
C Update term.
TERM=TERM*X/COEF
C Store current value of sum.
SUMOLD=SUM
C Update sum.
SUM=SUM+TERM
C Check relative size of current term.
IF(SUM.EQ.SUMOLD) GO TO 13

12 CONTINUE

PAUSE 'A too small or X too large.'

13 SEGI=SUM*DEXP(-X+A*DLOG(X)-ALG)

RETURN

10 PAUSE 'X is negative or A is not positive.'

RETURN

11 SEGI=0.D0

RETURN
END

a , the incomplete gamma function is found from (2.405)
using the continued fraction expansion (2.407). Evaluation of a continued fraction
of the form

For y greater than 1 +

a 1

= +

f b 0 (2.410)

a 2

b 1 +

a 3

b 2 +

a 4

b 3 +

a 5

b 4 +

b 5 +···