! *********************************************************
      FUNCTION bessi(n,x)
      implicit double precision (a-h,o-z)
      INTEGER n,IACC
      PARAMETER (IACC=40,BIGNO=1.0e10,BIGNI=1.0e-10)
CU    USES bessi0
      INTEGER j,m
!
      if (n.lt.2) stop 'bad argument n in bessi'
! ---------------------------------------------------------
!  RCF modified  31 Oct 00; G. Penn 30 Jun 01
!      if (x.eq.0.) then
!      if( ABS(x) .lt. 1.e-10 ) then
      if( (ABS(x)/2.)**n .lt. 1.e-20 ) then
! ---------------------------------------------------------
        bessi=0.
      else
        tox=2.0/abs(x)
        bip=0.0
        bi=1.0
        bessi=0.
        m=2*((n+int(sqrt(float(IACC*n)))))
        do 11 j=m,1,-1
          bim=bip+float(j)*tox*bi
          bip=bi
          bi=bim
          if (abs(bi).gt.BIGNO) then
            bessi=bessi*BIGNI
            bi=bi*BIGNI
            bip=bip*BIGNI
          endif
          if (j.eq.n) bessi=bip
11      continue
        bessi=bessi*bessi0(x)/bi
        if (x.lt.0..and.mod(n,2).eq.1) bessi=-bessi
      endif
      return
      END
! *********************************************************
      FUNCTION bessi0(x)
      implicit double precision (a-h,o-z)
      common/FLAG/iflag   ! RCF modified 10 MAR 01
      SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
      DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,3.5156229d0,3.0899424d0,
     *1.2067492d0,0.2659732d0,0.360768d-1,0.45813d-2/
      DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,0.1328592d-1,
     *0.225319d-2,-0.157565d-2,0.916281d-2,-0.2057706d-1,0.2635537d-1,
     *-0.1647633d-1,0.392377d-2/
! ---------------------------------------------------------
!  RCF modified  31 Oct 00; 10 Mar 01
      if( ABS(x) .lt. 1.e-10 ) then
        bessi0 = 1.
        return
      else if( ABS(x) .gt. 100. ) then
        bessi0 = 0.
        iflag = -74
        return
      end if
! ---------------------------------------------------------
      if (abs(x).lt.3.75) then
        y=(x/3.75)**2
        bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))
      else
        ax=abs(x)
        y=3.75/ax
        bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
     *(q7+y*(q8+y*q9))))))))
      endif
      return
      END
c  ********************************************************
      FUNCTION bessi1(x)
      implicit double precision (a-h,o-z)
      common/FLAG/iflag   ! RCF modified 10 Mar 01
      SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
      DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,
     *0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/
      DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,-0.3988024d-1,
     *-0.362018d-2,0.163801d-2,-0.1031555d-1,0.2282967d-1,-0.2895312d-1,
     *0.1787654d-1,-0.420059d-2/
! ---------------------------------------------------------
!  RCF modified  31 Oct 00; 10 Mar 01; 26 Apr 01
      if( ABS(x) .lt. 1.e-10 ) then
        bessi1 = 0.
        return
      else if( ABS(x) .gt. 100. ) then
        bessi1 = 0.
        iflag = -74
        return
      end if
! ---------------------------------------------------------
      if (abs(x).lt.3.75) then
        y=(x/3.75)**2
        bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
      else
        ax=abs(x)
        y=3.75/ax
        bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
     *(q7+y*(q8+y*q9))))))))
        if(x.lt.0.)bessi1=-bessi1
      endif
      return
      END
c  ********************************************************
      FUNCTION bessj0(x)
      implicit double precision (a-h,o-z)
      SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
     *s5,s6
      DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
     *-.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1,
     *.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
      DATA r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,
     *651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/,s1,s2,
     *s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0,
     *59272.64853d0,267.8532712d0,1.d0/
!
      if(abs(x).lt.8.)then
        y=x**2
        bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y*
     *(s4+y*(s5+y*s6)))))
      else
        ax=abs(x)
        z=8./ax
        y=z**2
        xx=ax-.785398164
        bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
     *p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
      endif
      return
      END
c  ********************************************************
      FUNCTION bessj1(x)
      implicit double precision (a-h,o-z)
      SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
     *s5,s6
      DATA r1,r2,r3,r4,r5,r6/72362614232.d0,-7895059235.d0,
     *242396853.1d0,-2972611.439d0,15704.48260d0,-30.16036606d0/,s1,s2,
     *s3,s4,s5,s6/144725228442.d0,2300535178.d0,18583304.74d0,
     *99447.43394d0,376.9991397d0,1.d0/
      DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4,
     *.2457520174d-5,-.240337019d-6/, q1,q2,q3,q4,q5/.04687499995d0,
     *-.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/
!
      if(abs(x).lt.8.)then
        y=x**2
        bessj1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+
     *y*(s4+y*(s5+y*s6)))))
      else
        ax=abs(x)
        z=8./ax
        y=z**2
        xx=ax-2.356194491
        bessj1=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
     *p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))*sign(1.d0,x)
      endif
      return
      END
c  ********************************************************
      FUNCTION bessk(n,x)
      implicit double precision (a-h,o-z)
      INTEGER n
CU    USES bessk0,bessk1
      INTEGER j
