C++*********************************************************************
C
C GAMMLN.F
C
C **********************************************************************
C *                                                                        *
C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM.   AUTHOR: J.FRANK  *
C=* Copyright (C) 1985-2005  Health Research Inc.                      *
C=*                                                                    *
C=* HEALTH RESEARCH INCORPORATED (HRI),                                *   
C=* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455.                   *
C=*                                                                    *
C=* Email:  spider@wadsworth.org                                       *
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 the *
C=* 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 GNU  *
C=* General Public License for more details.                           *
C=*                                                                    *
C=* You should have received a copy of the GNU General Public License  *
C=* along with this program; if not, write to the                      *
C=* Free Software Foundation, Inc.,                                    *
C=* 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.      *
C=*                                                                    *
C **********************************************************************
C
C	RETURNS THE VALUE ln(XX) FOR XX > 0. FULL ACCURACY IS 
C	OBTAINED FOR  0 < XX < 1, THE REFLECTION FORMULA CAN BE USED
C       FIRST.
C
C IMAGE_PROCESSING_ROUTINE
C
C--*********************************************************************

        DOUBLE PRECISION FUNCTION GAMMLN(XX)

        DOUBLE PRECISION COF(6), STP, HALF, ONE, FPF, TMP, SER
        DOUBLE PRECISION X, XX
        INTEGER J


        DATA COF, STP/76.18009173D0, -86.50532033D0, 24.01409822D0,
     &          -1.231739516D0, .120858003D-2, -.536382D-5, 
     &          2.50662827465D0/
        DATA HALF, ONE, FPF/0.5D0, 1.0D0, 5.5D0/

        X = XX - ONE
        TMP = X + FPF
        TMP = (X + HALF) * DLOG(TMP) - TMP
        SER = ONE
        DO J = 1, 6
          X = X + ONE
          SER = SER + COF(J) / X
        END DO
        GAMMLN = TMP + DLOG(STP * SER)
        END