bief_gcf.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\bief_gcf.f
00002 !
00045                          SUBROUTINE BIEF_GCF
00046 !                        *******************
00047 !
00048      &(GAMMCF,A,X,GLN)
00049 !
00050 !***********************************************************************
00051 ! BIEF   V6P3                                   07/03/2013
00052 !***********************************************************************
00053 !
00054 !
00055 !
00056 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00057 !| A              |-->| PARAMETER
00058 !| GAMMCF         |<->| Incomplete gamma function Q(a,x)
00059 !| GLN            |<->| gamma(a)
00060 !| X              |-->| OPERAND
00061 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00062 !
00063 !     USE BIEF
00064 !
00065       IMPLICIT NONE
00066       INTEGER LNG,LU
00067       COMMON/INFO/LNG,LU
00068 !
00069 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00070 !
00071       DOUBLE PRECISION, INTENT(IN)    :: A,X
00072       DOUBLE PRECISION, INTENT(INOUT) :: GAMMCF,GLN
00073 !
00074 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00075 !
00076       INTEGER ITMAX,N
00077       DOUBLE PRECISION EPS
00078       DOUBLE PRECISION GOLD,A0,A1,B0,B1,FAC,AN,ANA,ANF,G
00079 !
00080       PARAMETER(ITMAX=100,EPS=3.E-7)
00081 !
00082       DOUBLE PRECISION BIEF_GAMMLN
00083       EXTERNAL         BIEF_GAMMLN
00084 !
00085       INTRINSIC EXP,LOG,FLOAT,ABS
00086 !
00087 !-----------------------------------------------------------------------
00088 !
00089       GLN = BIEF_GAMMLN(A)
00090 !     This is the previous value, tested against for convergence.
00091       GOLD = 0.D0
00092       A0   = 1.D0
00093 !     We are here setting up the A's and B's of equation (5.2.4) for
00094 !     evaluating the continued fraction.
00095       A1     = X
00096       B0     = 0.D0
00097       B1     = 1.D0
00098 !     FAC is the renormalization factor for preventing overflow of the
00099 !     partial numerators and denominators.
00100       FAC    = 1.D0
00101       DO N = 1,ITMAX
00102         AN     = FLOAT(N)
00103         ANA    = AN - A
00104 !       One step of the recurrence (5.2.5).
00105         A0     = (A1+A0*ANA)*FAC
00106         B0     = (B1+B0*ANA)*FAC
00107         ANF    = AN*FAC
00108 !       The next step of the recurrence (5.2.5).
00109         A1     = X*A0 + ANF*A1
00110         B1     = X*B0 + ANF*B1
00111 !       Shall we renormalize?
00112         IF(A1.NE.0.D0) THEN
00113 !         Yes.  Set FAC so it happens.
00114           FAC    = 1.D0/A1
00115 !         New value of answer.
00116           G      = B1*FAC
00117 !         Converged? If so, exit.
00118           IF (ABS((G-GOLD)/G).LT.EPS) GO TO 110
00119 !         If not, save value.
00120           GOLD   = G
00121         ENDIF
00122       ENDDO
00123       IF(LNG.EQ.1) THEN
00124         WRITE(LU,*) 'BIEF_GCF : ERREUR'
00125       ENDIF
00126       IF(LNG.EQ.2) THEN
00127         WRITE(LU,*) 'BIEF_GCF: FATAL ERROR'
00128       ENDIF
00129       CALL PLANTE(1)
00130       STOP
00131 110   CONTINUE
00132 !
00133       GAMMCF = EXP(-X+A*LOG(X)-GLN)*G
00134 !
00135 !-----------------------------------------------------------------------
00136 !
00137       RETURN
00138       END

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