!
      if (n.lt.2) stop 'bad argument n in bessk'
      tox=2.0/x
      bkm=bessk0(x)
      bk=bessk1(x)
      do 11 j=1,n-1
        bkp=bkm+j*tox*bk
        bkm=bk
        bk=bkp
11    continue
      bessk=bk
      return
      END
c  ********************************************************
      FUNCTION bessk0(x)
      implicit double precision (a-h,o-z)
      common/FLAG/iflag   ! RCF modified 10 MAR 01
CU    USES bessi0
      DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y
      SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7
      DATA p1,p2,p3,p4,p5,p6,p7/-0.57721566d0,0.42278420d0,0.23069756d0,
     *0.3488590d-1,0.262698d-2,0.10750d-3,0.74d-5/
      DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,-0.7832358d-1,0.2189568d-1,
     *-0.1062446d-1,0.587872d-2,-0.251540d-2,0.53208d-3/
! ---------------------------------------------------------
!  RCF modified  10 Mar 01
      if( x .lt. 1.e-3 ) then
        bessk0 = 0.
        iflag = -75
        return
      else if( x .gt. 100. ) then
        bessk0 = 0.
        return
      end if
! ---------------------------------------------------------
      if (x.le.2.0) then
        y=x*x/4.0
        bessk0=(-log(x/2.0)*bessi0(x))+(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*
     *(p6+y*p7))))))
      else
        y=(2.0/x)
        bessk0=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
     *q7))))))
      endif
      return
      END
c  ********************************************************
      FUNCTION bessk1(x)
      implicit double precision (a-h,o-z)
      common/FLAG/iflag   ! RCF modified 10 MAR 01
CU    USES bessi1
      DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y
      SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7
      DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,0.15443144d0,-0.67278579d0,
     *-0.18156897d0,-0.1919402d-1,-0.110404d-2,-0.4686d-4/
      DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,0.23498619d0,-0.3655620d-1,
     *0.1504268d-1,-0.780353d-2,0.325614d-2,-0.68245d-3/
! ---------------------------------------------------------
!  RCF modified  10 Mar 01
      if( x .lt. 1.e-3 ) then
        bessk1 = 0.
        iflag = -75
        return
      else if( x .gt. 100. ) then
        bessk1 = 0.
        return
      end if
! ---------------------------------------------------------
      if (x.le.2.0) then
        y=x*x/4.0
        bessk1=(log(x/2.0)*bessi1(x))+(1.0/x)*(p1+y*(p2+y*(p3+y*(p4+y*
     *(p5+y*(p6+y*p7))))))
      else
        y=2.0/x
        bessk1=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
     *q7))))))
      endif
      return
      end
c  ********************************************************
      FUNCTION bessy0(x)
      implicit double precision (a-h,o-z)
