      subroutine test_t(nsamp1,samp1,nsamp2,samp2,alpha,
     &ntail,fmean1,std1,fmean2,std2,isitdifferent)
      implicit none
!----------------Inputs--------------------------
!nsamp: the number of samples
!samp: the sample values
!alpha: the significance level to consider (e.g. 0.05. 0.01)
!ntail: ntail = 1, one tail (one mean is significantly larger (smaller) than the other mean
!       ntail = 2, two tails (the two means are significantly different from each other)
!----------------Outputs---------------------------
!fmean and std: sample means and standard deviations
!isitdifferent=0, not different
!             =1, different
      integer nsamp1,nsamp2,ntail,isitdifferent,ndegfree
      double precision samp1(nsamp1),samp2(nsamp2),alpha,
     &fmean1,std1,fmean2,std2,Sign_Level,student_t,t,t0,s0

      call stdmean(nsamp1,samp1,std1,fmean1)
      call stdmean(nsamp2,samp2,std2,fmean2)
      ndegfree=nsamp1+nsamp2-2
      s0=dble((nsamp1-1))*std1*std1+
     &dble((nsamp2-1))*std2*std2
      s0=dsqrt(s0/dble(ndegfree))
      t=(fmean1-fmean2)/
     &(s0*dsqrt(1.0d0/dble(nsamp1)+1.0d0/dble(nsamp2)))      
      if(ntail.eq.1)then
        Sign_Level=1.0d0-alpha*2.0d0
      else
        Sign_Level=1.0d0-alpha/2.0d0
      endif
      t0=student_t(ndegfree,Sign_Level)
      isitdifferent=0
      if(dabs(t).gt.t0)isitdifferent=1
      return
      end
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      double precision function student_t(ndegfree,Sign_Level)
!
! the integration from -student_t to student_t =Sign_Level
!
! The student-t is calculated for a given degree of freedom
! and at a certain significance level.
! 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 solve the above equation for x (or student_t).
! Routines from Numerical Recipes are used for that.

      implicit none
! Input variables.
      integer ndegfree
! Degree of freedom
      double precision Sign_Level
! Significance level

! Functions and parameters.
      double precision zbrent,tobesolved
      double precision x1,x2,b,eps
      parameter(x1=0.0d0,x2=1.0d0,b=0.50d0,eps=1.0d-7)
  
!  Various parameters: x1, x2 bracket the root, given with 
!  accuracy eps.
           
! Local
      double precision Df,a
! Degrees of freedom
  
! a = 0.50 * Df

      external zbrent,tobesolved
     
      Df = dble(ndegfree)
      a = 0.50d0 * Df
      student_t=zbrent(tobesolved,a,b,Sign_Level,x1,x2,eps)
      student_t = dsqrt( Df/student_t  - Df)
      end function student_t


      double precision function tobesolved( a, b, c, x ) 

      implicit none
      double precision a, b, c, x
! a, b, c: parameters to the function
! x:       variable
      double precision betai
      external betai
! Incomplete beta function.
      tobesolved = betai(a,b,x) - 1.0d0 + 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. 
     
      double precision function zbrent(func,aaa,bbb,ccc,x1,x2,tol)
      implicit none
      integer ITMAX,iter
      double precision tol,x1,x2,func,EPS,
     &       a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm,
     &       aaa,bbb,ccc
      parameter(ITMAX=5000)
      parameter(EPS=3.0d-8)
      external func
      
      a=x1
      b=x2
      
      fa=func(aaa,bbb,ccc,a)
      fb=func(aaa,bbb,ccc,b)
      if((fa.gt.0.0d0.and.fb.gt.0.0d0).or.
     &   (fa.lt.0.0d0.and.fb.lt.0.0d0))then
	write(*,*) 'root must be bracketed for zbrent'
      endif
      c=b
      fc=fb
      do 11 iter=1,ITMAX
        if((fb.gt.0.0d0.and.fc.gt.0.0d0).or.
     &      (fb.lt.0.0d0.and.fc.lt.0.0d0))then
          c=a
          fc=fa
          d=b-a
          e=d
        endif
        if(dabs(fc).lt.dabs(fb)) then
          a=b
          b=c
          c=a
          fa=fb
          fb=fc
          fc=fa
        endif
        tol1=2.0d0*EPS*dabs(b)+0.5d0*tol
        xm=0.5d0*(c-b)
        if(dabs(xm).le.tol1.or.fb.eq.0.0d0)then
          zbrent=b
          return
        endif
        if(dabs(e).ge.tol1.and.dabs(fa).gt.dabs(fb))then
          s=fb/fa
          if(a.eq.c) then
            p=2.0d0*xm*s
            q=1.0d0-s
          else
            q=fa/fc
            r=fb/fc
            p=s*(2.0d0*xm*q*(q-r)-(b-a)*(r-1.0d0))
            q=(q-1.0d0)*(r-1.0d0)*(s-1.0d0)
          endif
          if(p.gt.0.0d0)q=-q
          p=dabs(p)
          if(2.0d0*p.lt.dmin1(3.0d0*xm*q-dabs(tol1*q),dabs(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(dabs(d).gt.tol1)then
          b=b+d
        else
          b=b+dsign(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.
      double precision function betai(a,b,x)
      double precision a,b,x
!U    USES betacf,gammln
      double precision bt
      double precision betacf,gammln
      external betacf,gammln

      if(x.lt.0.0d0.or.x.gt.1.0d0)then
        write(*,*) 'bad argument x in betai'
      endif
      if(x.eq.0.0d0.or.x.eq.1.0d0)then
        bt=0.0d0
      else
        bt=dexp(gammln(a+b)-gammln(a)-gammln(b)+
     &          a*dlog(x)+b*dlog(1.0d0-x))
      endif
      if(x.lt.(a+1.0d0)/(a+b+2.0d0))then
        betai=bt*betacf(a,b,x)/a
        return
      else
        betai=1.0d0-bt*betacf(b,a,1.0d0-x)/b
        return
      endif
      end function betai


! Continued fraction evaluation.
! Used by routine betai.
 ! Numerical Recipes, chapter 6.4.
      double precision function betacf(a,b,x)
      implicit none
      integer MAXIT,m,m2
      double precision a,b,x,EPS,FPMIN
      double precision aa,c,d,del,h,qab,qam,qap
      parameter(MAXIT = 100)
      parameter(EPS=3.0d-7,FPMIN=1.0d-30)

      qab=a+b
      qap=a+1.0d0
      qam=a-1.0d0
      c=1.0d0
      d=1.0d0-qab*x/qap
      if(dabs(d).lt.FPMIN)d=FPMIN
      d=1.0d0/d
      h=d
      do 11 m=1,MAXIT
        m2=2*m
        aa=dble(m)*(b-dble(m))*x/((qam+dble(m2))*(a+dble(m2)))
        d=1.0d0+aa*d
        if(dabs(d).lt.FPMIN)d=FPMIN
        c=1.0d0+aa/c
        if(dabs(c).lt.FPMIN)c=FPMIN
        d=1.0d0/d
        h=h*d*c
        aa=-(a+dble(m))*(qab+dble(m))*x/((a+dble(m2))*(qap+dble(m2)))
        d=1.0d0+aa*d
        if(dabs(d).lt.FPMIN)d=FPMIN
        c=1.0d0+aa/c
        if(dabs(c).lt.FPMIN)c=FPMIN
        d=1.0d0/d
        del=d*c
        h=h*del
        if(dabs(del-1.0d0).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
