C @(#)dgamma.for	17.1.1.1 (ES0-DMD) 01/25/02 17:10:41
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
       DOUBLE PRECISION FUNCTION DGAMMA(W)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C M. Rosa ST-ECF ESO Garching 11-12-87
C      IMPLICIT NONE
C
C Double Precision Gamma-Function (15 digits) using modified version
C of Comm. ACM algorithm No 309 "Gamma function with arbitrary precision"
C by Souza and Schwachheim, CAMC Vol. 10, p. 511, 1967.
C
C Def.   Gamma(x+1) = x!  fuer alle x, X.NEQ.-1 ;  0! = 1.0
C
C  Note: 1. The Coefficients C(i) are given as irreducible fractions so that
C          accuracy is automatically machine appropriate.
C           (The original article uses 20 coefficients to obtain an accuracy
C          of up to  50 digits.
C
C       2. The routine evaluates LOG(Gamma(w)). It might be preferred to
C          pass LOGGAMMA instead of exp(LOGGAMMA) to avoid additional
C          inaccuracies. (For example if 1/DGAMMA is needed, EXP(-LOGGAMMA)
C          would do at once.
C
C ---------------------------------------------------------------------
C
      INTEGER I
      REAL W2
      DOUBLE PRECISION W
      DOUBLE PRECISION F,DEN,S,C0,C(10),PS
      DATA C0/0.9189385332046727417803D0/
C
      C(1)   = 1.D0/12.D0
      C(2)   = -1.D0/360.D0
      C(3)   = 1.D0/1260.D0
      C(4)   = -1.D0/1680.D0
      C(5)   = 1.D0/1188.D0
      C(6)   = -691.D0/360360.D0
      C(7)   = 1.D0/156.D0
      C(8)   = -3617.D0/122400.D0
      C(9)   = 43867.D0/244188.D0
      C(10)  = -174611.D0/125400.D0
C
      IF (W.LE.7.0D0) THEN
          F      = W
   10     W      = W + 1.0D0
          IF (W.LT.7.0D0) THEN
              F      = F*W
              GO TO 10

          END IF

          F      = -DLOG(F)

      ELSE
          F      = 0.0D0
      END IF
C
      DEN    = W
      W2     = W*W
      PS     = (W-0.5D0)*DLOG(W) - W + C0
      DO 20 I = 1,10
          S      = PS + C(I)/DEN
          IF (S.EQ.PS) GO TO 30
          DEN    = DEN*W2
          PS     = S
   20 CONTINUE
   30 S      = S + F
C
      DGAMMA = DEXP(S)
      END