CU    USES bessj0
      SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
     *s5,s6
      DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
     *-.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1,
     *.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
      DATA r1,r2,r3,r4,r5,r6/-2957821389.d0,7062834065.d0,
     *-512359803.6d0,10879881.29d0,-86327.92757d0,228.4622733d0/,s1,s2,
     *s3,s4,s5,s6/40076544269.d0,745249964.8d0,7189466.438d0,
     *47447.26470d0,226.1030244d0,1.d0/
      if(x.lt.8.)then
        y=x**2
        bessy0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y*
     *(s4+y*(s5+y*s6)))))+.636619772*bessj0(x)*log(x)
      else
        z=8./x
        y=z**2
        xx=x-.785398164
        bessy0=sqrt(.636619772/x)*(sin(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
     *p5))))+z*cos(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
      endif
      return
      END
c  ********************************************************
      FUNCTION bessy1(x)
      implicit double precision (a-h,o-z)
CU    USES bessj1
      SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
     *s5,s6,s7
      DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4,
     *.2457520174d-5,-.240337019d-6/, q1,q2,q3,q4,q5/.04687499995d0,
     *-.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/
      DATA r1,r2,r3,r4,r5,r6/-.4900604943d13,.1275274390d13,
     *-.5153438139d11,.7349264551d9,-.4237922726d7,.8511937935d4/,s1,s2,
     *s3,s4,s5,s6,s7/.2499580570d14,.4244419664d12,.3733650367d10,
     *.2245904002d8,.1020426050d6,.3549632885d3,1.d0/
      if(x.lt.8.)then
        y=x**2
        bessy1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+
     *y*(s4+y*(s5+y*(s6+y*s7))))))+.636619772*(bessj1(x)*log(x)-1./x)
      else
        z=8./x
        y=z**2
        xx=x-2.356194491
        bessy1=sqrt(.636619772/x)*(sin(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
     *p5))))+z*cos(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
      endif
      return
      END
! *********************************************************
C
C     THREE DIMENSIONAL BSPLINE ROUTINES
C
c     Report all suspected bugs to J. Scott Berg <jsberg@bnl.gov>
c     or <jsberg@mailaps.org>
c
c     $Date: 2003/03/28 17:51:07 $
c     $Revision: 1.3 $
c     $Source: /home/jsberg/cvsroot/bsplinef/bsplinef3.f,v $
c
C     THESE ARE LIKE THE ONE-DIMENSIONAL ROUTINES, EXCEPT THAT THREE
C     SETS OF ABCISSAS ARE SUPPLIED, NX, N, AND K ARE THREE-ELEMENT
C     ARRAYS, AND THE VALUES TO INTERPOLATE ARE SUPPLIED IN AN ARRAY
C     WITH THE FIRST VARIABLE VARYING MOST RAPIDLY.
C     
C     DOUBLE PRECISION FUNCTION BS3EV(X,C,N,K,T,S)
C     
C     INPUTS:
C
C     X(3)              (DP) = POINT AT WHICH TO EVALUATE SPLINE SERIES
C     C(N(1)*N(2)*N(3)) (DP) = SPLINE COEFFICIENTS
C     N(3)              (I)  = NUMBERS OF SPLINE BASIS FUNCTIONS
C     K(3)              (I)  = SPLINE ORDERS
C     T(N(1)+K(1)+N(2)+K(2)+N(3)+K(3)) (DP) = KNOT SEQUENCES
C     S(K(2)*K(3)+K(1)) (DP) = SCRATCH ARRAY
C
! *********************************************************
      DOUBLE PRECISION FUNCTION BS3EV(X,C,N,K,T,S)
C
C     EVALUATE INTERPOLATION OVER 3 DIMENSIONS.
C     X IS POINT TO EVALUATE AT
C     C IS COEFFICIENT ARRAY
C     N ARE THE NUMBER OF COEFFICIENTS IN EACH DIMENSION
C     K IS BSPLINE ORDER
C     T ARE KNOT SEQUENCES, STACKED SEQUENTIALLY
C     S IS SCRATCH ARRAY, SIZE K(2)*K(3)+K(1)
C

      IMPLICIT NONE
      INTEGER K(3),N(3)
      DOUBLE PRECISION C(*),S(*),T(*),X(3)
      INTEGER I1,I2,I3,J1,J2,J3,KK,NS,BSIVL,IC

      I1 = BSIVL(X(1),T(1),N(1),K(1))
      I2 = BSIVL(X(2),T(N(1)+K(1)+1),N(2),K(2))
      I3 = BSIVL(X(3),T(N(1)+K(1)+N(2)+K(2)+1),N(3),K(3))
      NS = K(2)*K(3)

      DO 100 J3=1,K(3)
         DO 100 J2=1,K(2)
            IC=N(1)*(I2+J2-K(2)-1+N(2)*(I3+J3-K(3)-1))
            DO 120 J1=1,K(1)
 120           S(J1+NS)=C(I1+J1-K(1)+IC)
            DO 130 KK=K(1)-1,1,-1
               DO 130 J1=I1,I1-KK+1,-1
 130              S(J1-I1+K(1)+NS) = (S(J1-I1+K(1)+NS)*(X(1)-T(J1))+
     $                 S(J1-I1+K(1)-1+NS)*(T(J1+KK)-X(1)))/
     $                 (T(J1+KK)-T(J1))
 100        S(J2+K(2)*(J3-1)) = S(K(1)+NS)

      NS = N(1)+K(1)

      DO 200 J3=1,K(3)
         J1=K(2)*J3-I2
         DO 200 KK=K(2)-1,1,-1
            DO 200 J2=I2,I2-KK+1,-1
 200           S(J2+J1) = (S(J2+J1)*(X(2)-T(J2+NS))+
     $              S(J2+J1-1)*(T(J2+KK+NS)-X(2)))
     $              /(T(J2+KK+NS)-T(J2+NS))

      NS = NS+N(2)+K(2)
      J2 = K(3)-I3

      DO 300 KK=K(3)-1,1,-1
         DO 300 J3=I3,I3-KK+1,-1
 300        S((J3+J2)*K(2)) = (S((J3+J2)*K(2))*(X(3)-T(J3+NS))+
     $           S((J3+J2-1)*K(2))*(T(J3+KK+NS)-X(3)))/
     $           (T(J3+KK+NS)-T(J3+NS))

      BS3EV = S(K(2)*K(3))

      END
! *********************************************************
      INTEGER FUNCTION BSIVL(X,T,N,K)
C
C     RETURN I SUCH THAT T(I) <= X < T(I+1)
C

      IMPLICIT NONE
      INTEGER K,N
      DOUBLE PRECISION X,T(N+K)
      INTEGER I,J,L

      I = K
      J = N+1

 100  IF (I+1.EQ.J) GOTO 300
      L = (I+J)/2
      IF (X.LT.T(L)) GOTO 200
      I = L
      GOTO 100
 200  J = L
      GOTO 100

 300  BSIVL = I
      END
c  ********************************************************
      FUNCTION dfridr(func,x,h,err)
      implicit double precision (a-h,o-z)
      INTEGER NTAB
      PARAMETER (CON=1.4,CON2=CON*CON,BIG=1.E30,NTAB=10,SAFE=2.)
      EXTERNAL func
CU    USES func
      INTEGER i,j
      REAL*8 a(NTAB,NTAB)
!
      if(h.eq.0.) stop 'h must be nonzero in dfridr'
      hh=h
      a(1,1)=(func(x+hh)-func(x-hh))/(2.0*hh)
      err=BIG
      do 12 i=2,NTAB
        hh=hh/CON
        a(1,i)=(func(x+hh)-func(x-hh))/(2.0*hh)
        fac=CON2
        do 11 j=2,i
          a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1.)
          fac=CON2*fac
          errt=max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
          if (errt.le.err) then
            err=errt
            dfridr=a(j,i)
          endif
11      continue
        if(abs(a(i,i)-a(i-1,i-1)).ge.SAFE*err)return
