C @(#)distri.for	19.1 (ES0-DMD) 02/25/03 13:21:22
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

        REAL FUNCTION GAMMAQ(A,X)
        REAL X,A
        REAL GLN,GAMSER,GAMMCF
        IF (X.LT.A+1)THEN
           CALL GSER(GAMSER,A,X,GLN)
           GAMMAQ=1-GAMSER
        ELSE
           CALL GCF(GAMMCF,A,X,GLN)
           GAMMAQ=GAMMCF
        ENDIF
        RETURN
        END

        REAL FUNCTION GAMMAP(A,X)
        REAL X,A
        REAL GLN,GAMSER,GAMMCF
        IF (X.LT.A+.1)THEN
           CALL GSER(GAMSER,A,X,GLN)
           GAMMAP=GAMSER
        ELSE
           CALL GCF(GAMMCF,A,X,GLN)
           GAMMAP=1.-GAMMCF
        ENDIF
        RETURN
        END



               
        SUBROUTINE GSER(GAMSER,A,X,GLN)
        REAL GAMSER,X,GLN,SUM,DEL,EPS,A
        REAL GAMMLN,AP
        INTEGER ITMAX,N
        PARAMETER (ITMAX=100,EPS=3.E-7)
        GLN=GAMMLN(A)
        IF (X.LE.0)THEN
           IF(X.LT.0)RETURN
           GAMSER=0
           RETURN
        ENDIF
        AP = A
       SUM = 1./A
        DEL = SUM
        DO 10 N=1,ITMAX
           AP = AP+1
           DEL = DEL*X/AP
           SUM = SUM+DEL
           IF (ABS(DEL).LT.ABS(SUM)*EPS)GOTO 20
10      CONTINUE
20       GAMSER=SUM*EXP(-X+A*LOG(X)-GLN)
        RETURN
       END

        SUBROUTINE GCF(GAMMCF,A,X,GLN)
        REAL GAMMCF,X,GLN,GOLD,A0,A1,B0,B1,FAC,EPS,ANA,G,A
        REAL GAMMLN,ANF
        INTEGER N,ITMAX
        PARAMETER (ITMAX=100,EPS=3.E-7)
        GLN=GAMMLN(A)
        A0 = 1.
        A1 = X
        B0 = 0
        B1 = 1
        FAC = 1
        DO 10 N = 1,ITMAX
           ANA = N-A
           A0 = (A1+A0*ANA)*FAC
           B0 = (B1+B0*ANA)*FAC
           ANF = N*FAC
           A1 = X*A0+ANF*A1
           B1 = X*B0+ANF*B1
           IF(A1.NE.0)THEN
             FAC = 1./A1
             G=B1*FAC
             IF(ABS((G-GOLD)/G).LT.EPS)GOTO 20
             GOLD=G
           ENDIF
10      CONTINUE
20      GAMMCF = EXP(-X+A*ALOG(X)-GLN)*G
        RETURN
        END        

  
        REAL FUNCTION GAMMLN(XX)
        REAL XX
        INTEGER J
        DOUBLE PRECISION COF(6),STP,HALF,ONE,FPF,X,TMP,SER
        DATA COF/76.18009173,-86.50532033,24.01409822,
     1        -1.231739516,0.120858003D-2,-0.536382D-5/
        DATA STP/2.50662827465/
        DATA HALF,ONE,FPF/0.5,1.,5.5/
        X = XX - ONE
        TMP = X + FPF
        TMP = (X+HALF)*LOG(TMP)-TMP
        SER = ONE
        DO 10 J=1,6
           X = X + ONE
           SER = SER + COF(J)/X
10      CONTINUE
        GAMMLN = TMP + LOG(STP*SER)
        RETURN
        END


       REAL FUNCTION ERF(X)
       REAL X,GAMMAP
c       returns the error function erf(X)
       IF(X.LT.0) THEN
           ERF = -1*GAMMAP(0.5,X**2)
       ELSE
           ERF = GAMMAP(0.5,X**2)
       ENDIF
       RETURN
       END

       REAL FUNCTION BETAI(A,B,X)
C       returns the imcomplete beta function
       REAL*4 A,B,X,GAMMLN,BT,BETACF       
       IF (X.EQ.0.OR.X.EQ.1.)THEN
           BT = 0.
       ELSE
           BT = EXP(GAMMLN(A+B)-GAMMLN(A)-GAMMLN(B)+
     1            A*LOG(X)+B*LOG(1.-X))
       ENDIF
       IF (X.LT.(A+1.)/(A+B+2.))THEN      
           BETAI = BT*BETACF(A,B,X)/A
       ELSE
           BETAI = 1.-1*(BT*BETACF(B,A,1.-X)/B)
       ENDIF
       END
       
       REAL FUNCTION BETACF(A,B,X)
       REAL A,B,X,TEM,APP,BPP,AOLD
       REAL EPS,AM,BM,AZ,QAB,QAM,QAP,BZ,EM,D,AP,BP
       INTEGER ITMAX,M
       PARAMETER (ITMAX=100,EPS=3.E-7)
       AM = 1.
       BM = 1.
       AZ = 1.
       QAB = A+B
       QAP = A+1.
       QAM = A-1.
       BZ = 1.-QAB*X/QAP
       DO 10 M = 1,ITMAX
           EM = M
           TEM = EM+EM
           D = EM*(B-M)*X/((QAM+TEM)*(A+TEM))
           AP = AZ+D*AM
           BP = BZ+D*BM
           D = -(A+EM)*(QAB+EM)*X/((A+TEM)*(QAP+TEM))
           APP = AP+D*AZ
           BPP = BP+D*BZ
           AOLD = AZ
           AM = AP/BPP
           BM = BP/BPP
           AZ = APP/BPP
           BZ = 1.
           IF(ABS(AZ-AOLD).LT.EPS*ABS(AZ))GOTO 20
10        CONTINUE
20       BETACF = AZ
        RETURN
        END                    

       REAL FUNCTION ERFCC(X)
       REAL X,Z,T
       Z = ABS(X)
       T = 1./(1.+0.5*Z)
       ERFCC = T*EXP(-Z*Z-1.26551223+T*(1.00002368+T*(0.37409196+
     1        T*(0.09678418+T*(-0.18628806+T*(0.27886807+T*(-1.13520398+
     2     T*(1.48851587+T*(-0.82215223+T*0.17087277)))))))))
       IF (X.LT.0.) ERFCC = 2.-ERFCC
       RETURN
       END


       REAL FUNCTION UNIF(X,PARA)
       REAL X,PARA(2),A,B
       A = PARA(1)
       B = PARA(2)
       IF (X.LT.B.AND.X.GT.A)THEN
         UNIF = (B-X)/(B-A)
       ELSE
         UNIF = 0
       ENDIF
       END


       REAL FUNCTION BINOM(NOBS,N,PROB)
       REAL FACT
       INTEGER NOBS,N,I
       REAL PROB
       BINOM = 0.
       DO 10 I = 0,NOBS
         BINOM = BINOM+FACT(N)/FACT(I)/FACT(N-I)*(PROB**I)*
     1              ((1-PROB)**(N-I))
10     CONTINUE
       RETURN
       END

       REAL FUNCTION FACT(N)
       REAL GAMMLN
       INTEGER N,J
       FACT = 1.
       IF (N.LT.0) THEN
          RETURN
       ELSE IF (N.LT.32) THEN
          DO 11 J=1,N
              FACT = FACT*J
11        CONTINUE          
       ELSE
         FACT = EXP(GAMMLN(N+1.))
       ENDIF
       RETURN
       END

       
       
