! Collected Special Functions

!----------------------------------------
FUNCTION betai(a,b,x)
REAL betai,a,b,x
REAL bt,betacf,gammln

if(x.lt.0..or.x.gt.1.)pause "bad argument x in betai"
if(x.eq.0..or.x.eq.1.)then
   bt=0.
else
   bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.-x))
endif
if (x.lt.(a+1.)/(a+b+2.))then
   betai=bt*betacf(a,b,x)/a
   return
else
   betai=1.-bt*betacf(b,a,1.-x)/b
   return
endif
END FUNCTION betai
!------------------------------------
FUNCTION betacf(a,b,x)
INTEGER MAXIT
REAL betacf,a,b,x,EPS,FPMIN
PARAMETER (MAXIT=100,EPS=3.e-7,FPMIN=1.e-30)
 INTEGER m,m2
REAL aa,c,d,del,h,qab,qam,qap
qab=a+b
qap=a+1.
qam=a-1.
c=1.
d=1.-qab*x/qap
if(abs(d).lt.FPMIN)d=FPMIN
d=1./d
h=d
Do m=1,MAXIT
   m2=2*m
   aa=m*(b-m)*x/((qam+m2)*(a+m2))
   d=1.+aa*d
   if(abs(d).lt.FPMIN)d=FPMIN
   c=1.+aa/c
   if(abs(c).lt.FPMIN)c=FPMIN
   d=1./d
   h=h*d*c
   aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2))
   d=1.+aa*d
   if(abs(d).lt.FPMIN)d=FPMIN
   c=1.+aa/c
   if(abs(c).lt.FPMIN)c=FPMIN
   d=1./d
   del=d*c
   h=h*del
   if(abs(del-1.).lt.EPS) goto 1
Enddo
1  betacf=h
return
END FUNCTION betacf
!------------------------------------
FUNCTION gammln(xx)
REAL gammln,xx
INTEGER j
DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
SAVE cof,stp
DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, &
     24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, &
     -.5395239384953d-5,2.5066282746310005d0/
x=xx
y=x
tmp=x+5.5d0
tmp=(x+0.5d0)*log(tmp)-tmp
ser=1.000000000190015d0
Do j=1,6
   y=y+1.d0
   ser=ser+cof(j)/y
Enddo
gammln=tmp+log(stp*ser/x)
return
END FUNCTION gammln
!-------------------------------------------
Real FUNCTION FINV(P,N1,N2)
Real, PARAMETER :: EPS=1.0E-6
REAL :: P,AP,N1,N2, A,B,C,P0,SS,SS1
IF (P.LT.0.0.OR.P.GT.1.0) THEN
   WRITE(*,*)' INVALID NUMERIC INPUT IN FINV FUNCTION'
   STOP
ENDIF
IF (P.EQ.0.0) THEN
   WRITE(*,*)' THE P VALUE TOO SMALL IN FINV FUNCTION'
   STOP
ELSE IF (P.EQ.1) THEN
   FINV=0.0
   RETURN
ENDIF
!
A=0.0
B=9999999.0
C=A
AP=P
 10 P0=FDIST(B,N1,N2)
    SS=ABS((P0-AP)/AP)
    SS1=ABS((B-C)/B)
    IF (SS.GE.EPS.AND.SS1.GE.EPS) THEN
       IF (AP.GT.P0) THEN
          C=B
          B=(A+B)/2.0
          GOTO 10
       ELSE IF (AP.LT.P0) THEN
          B=(C+B)/2.0
          GOTO 10
       ENDIF
    ENDIF
    A=(C+B)/2.0
    FINV=A
END FUNCTION FINV
!-------------------------------------------
Real FUNCTION FDIST(X0,N1,N2)
REAL :: X0,N1,N2,X
X=N2/(N2+N1*X0)
FDIST=BETAI(0.5*N2,0.5*N1,X)
END FUNCTION FDIST
!-------------------------------------------
Real FUNCTION erfcc (x)
REAL :: x, tmp
REAL :: t,z
z = abs(x)
t = 1./(1.+0.5*z)
 tmp  = t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+       &
        t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+ &
        t*(1.48851587+t*(-.82215223+t*.17087277)))))))))
if (x.lt.0.) tmp=2.-tmp
erfcc=tmp
END FUNCTION erfcc
!-----------------------------------------------------------------------
      Real function percentile(x, length, per, RMiss)
      integer length
      real x(length), per, RMiss
      real xtos(length),bb,cc,oout
      integer nn, i

      nn=0
      do i=1, length
        if(x(i) /= RMiss)then
          nn=nn+1
          xtos(nn)=x(i)
        endif
      enddo

      if(nn.eq.0) then
        oout=RMiss
      else
        call sort(nn,xtos)
          bb=nn*per+per/3.+1/3.
          cc=real(int(bb))
          if(int(cc).ge.nn) then
            oout=xtos(nn)
          else
            oout=xtos(int(cc))+(bb-cc)*(xtos(int(cc)+1)-xtos(int(cc)))
          endif
      endif
      percentile = oout
      end function percentile
!-----------------------------------------------------------------------
      SUBROUTINE sort(n,arr)
      INTEGER n,M,NSTACK
      REAL arr(n)
      PARAMETER (M=7,NSTACK=50)
      INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
      REAL a,temp
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
        do 12 j=l+1,ir
          a=arr(j)
          do 11 i=j-1,1,-1
            if(arr(i).le.a)goto 2
            arr(i+1)=arr(i)
11        continue
          i=0
2         arr(i+1)=a
12      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        temp=arr(k)
        arr(k)=arr(l+1)
        arr(l+1)=temp
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          arr(l+1)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l).gt.arr(ir))then
          temp=arr(l)
          arr(l)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l+1).gt.arr(l))then
          temp=arr(l+1)
          arr(l+1)=arr(l)
          arr(l)=temp
        endif
        i=l+1
        j=ir
        a=arr(l)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        arr(i)=arr(j)
        arr(j)=temp
        goto 3
5       arr(l)=arr(j)
        arr(j)=a
        jstack=jstack+2
        if(jstack.gt.NSTACK)pause 'NSTACK too small in sort'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1
      END SUBROUTINE sort
!--------------------------End of this file--------------------