12    continue
      return
      END
c  ********************************************************
      FUNCTION elle(phi,ak)
      implicit double precision (a-h,o-z)
CU    USES rd,rf
!
      s=sin(phi)
      cc=cos(phi)**2
      q=(1.-s*ak)*(1.+s*ak)
      elle=s*(rf(cc,q,1.d0)-((s*ak)**2)*rd(cc,q,1.d0)/3.)
      return
      END
c  ********************************************************
      FUNCTION ellf(phi,ak)
      implicit double precision (a-h,o-z)
CU    USES rf
!
      s=sin(phi)
      ellf=s*rf(cos(phi)**2,(1.-s*ak)*(1.+s*ak),1.d0)
      return
      END
c  ********************************************************
      FUNCTION ellpi(phi,en,ak)
      implicit double precision (a-h,o-z)
CU    USES rf,rj
!
      s=sin(phi)
      enss=en*s*s
      cc=cos(phi)**2
      q=(1.-s*ak)*(1.+s*ak)
      ellpi=s*(rf(cc,q,1.d0)-enss*rj(cc,q,1.d0,1.+enss)/3.)
      return
      END
c  ********************************************************
      FUNCTION factrl(n)
      implicit double precision (a-h,o-z)
      INTEGER n
CU    USES gammln
      INTEGER j,ntop
      REAL*8 a(33)
      SAVE ntop,a
      DATA ntop,a(1)/0,1./
!
      if (n.lt.0) then
        stop 'negative factorial in factrl'
      else if (n.le.ntop) then
        factrl=a(n+1)
      else if (n.le.32) then
        do 11 j=ntop+1,n
          a(j+1)=j*a(j)
11      continue
        ntop=n
        factrl=a(n+1)
      else
        factrl=exp(gammln(n+1.d0))
      endif
      return
      END
c  ********************************************************
      FUNCTION gammln(xx)
      implicit double precision (a-h,o-z)
      INTEGER j
      real*8 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 11 j=1,6
        y=y+1.d0
        ser=ser+cof(j)/y
11    continue
      gammln=tmp+log(stp*ser/x)
      return
      END
c  ********************************************************
      FUNCTION gasdev(idum)
      implicit double precision (a-h,o-z)
      INTEGER idum
CU    USES ran1
      INTEGER iset
      SAVE iset,gset
      DATA iset/0/
!
      if (iset.eq.0) then
1       v1=2.*ran1(idum)-1.
        v2=2.*ran1(idum)-1.
        rsq=v1**2+v2**2
        if(rsq.ge.1..or.rsq.eq.0.)goto 1
        fac=sqrt(-2.*log(rsq)/rsq)
        gset=v1*fac
        gasdev=v2*fac
        iset=1
      else
        gasdev=gset
        iset=0
      endif
      return
      END
c  ********************************************************
      SUBROUTINE hunt(xx,n,x,jlo)
      implicit double precision (a-h,o-z)
      INTEGER jlo,n
      REAL*8 xx(n)
      INTEGER inc,jhi,jm
      LOGICAL ascnd
!
      ascnd=xx(n).gt.xx(1)
      if(jlo.le.0.or.jlo.gt.n)then
        jlo=0
        jhi=n+1
        goto 3
      endif
      inc=1
      if(x.ge.xx(jlo).eqv.ascnd)then
1       jhi=jlo+inc
        if(jhi.gt.n)then
          jhi=n+1
        else if(x.ge.xx(jhi).eqv.ascnd)then
          jlo=jhi
          inc=inc+inc
          goto 1
        endif
      else
        jhi=jlo
2       jlo=jhi-inc
        if(jlo.lt.1)then
          jlo=0
        else if(x.lt.xx(jlo).eqv.ascnd)then
          jhi=jlo
          inc=inc+inc
          goto 2
        endif
      endif
3     if(jhi-jlo.eq.1)return
      jm=(jhi+jlo)/2
      if(x.gt.xx(jm).eqv.ascnd)then
        jlo=jm
      else
        jhi=jm
      endif
      goto 3
      END
c  ********************************************************
      integer function IMPMID(x0,v)
C
C Implicit midpoint integrator.  Assumes step size 1.
C Takes a 6-D vector x0 which is modified, and a subroutine v
C The subroutine v is the vector field, and takes two arguments:
C a 6-D vector where the vector field is evaluated, and a 6-D
C vector in which the vector field is returned.
C Returns 0 on success, 1 if we don't converge in maxit iterations
C
!  Ref.  J.S. Berg            22 Nov 2004
!
      double precision x0(6)
      double precision e0(6),e1(6),v0(6),v1(6),x1(6)
      logical b(6),ball
      integer i,j,maxit,rc,v
      parameter (maxit=100)
      external v

      rc = v(x0,v0)
      if (rc.ne.0) then
         impmid = rc
         return
      end if
      do 100 i=1,6
         x1(i) = x0(i)+v0(i)
         b(i) = .false.
 100  continue
      rc = v(x1,v1)
      if (rc.ne.0) then
         impmid = rc
         return
      end if
      do 200 i=1,6
         e1(i) = abs(v0(i))+abs(v1(i))
         v1(i) = 0.5d0*(v0(i)+v1(i))
 200  continue

      j = 0
 300  do 400 i=1,6
         x1(i) = x0(i)+0.5d0*v1(i)
         v0(i) = v1(i)
         e0(i) = e1(i)
 400  continue
      rc = v(x1,v1)
      if (rc.ne.0) then
         impmid = rc
         return
      end if
      ball = .true.
      do 500 i=1,6
         e1(i) = abs(v1(i)-v0(i))
         b(i) = b(i) .or. e1(i).eq.0d0 .or. e1(i).ge.e0(i)
         ball = ball .and. b(i)
 500  continue
      j = j + 1
      if (.not.ball.and.j.lt.maxit) goto 300

      if (j.lt.maxit) then
         impmid = 0
      else
         impmid = -1
      end if

      do 600 i=1,6
         x0(i) = x0(i) + v1(i)
 600  continue

      end
