bief_gser.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\bief_gser.f
00002 !
00047                          SUBROUTINE BIEF_GSER
00048 !                        ********************
00049 !
00050      &(GAMSER,A,X,GLN)
00051 !
00052 !***********************************************************************
00053 ! BIEF   V6P3                                   07/03/2013
00054 !***********************************************************************
00055 !
00056 !
00057 !
00058 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00059 !| A              |-->| PARAMETER
00060 !| GAMSER         |<->| Incomplete gamma function P(a,x)
00061 !| GLN            |<->| gamma(a)
00062 !| X              |-->| OPERAND
00063 !| IFAIL          |<->| ERROR TAG: 1 IS OK, 0 IS ERROR
00064 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00065 !
00066 !     USE BIEF
00067 !
00068       IMPLICIT NONE
00069       INTEGER LNG,LU
00070       COMMON/INFO/LNG,LU
00071 !
00072 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00073 !
00074       DOUBLE PRECISION, INTENT(IN)    :: X,A
00075       DOUBLE PRECISION, INTENT(INOUT) :: GAMSER,GLN
00076 !
00077 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00078 !
00079       INTEGER ITMAX,N
00080       DOUBLE PRECISION AP,SUM,DEL,EPS
00081 !
00082       PARAMETER (ITMAX=100,EPS=3.D-7)
00083 !
00084       DOUBLE PRECISION BIEF_GAMMLN
00085       EXTERNAL         BIEF_GAMMLN
00086 !
00087       INTRINSIC ABS,EXP,LOG
00088 !
00089 !-----------------------------------------------------------------------
00090 !
00091       GLN = BIEF_GAMMLN(A)
00092 !
00093       IF(X.LT.0.D0) THEN
00094 !       Cannot occur when used by ERF or ERFC
00095         IF(LNG.EQ.1) WRITE(LU,*) 'BIEF_GSER : X NEGATIF NON PREVU'
00096         IF(LNG.EQ.2) WRITE(LU,*) 'BIEF_GSER: UNEXPECTED NEGATIVE X'
00097         CALL PLANTE(1)
00098         STOP
00099       ELSEIF(X.EQ.0.D0) THEN
00100         GAMSER = 0.D0
00101       ELSE
00102         AP  = A
00103         SUM = 1.D0/A
00104         DEL = SUM
00105         DO N = 1,ITMAX
00106           AP  = AP + 1.D0
00107           DEL = DEL*X/AP
00108           SUM = SUM + DEL
00109           IF(ABS(DEL).LT.ABS(SUM)*EPS) GO TO 110
00110         ENDDO
00111         IF(LNG.EQ.1) THEN
00112           WRITE(LU,*) 'BIEF_GSER : ITERATION MAXIMUM DEPASSEE : ',ITMAX
00113         ENDIF
00114         IF(LNG.EQ.2) THEN
00115           WRITE(LU,*) 'BIEF_GSER: MAXIMUM OF ITERATIONS REACHED:',ITMAX
00116         ENDIF
00117         CALL PLANTE(1)
00118         STOP
00119 110     CONTINUE
00120         GAMSER = SUM*EXP(-X+A*LOG(X)-GLN)
00121       ENDIF
00122 !
00123 !-----------------------------------------------------------------------
00124 !
00125       RETURN
00126       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0