real function student_t ( Samples, Sign_Level ) 
! The student-t is calculated for a given number of samples
! and at a certain significance level.
! The degrees of freedom are calculated in the routine as
! Df = real( Samples - 1 ).
! The following relation holds:
! Sign_Level = 1 - IncompleteBetaFunction( x, a, b ) 
!                             x = Df / ( Df + student_t^2 )
!                             a = Df / 2
!                             b = 0.50 
! We need to sove the above equation for x (or student_t).
! Routines from Numerical Recipes are used for that.

implicit none

! Input variables.
integer, intent(in)		:: Samples     ! Number of samples
real, intent(in)		:: Sign_Level  ! Significance level

! Functions and parameters.
real, external			::  zbrent, tobesolved       ! Functions
real					::  x1=0.0, x2=1.0, b=0.50, eps=1.0e-7
!  Various parameters: x1, x2 bracket the root, given with 
!  accuracy eps.
           
! Local
real					::  Df  ! Degrees of freedom
real					::  a   ! = 0.50 * real( Samples - 1 )


Df = real( Samples - 1 )
a = 0.50 * Df
student_t = zbrent(tobesolved,a,b,Sign_Level, x1, x2, eps )
student_t = sqrt( Df/student_t  - Df )

end function student_t



real function tobesolved( a, b, c, x ) 
real, intent(in)	:: a, b, c, x
! a, b, c: parameters to the function
! x:       variable
real, external		:: betai
! Incomplete beta function.
tobesolved = betai(a,b,x) - 1.0 + c
end function tobesolved




! The rest of this file comes from Numerical Recipes.
! Function zbrent has been modified slightly 
! (variables aaa, bbb, ccc have been intoduced).

  	  ! Brent's method for solving the equation
	  ! func(a,b,c,x)=0 for x, where a,b,c parameters.
	  ! Root is bracketed by x1 and x2.
	  ! Root is returned to varable zbrent with
	  ! accuracy tol. 
	  real function zbrent(func,aaa,bbb,ccc,x1,x2,tol)
      integer, parameter	:: ITMAX = 500
      real, intent(in)		:: tol,x1,x2
	  real, external		:: func
      real, parameter		:: EPS=3.0e-8
      integer				:: iter
      real					:: a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm
	  real, intent(in)		:: aaa, bbb, ccc
      a=x1
      b=x2
      fa=func(aaa,bbb,ccc,a)
      fb=func(aaa,bbb,ccc,b)
      if((fa.gt.0.0.and.fb.gt.0.0).or.(fa.lt.0.0.and.fb.lt.0.0))  &
	     write(*,*) 'root must be bracketed for zbrent'
      c=b
      fc=fb
      do 11 iter=1,ITMAX
        if((fb.gt.0.0.and.fc.gt.0.0).or.(fb.lt.0.0.and.fc.lt.0.0))then
          c=a
          fc=fa
          d=b-a
          e=d
        endif
        if(abs(fc).lt.abs(fb)) then
          a=b
          b=c
          c=a
          fa=fb
          fb=fc
          fc=fa
        endif
        tol1=2.0*EPS*abs(b)+0.5*tol
        xm=0.5*(c-b)
        if(abs(xm).le.tol1 .or. fb.eq.0.0)then
          zbrent=b
          return
        endif
        if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then
          s=fb/fa
          if(a.eq.c) then
            p=2.0*xm*s
            q=1.0-s
          else
            q=fa/fc
            r=fb/fc
            p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0))
            q=(q-1.0)*(r-1.0)*(s-1.0)
          endif
          if(p.gt.0.0) q=-q
          p=abs(p)
          if(2.0*p .lt. min(3.0*xm*q-abs(tol1*q),abs(e*q))) then
            e=d
            d=p/q
          else
            d=xm
            e=d
          endif
        else
          d=xm
          e=d
        endif
        a=b
        fa=fb
        if(abs(d) .gt. tol1) then
          b=b+d
        else
          b=b+sign(tol1,xm)
        endif
        fb=func(aaa,bbb,ccc,b)
11    continue
      write(*,*) 'zbrent exceeding maximum iterations'
      zbrent=b
      return
      
	  end function zbrent




      ! Incomplete beta function.
      real function betai(a,b,x)
      real, intent(in)	:: a,b,x
!U    USES betacf,gammln
      real				:: bt
	  real, external	:: betacf,gammln
      if(x.lt.0.0.or.x.gt.1.0)write(*,*) 'bad argument x in betai'
      if(x.eq.0.0.or.x.eq.1.0)then
        bt=0.0
      else
        bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x))
      endif
      if(x.lt.(a+1.0)/(a+b+2.0))then
        betai=bt*betacf(a,b,x)/a
        return
      else
        betai=1.0-bt*betacf(b,a,1.0-x)/b
        return
      endif
      
	  end function betai



	  ! Continued fraction evaluation.
	  ! Used by routine betai.
	  ! Numerical Recipes, chapter 6.4.
      real function betacf(a,b,x)
      integer, parameter	:: MAXIT = 100
      real, intent(in)		:: a,b,x
      real, parameter		:: EPS=3.0e-7,FPMIN=1.0e-30
      integer				:: m,m2
      real					:: aa,c,d,del,h,qab,qam,qap
      qab=a+b
      qap=a+1.0
      qam=a-1.0
      c=1.0
      d=1.0-qab*x/qap
      if(abs(d).lt.FPMIN)d=FPMIN
      d=1.0/d
      h=d
      do 11 m=1,MAXIT
        m2=2*m
        aa=m*(b-m)*x/((qam+m2)*(a+m2))
        d=1.0+aa*d
        if(abs(d).lt.FPMIN)d=FPMIN
        c=1.0+aa/c
        if(abs(c).lt.FPMIN)c=FPMIN
        d=1.0/d
        h=h*d*c
        aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2))
        d=1.0+aa*d
        if(abs(d).lt.FPMIN)d=FPMIN
        c=1.0+aa/c
        if(abs(c).lt.FPMIN)c=FPMIN
        d=1.0/d
        del=d*c
        h=h*del
        if(abs(del-1.0).lt.EPS)goto 1
11    continue
      write(*,*) 'a or b too big, or MAXIT too small in betacf'
1     betacf=h
      return
      
	  end function betacf



	  ! Logarithm of gamma function.
	  ! Used by routine betai.
	  ! Numerical Recipes, chapter 6.1.
      real function gammln(xx)
      real					:: xx
      integer				:: j
      real, dimension(6)	:: cof(6)
	  real					:: ser,stp,tmp,x,y
      save cof,stp
      data cof,stp/76.18009172947146,-86.50532032941677,  &
      24.01409824083091,-1.231739572450155,.1208650973866179e-2,  &
      -0.5395239384953e-5,2.5066282746310005/
      x=xx
      y=x
      tmp=x+5.5
      tmp=(x+0.5)*log(tmp)-tmp
      ser=1.000000000190015
      do 11 j=1,6
        y=y+1.0
        ser=ser+cof(j)/y
11    continue
      gammln=tmp+log(stp*ser/x)
      return
      
	  end function gammln