!  ********************************************************
      subroutine LOCATE2(xx,n,x,j)
!
! Locate position J of variable X in the array XX(1:N)
!  Separate array into regions [1) [2) [3) ... [N)
!  x=xx(1) => j=1          x=xx(n) => j=n
!  Return j=0 if x is out of the array range
!
      implicit real*8(a-h,o-z)
      real*8 xx(n)
!
      if( x.lt.xx(1) .or. x.gt.xx(n) ) then
        j = 0
        go to 900
      end if
!
      jl = 0            ! lower bracket
      jh = n + 1        ! upper bracket
!
10    continue
      if( jh-jl .gt. 1 ) then       ! continue working
        jm = (jl + jh) / 2          ! bisect the bracket
        if( x .ge. xx(jm) ) then
          jl = jm
        else
          jh = jm
        end if
        go to 10
      end if
!
      j = jl
!
900   continue
      end
c  ********************************************************
      subroutine LUDCMP(a,n,np,indx,d)
! 
! Make LU decomposition of matrix
!
      implicit double precision (a-h,o-z)
      INTEGER n,np,indx(n),NMAX
      REAL*8 a(np,np)
      PARAMETER (NMAX=500,TINY=1.0e-20)
      INTEGER i,imax,j,k
      REAL*8 vv(NMAX)
!
      d=1.
      do 12 i=1,n
        aamax=0.
        do 11 j=1,n
          if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
11      continue
        if (aamax.eq.0.) then ! singular matrix in ludcmp
          d = 0.
          go to 900
        end if
        vv(i)=1./aamax
12    continue
      do 19 j=1,n
        do 14 i=1,j-1
          sum=a(i,j)
          do 13 k=1,i-1
            sum=sum-a(i,k)*a(k,j)
13        continue
          a(i,j)=sum
14      continue
        aamax=0.
        do 16 i=j,n
          sum=a(i,j)
          do 15 k=1,j-1
            sum=sum-a(i,k)*a(k,j)
15        continue
          a(i,j)=sum
          dum=vv(i)*abs(sum)
          if (dum.ge.aamax) then
            imax=i
            aamax=dum
          endif
16      continue
        if (j.ne.imax)then
          do 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
17        continue
          d=-d
          vv(imax)=vv(j)
        endif
        indx(j)=imax
        if(a(j,j).eq.0.)a(j,j)=TINY
        if(j.ne.n)then
          dum=1./a(j,j)
          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
18        continue
        endif
19    continue
!
900   continue
      return
      END
c  ********************************************************
      FUNCTION poidev(xm,idum)
      implicit double precision (a-h,o-z)
      INTEGER idum
      PARAMETER (PI=3.141592654)
CU    USES gammln,ran1
      SAVE alxm,g,oldm,sq
      DATA oldm /-1./
!
      if (xm.lt.12.)then
        if (xm.ne.oldm) then
          oldm=xm
          g=exp(-xm)
        endif
        em=-1
        t=1.
2       em=em+1.
        t=t*ran1(idum)
        if (t.gt.g) goto 2
      else
        if (xm.ne.oldm) then
          oldm=xm
          sq=sqrt(2.*xm)
          alxm=log(xm)
          g=xm*alxm-gammln(xm+1.)
        endif
1       y=tan(PI*ran1(idum))
        em=sq*y+xm
        if (em.lt.0.) goto 1
! ---------------------------------------------------------
!  RCF modified 28 Sep 99
!   prevent large value of em (e.g. 559) from causing
!   a floating point underflow in the calculation of t, below
        if( em .gt. 30. ) go to 1
! ---------------------------------------------------------              
        em=int(em)
        t=0.9*(1.+y**2)*exp(em*alxm-gammln(em+1.)-g)
        if (ran1(idum).gt.t) goto 1
      endif
      poidev=em
      return
      END
c  ********************************************************
      SUBROUTINE polcof(xa,ya,n,cof)
      implicit double precision (a-h,o-z)
      INTEGER n,NMAX
      REAL*8 cof(n),xa(n),ya(n)
      PARAMETER (NMAX=15)
CU    USES polint
      INTEGER i,j,k
      REAL*8 x(NMAX),y(NMAX)
!
      do 11 j=1,n
        x(j)=xa(j)
        y(j)=ya(j)
11    continue
      do 14 j=1,n
        call polint(x,y,n+1-j,0.d0,cof(j),dy)
        xmin=1.e38
        k=0
        do 12 i=1,n+1-j
          if (abs(x(i)).lt.xmin)then
            xmin=abs(x(i))
            k=i
          endif
          if(x(i).ne.0.)y(i)=(y(i)-cof(j))/x(i)
