!******************************************************************
!*      Purpose: This program computes the beta function          * 
!*               B(p,q) for p > 0 and q > 0 using                 *
!*               subroutine BETA                                  *
!*      Input :  p  --- Parameter  ( p > 0 )                      *
!*               q  --- Parameter  ( q > 0 )                      *
!*      Output:  BT --- B(p,q)                                    *
!*      Examples:                                                 *
!*                p       q           B(p,q)                      *
!*              ---------------------------------                 *
!*               1.5     2.0     .2666666667D+00                  *
!*               2.5     2.0     .1142857143D+00                  *
!*               1.5     3.0     .1523809524D+00                  *
!* -------------------------------------------------------------- *
!* REFERENCE: "Fortran Routines for Computation of Special        *
!*             Functions jin.ece.uiuc.edu/routines/routines.html" *
!*                                                                *
!*                              F90 Release By J-P Moreau, Paris. *
!*                                      (www.jpmoreau.fr)         *
!******************************************************************
        PROGRAM MBETA
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        WRITE(*,*)
        WRITE(*,*)'    p       q           B(p,q)'
        WRITE(*,*)'  ---------------------------------'
	  P=1.5D0
	  Q=2.0D0
        CALL BETA(P,Q,BT)
        WRITE(*,10)P,Q,BT
	  P=2.5D0
	  Q=2.0D0
        CALL BETA(P,Q,BT)
        WRITE(*,10)P,Q,BT
	  P=1.5D0
	  Q=3.0D0
        CALL BETA(P,Q,BT)
        WRITE(*,10)P,Q,BT
	  WRITE(*,*)
10      FORMAT(2X,F5.1,3X,F5.1,D20.10)
        END


        SUBROUTINE BETA(P,Q,BT)

!       ==========================================
!       Purpose: Compute the beta function B(p,q)
!       Input :  p  --- Parameter  ( p > 0 )
!                q  --- Parameter  ( q > 0 )
!       Output:  BT --- B(p,q)
!       Routine called: GAMMA for computing â(x)
!       ==========================================

        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        CALL GAMMA(P,GP)
        CALL GAMMA(Q,GQ)
        PPQ=P+Q
        CALL GAMMA(PPQ,GPQ)
        BT=GP*GQ/GPQ
        RETURN
        END


        SUBROUTINE GAMMA(X,GA)

!       ==================================================
!       Purpose: Compute gamma function â(x)
!       Input :  x  --- Argument of â(x)
!                       ( x is not equal to 0,-1,-2,úúú)
!       Output:  GA --- â(x)
!       ==================================================

        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION G(26)
        PI=3.141592653589793D0
        IF (X.EQ.INT(X)) THEN
           IF (X.GT.0.0D0) THEN
              GA=1.0D0
              M1=X-1
              DO 10 K=2,M1
10               GA=GA*K
           ELSE
              GA=1.0D+300
           ENDIF
        ELSE
           IF (DABS(X).GT.1.0D0) THEN
              Z=DABS(X)
              M=INT(Z)
              R=1.0D0
              DO 15 K=1,M
15               R=R*(Z-K)
              Z=Z-M
           ELSE
              Z=X
           ENDIF
           DATA G/1.0D0,0.5772156649015329D0,  &
                -0.6558780715202538D0, -0.420026350340952D-1, &
                0.1665386113822915D0,-.421977345555443D-1, &
                -.96219715278770D-2, .72189432466630D-2, &
                -.11651675918591D-2, -.2152416741149D-3, &
                .1280502823882D-3, -.201348547807D-4, &
                -.12504934821D-5, .11330272320D-5, &
                -.2056338417D-6, .61160950D-8, &
                .50020075D-8, -.11812746D-8, &
                .1043427D-9, .77823D-11, &
                -.36968D-11, .51D-12, &
                -.206D-13, -.54D-14, .14D-14, .1D-15/
           GR=G(26)
           DO 20 K=25,1,-1
20            GR=GR*Z+G(K)
           GA=1.0D0/(GR*Z)
           IF (DABS(X).GT.1.0D0) THEN
              GA=GA*R
              IF (X.LT.0.0D0) GA=-PI/(X*GA*DSIN(PI*X))
           ENDIF
        ENDIF
        RETURN
        END

!end of file mbeta.f90