!*********************************************************** 
!* This program calculates the F Distribution Q(F|nu1,nu2) *
!* for F, nu1, nu2 given.                                  *
!* ------------------------------------------------------- *
!* Ref.: "Numerical Recipes  By W.H. Press, B.P. Flannery, *
!* S.A. Teukolsky and T. Vetterling, Cambridge University  *
!* Press, 1986, page 166" [BIBLI 08].                      *
!* ------------------------------------------------------- *
!* SAMPLE RUN:                                             *
!*                                                         *
!* Calculate F Distribution Function Q(F|nu1,nu2) for      * 
!* F, nu1 and nu2 given.                                   *
!*                                                         *
!* Ratio F of dispersion of 1st sample / dispersion of     *
!* 2nd sample: 1.25                                        *
!*                                                         *
!* Degree of freedom nu1 of first  sample: 25              *
!* Degree of freedom nu2 of second sample: 15              *
!*                                                         *
!*  Number of iterations:        6                         *
!*                                                         *
!*                                                         *
!*  Q =   0.3323728                                        *
!*                                                         *
!* ------------------------------------------------------- *
!*                       F90 version by J-P Moreau, Paris. *
!*                              (www.jpmoreau.fr)          *
!***********************************************************
PROGRAM FDISTRI
!**********************************************************************
!we use the formula: Q(F|nu1,nu2) = BetaI(nu2/2,nu1/2,nu2/(nu2+nu1*F))
!**********************************************************************
print *,''
print *,' Calculate F Distribution Function Q(F|nu1,nu2) for'
print *,' F, nu1 and nu2 given.'
print *,''
write(*,10,advance='no'); read *, f
print *,''
write(*,20,advance='no'); read *, xnu1
write(*,30,advance='no'); read *, xnu2
print *,''
tmp=xnu2/(xnu2+xnu1*f)

q=BETAI(xnu2/2.,xnu1/2.,tmp)

print *,''
print *,''
print *,' Q = ',q
print *,''
stop

10 format(' Ratio F of dispersion of 1st sample / dispersion of 2nd sample: ') 
20 format(' Degree of freedom of first  sample: ')
30 format(' Degree of freedom of second sample: ')

END


!*******************************************************************
!Returns the value ln(Gamma(xx)) for xx>0. Full accuracy is obtained
!for xx > 1. For 0<xx<1, the reflection formula can be used first:
!
!    Gamma(1-z) = pi/Gamma(z)/sin(pi*z) = pi*z/Gamma(1+z)/sin(pi*z)
!******************************************************************* 
Function GAMMLN(xx)
real*8 cof(6),stp,half,one,fpf,x,tmp,ser !internal arithmetic in double precision
data cof,stp/76.18009173D0,-86.50532033D0,24.01409822D0,         &
     -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/
data half,one,fpf/0.5d0,1.d0,5.5d0/
  x=xx-one
  tmp=x+fpf
  tmp=(x+half)*log(tmp)-tmp
  ser=one
  do j=1,6
    x=x+one
	ser=ser+cof(j)/x
  end do
  GAMMLN=tmp+log(stp*ser)
  return
End

!********************************************
!Returns the incomplete beta function Ix(a,b)
!********************************************
Function BETAI(a,b,x)
  if(x.lt.0..or.x.gt.1.) Pause ' Bad argument x in function BetaI.'
  if(x.eq.0..or.x.eq.1.) then
    bt=0.
  else
    tmp=gammln(a+b)-gammln(a)-gammln(b)+a*alog(x)+b*alog(1.-x)
    bt=exp(gammln(a+b)-gammln(a)-gammln(b) &
	   +a*alog(x)+b*alog(1.-x))
  end if
  if (x.lt.(a+1.)/(a+b+2.)) then !use continued fraction directly
    BETAI=bt*BETACF(a,b,x)/a
	return
  else
    BETAI=1.-bt*BETACF(b,a,1.-x)/b !use continued fraction after
	return						   !making the symmetry transformation
  end if
End

!******************************************************************
!  Continued fraction for incomplete beta function (used by BETAI)
!******************************************************************
Function BETACF(a,b,x)
Parameter(MAXIT=100,EPS=3e-7,FPMIN=1e-30)
  qab=a+b; qap=a+1.0; qam=a-1.0
  c=1.0
  d=1.0-qab*x/qap
  if (abs(d) < FPMIN) d=FPMIN
  d=1.0/d; h=d
  do m=1, MAXIT
    m2=2*m
	aa=m*(b-m)*x/((qam+m2)*(a+m2))
    d=1.0+aa*d
    if (abs(d) < FPMIN) d=FPMIN
    c=1.0+aa/c
	if (abs(c) < FPMIN) c=FPMIN
	d=1.0/d
	h=h*d*c
    aa=-(a+m)*(qab+m)*x/((qap+m2)*(a+m2))
    d=1.0+aa*d
    if (abs(d) < FPMIN) d=FPMIN
	c=1.0+aa/c
	if (abs(c) < FPMIN) c=FPMIN
	d=1.0/d
	del=d*c
	h=h*del
    if (abs(del-1.0) < EPS) goto 10
  end do !m loop
10 if (m>MAXIT) then 
	print *,' a or b too big, or MAXIT too small.'
  else 
    print *,' Number of iterations: ', m
  end if
  BETACF=h
  return
END

!end of file fdistri.f90