12      continue
        do 13 i=k+1,n+1-j
          y(i-1)=y(i)
          x(i-1)=x(i)
13      continue
14    continue
      return
      END
c  ********************************************************
      SUBROUTINE polint(xa,ya,n,x,y,dy)
      implicit double precision (a-h,o-z)
      INTEGER n,NMAX
      REAL*8 xa(n),ya(n)
      PARAMETER (NMAX=10)
      INTEGER i,m,ns
      REAL*8 c(NMAX),d(NMAX)
!
      ns=1
      dif=abs(x-xa(1))
      do 11 i=1,n
        dift=abs(x-xa(i))
        if (dift.lt.dif) then
          ns=i
          dif=dift
        endif
        c(i)=ya(i)
        d(i)=ya(i)
11    continue
      y=ya(ns)
      ns=ns-1
      do 13 m=1,n-1
        do 12 i=1,n-m
          ho=xa(i)-x
          hp=xa(i+m)-x
          w=c(i+1)-d(i)
          den=ho-hp
          if(den.eq.0.)stop 'failure in polint'
          den=w/den
          d(i)=hp*den
          c(i)=ho*den
12      continue
        if (2*ns.lt.n-m)then
          dy=c(ns+1)
        else
          dy=d(ns)
          ns=ns-1
        endif
        y=y+dy
13    continue
      return
      END
! *********************************************************
      SUBROUTINE qgaus(func,a,b,ss)
      implicit double precision (a-h,o-z)
!      REAL a,b,ss,func
      EXTERNAL func
      INTEGER j
!      REAL dx,xm,xr,w(5),x(5)
      dimension w(5),x(5)
      SAVE w,x
      DATA w/.2955242247,.2692667193,.2190863625,.1494513491,
     *.0666713443/
      DATA x/.1488743389,.4333953941,.6794095682,.8650633666,
     *.9739065285/
      xm=0.5*(b+a)
      xr=0.5*(b-a)
      ss=0
      do 11 j=1,5
        dx=xr*x(j)
        ss=ss+w(j)*(func(xm+dx)+func(xm-dx))
11    continue
      ss=xr*ss
      return
      END
! *********************************************************
      SUBROUTINE qromb(func,a,b,ss)
      implicit double precision (a-h,o-z)
      INTEGER JMAX,JMAXP,K,KM
      EXTERNAL func
! ---------------------------------------------------------
!  RCF modified 5 Oct 00
!      PARAMETER (EPS=1.e-6, JMAX=20, JMAXP=JMAX+1, K=5, KM=K-1)
      PARAMETER (EPS=5.e-6, JMAX=20, JMAXP=JMAX+1, K=5, KM=K-1)
! ---------------------------------------------------------
CU    USES polint,trapzd
      INTEGER j
      REAL*8 h(JMAXP),s(JMAXP)
!
      h(1)=1.
      do 11 j=1,JMAX
        call trapzd(func,a,b,s(j),j)
        if (j.ge.K) then
          call polint(h(j-KM),s(j-KM),K,0.d0,ss,dss)
          if (abs(dss).le.EPS*abs(ss)) return
! ---------------------------------------------------------
!  RCF modified 5 Oct 00
! Doesn't converge for symmetric case with area ~0
          if( ABS(ss) .le. eps ) return
! ---------------------------------------------------------
        endif
        s(j+1)=s(j)
        h(j+1)=0.25*h(j)
11    continue
      stop 'too many steps in qromb'
      END
c  ********************************************************
      FUNCTION ran1(idum)
      implicit double precision (a-h,o-z)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      PARAMETER (IA=16807,IM=2147483647,AM=1.d0/IM,IQ=127773,IR=2836,
     *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=2.3d-16,RNMX=1.d0-EPS)
      INTEGER j,k,iv(NTAB),iy
      SAVE iv,iy
      DATA iv /NTAB*0/, iy /0/
!
      if (idum.le.0.or.iy.eq.0) then
        idum=max(-idum,1)
        do 11 j=NTAB+8,1,-1
          k=idum/IQ
          idum=IA*(idum-k*IQ)-IR*k
          if (idum.lt.0) idum=idum+IM
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM
      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)
      return
      END
c  ********************************************************
      FUNCTION rc(x,y)
      implicit double precision (a-h,o-z)
      PARAMETER (ERRTOL=0.0012,TINY=1.69e-38,SQRTNY=1.3e-19,BIG=3.E37,
     *TNBG=TINY*BIG,COMP1=2.236/SQRTNY,COMP2=TNBG*TNBG/25.,THIRD=1./3.,
     *C1=.3,C2=1./7.,C3=.375,C4=9./22.)
!
      if(x.lt.0..or.y.eq.0..or.(x+abs(y)).lt.TINY.or.(x+
     *abs(y)).gt.BIG.or.(y.lt.-COMP1.and.x.gt.0..and.x.lt.COMP2))stop 
     *'invalid arguments in rc'
      if(y.gt.0.)then
        xt=x
        yt=y
        w=1.
      else
        xt=x-y
        yt=-y
        w=sqrt(x)/sqrt(xt)
      endif
1     continue
        alamb=2.*sqrt(xt)*sqrt(yt)+yt
        xt=.25*(xt+alamb)
        yt=.25*(yt+alamb)
        ave=THIRD*(xt+yt+yt)
        s=(yt-ave)/ave
      if(abs(s).gt.ERRTOL)goto 1
      rc=w*(1.+s*s*(C1+s*(C2+s*(C3+s*C4))))/sqrt(ave)
      return
      END
c  ********************************************************
      FUNCTION rd(x,y,z)
      implicit double precision (a-h,o-z)
      PARAMETER (ERRTOL=0.0015,TINY=1.e-25,BIG=4.5E21,C1=3./14.,
     &  C2=1./6.,C3=9./22.,C4=3./26.,C5=.25*C3,C6=1.5*C4)
!
      if(min(x,y).lt.0..or.min(x+y,z).lt.TINY.or.max(x,y,
     *z).gt.BIG)stop 'invalid arguments in rd'
      xt=x
      yt=y
      zt=z
      sum=0.
      fac=1.
1     continue
        sqrtx=sqrt(xt)
        sqrty=sqrt(yt)
        sqrtz=sqrt(zt)
        alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz
        sum=sum+fac/(sqrtz*(zt+alamb))
        fac=.25*fac
        xt=.25*(xt+alamb)
        yt=.25*(yt+alamb)
        zt=.25*(zt+alamb)
        ave=.2*(xt+yt+3.*zt)
        delx=(ave-xt)/ave
        dely=(ave-yt)/ave
        delz=(ave-zt)/ave
      if(max(abs(delx),abs(dely),abs(delz)).gt.ERRTOL)goto 1
      ea=delx*dely
      eb=delz*delz
      ec=ea-eb
      ed=ea-6.*eb
      ee=ed+ec+ec
      rd=3.*sum+fac*(1.+ed*(-C1+C5*ed-C6*delz*ee)+delz*(C2*ee+delz*(-C3*
     *ec+delz*C4*ea)))/(ave*sqrt(ave))
      return
      END
c  ********************************************************
      FUNCTION rf(x,y,z)
      implicit double precision (a-h,o-z)
      PARAMETER (ERRTOL=0.0025,TINY=1.5e-38,BIG=3.E37,THIRD=1./3.,
     *C1=1./24.,C2=.1,C3=3./44.,C4=1./14.)
!
      if(min(x,y,z).lt.0..or.min(x+y,x+z,y+z).lt.TINY.or.max(x,y,
     *z).gt.BIG)stop 'invalid arguments in rf'
      xt=x
      yt=y
      zt=z
1     continue
        sqrtx=sqrt(xt)
        sqrty=sqrt(yt)
        sqrtz=sqrt(zt)
        alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz
        xt=.25*(xt+alamb)
        yt=.25*(yt+alamb)
        zt=.25*(zt+alamb)
        ave=THIRD*(xt+yt+zt)
        delx=(ave-xt)/ave
        dely=(ave-yt)/ave
        delz=(ave-zt)/ave
      if(max(abs(delx),abs(dely),abs(delz)).gt.ERRTOL)goto 1
      e2=delx*dely-delz**2
      e3=delx*dely*delz
      rf=(1.+(C1*e2-C2-C3*e3)*e2+C4*e3)/sqrt(ave)
      return
      END
c  ********************************************************
      FUNCTION rj(x,y,z,p)
      implicit double precision (a-h,o-z)
      PARAMETER (ERRTOL=0.0015,TINY=2.5e-13,BIG=9.E11,C1=3./14.,
     &  C2=1./3.,C3=3./22.,C4=3./26.,C5=.75*C3,C6=1.5*C4,C7=.5*C2,
     &  C8=C3+C3)
CU    USES rc,rf
!
      if(min(x,y,z).lt.0..or.min(x+y,x+z,y+z,abs(p)).lt.TINY.or.max(x,y,
     *z,abs(p)).gt.BIG) then
        stop 'invalid arguments in rj'
      end if
!
      sum=0.
      fac=1.
      if(p.gt.0.)then
        xt=x
        yt=y
        zt=z
        pt=p
      else
        xt=min(x,y,z)
        zt=max(x,y,z)
        yt=x+y+z-xt-zt
        a=1./(yt-p)
        b=a*(zt-yt)*(yt-xt)
        pt=yt+b
        rho=xt*zt/yt
        tau=p*pt/yt
        rcx=rc(rho,tau)
      endif
1     continue
        sqrtx=sqrt(xt)
        sqrty=sqrt(yt)
        sqrtz=sqrt(zt)
        alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz
        alpha=(pt*(sqrtx+sqrty+sqrtz)+sqrtx*sqrty*sqrtz)**2
        beta=pt*(pt+alamb)**2
        sum=sum+fac*rc(alpha,beta)
        fac=.25*fac
        xt=.25*(xt+alamb)
        yt=.25*(yt+alamb)
        zt=.25*(zt+alamb)
        pt=.25*(pt+alamb)
        ave=.2*(xt+yt+zt+pt+pt)
        delx=(ave-xt)/ave
        dely=(ave-yt)/ave
        delz=(ave-zt)/ave
        delp=(ave-pt)/ave
      if(max(abs(delx),abs(dely),abs(delz),abs(delp)).gt.ERRTOL)goto 1
      ea=delx*(dely+delz)+dely*delz
      eb=delx*dely*delz
      ec=delp**2
      ed=ea-3.*ec
      ee=eb+2.*delp*(ea-ec)
      rj=3.*sum+fac*(1.+ed*(-C1+C5*ed-C6*ee)+eb*(C7+delp*(-C8+delp*C4))+
     *delp*ea*(C2-delp*C3)-C2*delp*ec)/(ave*sqrt(ave))
      if (p.le.0.) rj=a*(b*rj+3.*(rcx-rf(xt,yt,zt)))
      return
      END
! *********************************************************
      FUNCTION rtsafe(funcd,x1,x2,xacc)
      implicit double precision (a-h,o-z)
      INTEGER MAXIT
      EXTERNAL funcd
      PARAMETER (MAXIT=100)
      INTEGER j
!
      call funcd(x1,fl,df)
      call funcd(x2,fh,df)
      if((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.))stop
     *'root must be bracketed in rtsafe'
      if(fl.eq.0.)then
        rtsafe=x1
        return
      else if(fh.eq.0.)then
        rtsafe=x2
        return
      else if(fl.lt.0.)then
        xl=x1
        xh=x2
      else
        xh=x1
        xl=x2
      endif
      rtsafe=.5*(x1+x2)
      dxold=abs(x2-x1)
      dx=dxold
      call funcd(rtsafe,f,df)
      do 11 j=1,MAXIT
        if(((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f).ge.0..or. abs(2.*
     *f).gt.abs(dxold*df) ) then
          dxold=dx
          dx=0.5*(xh-xl)
          rtsafe=xl+dx
          if(xl.eq.rtsafe)return
        else
          dxold=dx
          dx=f/df
          temp=rtsafe
          rtsafe=rtsafe-dx
          if(temp.eq.rtsafe)return
        endif
        if(abs(dx).lt.xacc) return
        call funcd(rtsafe,f,df)
        if(f.lt.0.) then
          xl=rtsafe
        else
          xh=rtsafe
        endif
11    continue
      stop 'rtsafe exceeding maximum iterations'
      return
      END
! *********************************************************
      SUBROUTINE trapzd(func,a,b,s,n)
      implicit double precision (a-h,o-z)
      INTEGER n
      EXTERNAL func
      INTEGER it,j
!
      if (n.eq.1) then
        s=0.5*(b-a)*(func(a)+func(b))
      else
        it=2**(n-2)
        tnm=it
        del=(b-a)/tnm
        x=a+0.5*del
        sum=0.
        do 11 j=1,it
          sum=sum+func(x)
          x=x+del
11      continue
        s=0.5*(s+(b-a)*sum/tnm)
      endif
      return
      END
!  *********************************************************
      subroutine TRILINEAR(u,v,w,jx,jy,jz,fg,mx,my,mz,f)
!
!  Do tri-linear interpolation 
!
      implicit double precision (a-h,o-z)
      dimension fg(mx,my,mz)
!
      f = (1.-u)*(1.-v)*(1.-w) * fg(jx,jy,jz)
     &    +  (1.-u)*(1.-v)*w * fg(jx,jy,jz+1)
     &    +  (1.-u)*v*w * fg(jx,jy+1,jz+1)
     &    +  (1.-u)*v*(1.-w) * fg(jx,jy+1,jz)
     &    +   u*(1.-v)*(1.-w) * fg(jx+1,jy,jz)
     &    +   u*(1.-v)*w * fg(jx+1,jy,jz+1)
     &    +   u*v*w * fg(jx+1,jy+1,jz+1)
     &    +   u*v*(1.-w) * fg(jx+1,jy+1,jz)
!
      return
      end
!  *********************************************************
      subroutine TRIPOLY(kx,ky,kz,xg,yg,zg,fg,mp,np,pp,kk,
     &  x,y,z,h,f,df)
!
!  Polynomial interpolation on 3-dimensional grid
!
!  cf. Num. Rec., 2nd ed, p.103,118; uses POLINT
!
!  kx,ky,kz: starting location in 3D grid arrays
!  xg(mp),yg(np),zg(pp): full coordinate arrays
!  fg(mp,np,pp): full field array
!  kk: polynomial order {1,2,3}
!  x,y,z: observation point
!  h: curvature
!  f: field value at observation point
!  df: error estimate for returned field
!
      implicit double precision (a-h,o-z)
      integer pp
      dimension fg(mp,np,pp),xg(mp),yg(np),zg(pp),zgh(pp)
      dimension f1tmp(4),f2tmp(4),f3tmp(4)
!
      do 300 k=1,kk+1       ! loop over z grid points
      do 200 j=1,kk+1       ! loop over y grid points
!
! Load temporary array f1tmp with x grid fields for fixed y and z        
      do 100 i=1,kk+1               
      f1tmp(i) = fg(kx+i-1,ky+j-1,kz+k-1)
100   continue
!
! f2tmp = field interpolated over x as a function of y for fixed z
      call POLINT(xg(kx),f1tmp,kk+1,x,f2tmp(j),df)
!
200   continue
!
! f3tmp = field interpolated over y as a function of z
      call POLINT(yg(ky),f2tmp,kk+1,y,f3tmp(k),df)
!
300   continue
!
! In cylidrical coordinate system we need to rescale the step sizes
!  for different radial positions
      do i=1,pp
      zgh(i) = zg(i) * (1. + h*x)
      end do
!
! f = field interpolated over z
      call POLINT(zgh(kz),f3tmp,kk+1,z,f,df)
!
      return
      end
