cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE IN(SCAL, X, NB, B, NCALC)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    calculates the i_n(z)=sqrt(pi/2/z)*i_(n+1/2)(z).
c
c  on input:
c    scal: the scaling factor to avoid underflow.
c    x: the parameter for i_n(x)
c    nb: number of terms for the subindex n, n=0:nb.
c
c  on output :
c    b: contains the i_n(x), n=0:nb.
c    ncalc: an error index. this is from the subroutine
c           for the calculation of bessel function.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NB,NCALC
      REAL *8 X,B(0:NB),SCAL
c
c-----local variables.
c
      INTEGER *4 I,IZE
      REAL *8 CONST,ENSIG,ENMTEN,HALFPI
      REAL *8 XSCAL,TERM1,TERM2,ALPHA
      DATA ENSIG, ENMTEN/1.0D-4, 1.0D-300/
c
c-----function called
c
      REAL *8 DABS
c
      IF (X.LT. 0.0D0) THEN
        PRINT *,' ERROR IN INPUT FOR IN.F'
      ENDIF
c
c-----if x.le. 10e-4, then use the 2 terms taylor expansion.
c
      IF (X.LE.ENSIG) THEN
        XSCAL=X/SCAL
        TERM1=1.0D0
        TERM2=0.5D0*X*X
        B(0)=TERM1*(1.0D0+TERM2/3.0D0)
        DO 1 I=1, NB
          TERM1=TERM1*XSCAL/DBLE(2*I+1)
          IF (TERM1.LE. ENMTEN) THEN
            TERM1=0.0D0
          ENDIF
          B(I)=TERM1*(1.0D0+TERM2/DBLE(2*I+3) )
1       CONTINUE
        NCALC=NB+1
c
c-------usually, scal should be less than one, however, in the
c         calculation of h_n, scal will be greater than one,
c         in this case, we should modify the program and replace every
c         small item with zero.
c
      ELSEIF (X.GT.1.0D+2) THEN
        DO I=0,NB
          B(I)=0.0D0
        ENDDO
        NCALC=NB+1
      ELSE
c
c-------otherwise, we will call besselj() and then scale the
c         jn by the scaling factor.
c
        HALFPI=DATAN(1.0D0)*2.0D0
        CONST=DSQRT(HALFPI/X)
c
c-------calculate i_(n+1/2).
c
        ALPHA=0.5D0
        IZE=1
        CALL RIBESL(X, ALPHA, NB+1, IZE, B, NCALC)
c
c-------now scale to i_n*scal**(-n).
c
        DO 2 I=0,NB
          B(I)=B(I)*CONST
          CONST=CONST/SCAL
          IF (DABS(B(I)) .LE. ENMTEN ) THEN
            CONST=0.0D0
          ENDIF
2       CONTINUE
      ENDIF
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE KN(SCAL,X,NB,BY,NCALC)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c purpose:
c  this routine calculates scaled sphecial bessel functions y_n (x)
c  for non-negative argument x, and n=0,...,nb.
c
c on input:
c   scal: the scaling factor.
c   x: working precision non-negative real argument for which
c      y's are to be calculated.
c   nb: integer number of functions to be calculated, nb .gt. 0.
c       the first function calculated is of order 0, and the
c       last is of order nb.
c   by: working precision output vector of length nb+1.  if the
c       routine terminates normally (ncalc=nb+1), the vector by
c       contains the functions y(0,x), ... , y(nb,x),
c       if (0 .lt. ncalc .lt. nb), by(i) contains correct function
c       values for i .le. ncalc, and contains the ratios
c       y(i-1,x)/y(i-2,x) for the rest of the array.
c   ncalc: integer output variable indicating possible errors.
c       before using the vector by, the user should check that
c       ncalc=nb, i.e., all orders have been calculated to
c       the desired accuracy.  see error returns below.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c explanation of machine-dependent constants
c
c   beta   = radix for the floating-point system
c   p      = number of significant base-beta digits in the
c            significand of a floating-point number
c   minexp = smallest representable power of beta
c   maxexp = smallest power of beta that overflows
c   eps    = beta ** (-p)
c   del    = machine number below which sin(x)/x = 1; approximately
c            sqrt(eps).
c   xmin   = smallest acceptable argument for rbesy; approximately
c            max(2*beta**minexp,2/xinf), rounded up
c   xinf   = largest positive machine number; approximately
c            beta**maxexp
c   thresh = lower bound for use of the asymptotic form; approximately
c            aint(-log10(eps/2.0))+1.0
c   xlarge = upper bound on x; approximately 1/del, because the sine
c            and cosine functions have lost about half of their
c            precision at that point.
c
c
c     approximate values for some important machines are:
c
c                        beta    p     minexp      maxexp      eps
c
c  cray-1        (s.p.)    2    48     -8193        8191    3.55e-15
c  cyber 180/185
c    under nos   (s.p.)    2    48      -975        1070    3.55e-15
c  ieee (ibm/xt,
c    sun, etc.)  (s.p.)    2    24      -126         128    5.96e-8
c  ieee (ibm/xt,
c    sun, etc.)  (d.p.)    2    53     -1022        1024    1.11d-16
c  ibm 3033      (d.p.)   16    14       -65          63    1.39d-17
c  vax           (s.p.)    2    24      -128         127    5.96e-8
c  vax d-format  (d.p.)    2    56      -128         127    1.39d-17
c  vax g-format  (d.p.)    2    53     -1024        1023    1.11d-16
c
c
c                         del      xmin      xinf     thresh  xlarge
c
c cray-1        (s.p.)  5.0e-8  3.67e-2466 5.45e+2465  15.0e0  2.0e7
c cyber 180/855
c   under nos   (s.p.)  5.0e-8  6.28e-294  1.26e+322   15.0e0  2.0e7
c ieee (ibm/xt,
c   sun, etc.)  (s.p.)  1.0e-4  2.36e-38   3.40e+38     8.0e0  1.0e4
c ieee (ibm/xt,
c   sun, etc.)  (d.p.)  1.0d-8  4.46d-308  1.79d+308   16.0d0  1.0d8
c ibm 3033      (d.p.)  1.0d-8  2.77d-76   7.23d+75    17.0d0  1.0d8
c vax           (s.p.)  1.0e-4  1.18e-38   1.70e+38     8.0e0  1.0e4
c vax d-format  (d.p.)  1.0d-9  1.18d-38   1.70d+38    17.0d0  1.0d9
c vax g-format  (d.p.)  1.0d-8  2.23d-308  8.98d+307   16.0d0  1.0d8
c
c*******************************************************************
c
c error returns
c
c  in case of an error, ncalc .ne. nb, and not all y's are
c  calculated to the desired accuracy.
c
c  ncalc .lt. -1:  an argument is out of range. for example,
c       nb .le. 0, ize is not 1 or 2, or ize=1 and abs(x) .ge.
c       xmax.  in this case, by(1) = 0.0, the remainder of the
c       by-vector is not calculated, and ncalc is set to
c       min0(nb,0)-2  so that ncalc .ne. nb.
c  ncalc = -1:  y(alpha,x) .ge. xinf.  the requested function
c       values are set to 0.0.
c  1 .lt. ncalc .lt. nb: not all requested function values could
c       be calculated accurately.  by(i) contains correct function
c       values for i .le. ncalc, and and the remaining nb-ncalc
c       array elements contain 0.0.
c
c
c intrinsic functions required are:
c
c     dble, exp, int, max, min, real, sqrt
c
c
c acknowledgement
c
c  this program draws heavily on temme's algol program for y(a,x)
c  and y(a+1,x) and on campbell's programs for y_nu(x).  temme's
c  scheme is used for  x < thresh, and campbell's scheme is used
c  in the asymptotic region.  segments of code from both sources
c  have been translated into fortran 77, merged, and heavily modified.
c  modifications include parameterization of machine dependencies,
c  use of a new approximation for ln(gamma(x)), and built-in
c  protection against over/underflow.
c
c references: "bessel functions j_nu(x) and y_nu(x) of real
c              order and real argument," campbell, j. b.,
c              comp. phy. comm. 18, 1979, pp. 133-142.
c
c             "on the numerical evaluation of the ordinary
c              bessel function of the second kind," temme,
c              n. m., j. comput. phys. 21, 1976, pp. 343-350.
c
c  latest modification: march 19, 1990
c
c  modified by: w. j. cody
c               applied mathematics division
c               argonne national laboratory
c               argonne, il  60439
c
c----------------------------------------------------------------------
      IMPLICIT REAL *8 (A-H,O-Z)
      INTEGER *4 I,NB,NCALC
      REAL *8 EX, BY(0:NB), P, SCAL, XINF, XMIN, XLARGE
      REAL *8 U1,U2
c----------------------------------------------------------------------
c  machine-dependent constants
c----------------------------------------------------------------------
      DATA XMIN,XINF/4.46D-308,1.79D308/
      DATA XLARGE/1.0D8/
c----------------------------------------------------------------------
      HALFPI=DATAN(1.0D0)*2.0D0
      EX = X
      IF ((NB.GE.0) .AND. (X.GE.XMIN) .AND. (EX.LT.XLARGE) )  THEN
c
c----------------------------------------------------------------------
c  now have first one or two y's
c----------------------------------------------------------------------
         P=DEXP(-EX)/EX*HALFPI
         BY(0) = P
         BY(1) = P*SCAL*(1.0D0+1.0D0/EX)
         NCALC = 1
         U1=SCAL/EX
         U2=SCAL*SCAL
         DO 400 I = 2, NB
           IF (DABS(BY(I-1))*U1 .GE. XINF/DBLE(2*I-1) ) GO TO 450
           BY(I) =  DBLE(2*I-1)*U1*BY(I-1)+U2*BY(I-2)
           NCALC = NCALC + 1
  400   CONTINUE
  450   DO 460 I = NCALC+1, NB
                BY(I) = ZERO
  460   CONTINUE
      ELSE
        BY(0)=0.0D0
        NCALC=MIN(NB+1,0)-1
      ENDIF
c
      RETURN
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE RIBESL(X,ALPHA,NB,IZE,B,NCALC)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c    this routine calculates bessel functions i sub(n+alpha) (x)
c    for non-negative argument x, and non-negative order n+alpha,
c    with or without exponential scaling.
c
c  on input:
c    x: working precision non-negative real argument for which
c       i's or exponentially scaled i's (i*exp(-x))
c       are to be calculated.  if i's are to be calculated,
c       x must be less than exparg (see below).
c    alpha: working precision fractional part of order for which
c           i's or exponentially scaled i's (i*exp(-x)) are
c           to be calculated.  0 .le. alpha .lt. 1.0.
c    nb: integer number of functions to be calculated, nb .gt. 0.
c        the first function calculated is of order alpha, and the
c        last is of order (nb - 1 + alpha).
c    ize: integer type.  ize = 1 if unscaled i's are to calculated,
c         and 2 if exponentially scaled i's are to be calculated.
c    b: working precision output vector of length nb.  if the routine
c       terminates normally (ncalc=nb), the vector b contains the
c       functions i(alpha,x) through i(nb-1+alpha,x), or the
c       corresponding exponentially scaled functions.
c    ncalc: integer output variable indicating possible errors.
c           before using the vector b, the user should check that
c           ncalc=nb, i.e., all orders have been calculated to
c           the desired accuracy.  see error returns below.
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c explanation of machine-dependent constants
c
c   beta   = radix for the floating-point system
c   minexp = smallest representable power of beta
c   maxexp = smallest power of beta that overflows
c   it     = number of bits in the mantissa of a working precision
c            variable
c   nsig   = decimal significance desired.  should be set to
c            int(log10(2)*it+1).  setting nsig lower will result
c            in decreased accuracy while setting nsig higher will
c            increase cpu time without increasing accuracy.  the
c            truncation error is limited to a relative error of
c            t=.5*10**(-nsig).
c   enten  = 10.0 ** k, where k is the largest integer such that
c            enten is machine-representable in working precision
c   ensig  = 10.0 ** nsig
c   rtnsig = 10.0 ** (-k) for the smallest integer k such that
c            k .ge. nsig/4
c   enmten = smallest abs(x) such that x/4 does not underflow
c   xlarge = upper limit on the magnitude of x when ize=2.  bear
c            in mind that if abs(x)=n, then at least n iterations
c            of the backward recursion will be executed.  the value
c            of 10.0 ** 4 is used on every machine.
c   exparg = largest working precision argument that the library
c            exp routine can handle and upper limit on the
c            magnitude of x when ize=1; approximately
c            log(beta**maxexp)
c
c
c     approximate values for some important machines are:
c
c                        beta       minexp      maxexp       it
c
c  cray-1        (s.p.)    2        -8193        8191        48
c  cyber 180/855
c    under nos   (s.p.)    2         -975        1070        48
c  ieee (ibm/xt,
c    sun, etc.)  (s.p.)    2         -126         128        24
c  ieee (ibm/xt,
c    sun, etc.)  (d.p.)    2        -1022        1024        53
c  ibm 3033      (d.p.)   16          -65          63        14
c  vax           (s.p.)    2         -128         127        24
c  vax d-format  (d.p.)    2         -128         127        56
c  vax g-format  (d.p.)    2        -1024        1023        53
c
c
c                        nsig       enten       ensig      rtnsig
c
c cray-1        (s.p.)    15       1.0e+2465   1.0e+15     1.0e-4
c cyber 180/855
c   under nos   (s.p.)    15       1.0e+322    1.0e+15     1.0e-4
c ieee (ibm/xt,
c   sun, etc.)  (s.p.)     8       1.0e+38     1.0e+8      1.0e-2
c ieee (ibm/xt,
c   sun, etc.)  (d.p.)    16       1.0d+308    1.0d+16     1.0d-4
c ibm 3033      (d.p.)     5       1.0d+75     1.0d+5      1.0d-2
c vax           (s.p.)     8       1.0e+38     1.0e+8      1.0e-2
c vax d-format  (d.p.)    17       1.0d+38     1.0d+17     1.0d-5
c vax g-format  (d.p.)    16       1.0d+307    1.0d+16     1.0d-4
c
c
c                         enmten      xlarge   exparg
c
c cray-1        (s.p.)   1.84e-2466   1.0e+4    5677
c cyber 180/855
c   under nos   (s.p.)   1.25e-293    1.0e+4     741
c ieee (ibm/xt,
c   sun, etc.)  (s.p.)   4.70e-38     1.0e+4      88
c ieee (ibm/xt,
c   sun, etc.)  (d.p.)   8.90d-308    1.0d+4     709
c ibm 3033      (d.p.)   2.16d-78     1.0d+4     174
c vax           (s.p.)   1.17e-38     1.0e+4      88
c vax d-format  (d.p.)   1.17d-38     1.0d+4      88
c vax g-format  (d.p.)   2.22d-308    1.0d+4     709
c
c*******************************************************************
c*******************************************************************
c
c error returns
c
c  in case of an error,  ncalc .ne. nb, and not all i's are
c  calculated to the desired accuracy.
c
c  ncalc .lt. 0:  an argument is out of range. for example,
c     nb .le. 0, ize is not 1 or 2, or ize=1 and abs(x) .ge. exparg.
c     in this case, the b-vector is not calculated, and ncalc is
c     set to min0(nb,0)-1 so that ncalc .ne. nb.
c
c  nb .gt. ncalc .gt. 0: not all requested function values could
c     be calculated accurately.  this usually occurs because nb is
c     much larger than abs(x).  in this case, b(n) is calculated
c     to the desired accuracy for n .le. ncalc, but precision
c     is lost for ncalc .lt. n .le. nb.  if b(n) does not vanish
c     for n .gt. ncalc (because it is too small to be represented),
c     and b(n)/b(ncalc) = 10**(-k), then only the first nsig-k
c     significant figures of b(n) can be trusted.
c
c
c intrinsic functions required are:
c
c     dble, exp, dgamma, gamma, int, max, min, real, sqrt
c
c
c acknowledgement
c
c  this program is based on a program written by david j.
c  sookne (2) that computes values of the bessel functions j or
c  i of real argument and integer order.  modifications include
c  the restriction of the computation to the i bessel function
c  of non-negative real argument, the extension of the computation
c  to arbitrary positive order, the inclusion of optional
c  exponential scaling, and the elimination of most underflow.
c  an earlier version was published in (3).
c
c references: "a note on backward recurrence algorithms," olver,
c              f. w. j., and sookne, d. j., math. comp. 26, 1972,
c              pp 941-947.
c
c             "bessel functions of real argument and integer order,"
c              sookne, d. j., nbs jour. of res. b. 77b, 1973, pp
c              125-132.
c
c             "algorithm 597, sequence of modified bessel functions
c              of the first kind," cody, w. j., trans. math. soft.,
c              1983, pp. 242-245.
c
c  latest modification: may 30, 1989
c
c  modified by: w. j. cody and l. stoltz
c               applied mathematics division
c               argonne national laboratory
c               argonne, il  60439
c
c-------------------------------------------------------------------
      IMPLICIT REAL *8 (A-H,O-Z)
      INTEGER *4 IZE,K,L,MAGX,N,NB,NBMX,NCALC,NEND,NSIG,NSTART
      DOUBLE PRECISION DGAMMA,
     1 ALPHA,B,CONST,CONV,EM,EMPAL,EMP2AL,EN,ENMTEN,ENSIG,
     2 ENTEN,EXPARG,FUNC,HALF,HALFX,ONE,P,PLAST,POLD,PSAVE,PSAVEL,
     3 RTNSIG,SUM,TEMPA,TEMPB,TEMPC,TEST,TOVER,TWO,X,XLARGE,ZERO
      DIMENSION B(NB)
c-------------------------------------------------------------------
c  mathematical constants
c-------------------------------------------------------------------
      DATA ONE,TWO,ZERO,HALF,CONST/1.0D0,2.0D0,0.0D0,0.5D0,1.585D0/
c-------------------------------------------------------------------
c  machine-dependent parameters
c-------------------------------------------------------------------
      DATA NSIG,XLARGE,EXPARG /16,1.0D4,709.0D0/
      DATA ENTEN,ENSIG,RTNSIG/1.0D308,1.0D16,1.0D-4/
      DATA ENMTEN/8.9D-308/
c-------------------------------------------------------------------
c  statement functions for conversion
c-------------------------------------------------------------------
      CONV(N) = DBLE(N)
      FUNC(X) = DGAMMA(X)
c-------------------------------------------------------------------
c check for x, nb, or ize out of range.
c-------------------------------------------------------------------
      IF ((NB.GT.0) .AND. (X .GE. ZERO) .AND.
     1    (ALPHA .GE. ZERO) .AND. (ALPHA .LT. ONE) .AND.
     2    (((IZE .EQ. 1) .AND. (X .LE. EXPARG)) .OR.
     3     ((IZE .EQ. 2) .AND. (X .LE. XLARGE)))) THEN
c-------------------------------------------------------------------
c use 2-term ascending series for small x
c-------------------------------------------------------------------
            NCALC = NB
            MAGX = INT(X)
            IF (X .GE. RTNSIG) THEN
c-------------------------------------------------------------------
c initialize the forward sweep, the p-sequence of olver
c-------------------------------------------------------------------
                  NBMX = NB-MAGX
                  N = MAGX+1
                  EN = CONV(N+N) + (ALPHA+ALPHA)
                  PLAST = ONE
                  P = EN / X
c-------------------------------------------------------------------
c calculate general significance test
c-------------------------------------------------------------------
                  TEST = ENSIG + ENSIG
                  IF (2*MAGX .GT. 5*NSIG) THEN
                        TEST = DSQRT(TEST*P)
                     ELSE
                        TEST = TEST / CONST**MAGX
                  END IF
                  IF (NBMX .GE. 3) THEN
c-------------------------------------------------------------------
c calculate p-sequence until n = nb-1.  check for possible overflow.
c-------------------------------------------------------------------
                     TOVER = ENTEN / ENSIG
                     NSTART = MAGX+2
                     NEND = NB - 1
                     DO 100 K = NSTART, NEND
                        N = K
                        EN = EN + TWO
                        POLD = PLAST
                        PLAST = P
                        P = EN * PLAST/X + POLD
                        IF (P .GT. TOVER) THEN
c-------------------------------------------------------------------
c to avoid overflow, divide p-sequence by tover.  calculate
c p-sequence until abs(p) .gt. 1.
c-------------------------------------------------------------------
                           TOVER = ENTEN
                           P = P / TOVER
                           PLAST = PLAST / TOVER
                           PSAVE = P
                           PSAVEL = PLAST
                           NSTART = N + 1
   60                      N = N + 1
                              EN = EN + TWO
                              POLD = PLAST
                              PLAST = P
                              P = EN * PLAST/X + POLD
                           IF (P .LE. ONE) GO TO 60
                           TEMPB = EN / X
c-------------------------------------------------------------------
c calculate backward test, and find ncalc, the highest n
c such that the test is passed.
c-------------------------------------------------------------------
                           TEST = POLD*PLAST / ENSIG
                           TEST = TEST*(HALF-HALF/(TEMPB*TEMPB))
                           P = PLAST * TOVER
                           N = N - 1
                           EN = EN - TWO
                           NEND = MIN0(NB,N)
                           DO 80 L = NSTART, NEND
                              NCALC = L
                              POLD = PSAVEL
                              PSAVEL = PSAVE
                              PSAVE = EN * PSAVEL/X + POLD
                              IF (PSAVE*PSAVEL .GT. TEST) GO TO 90
   80                      CONTINUE
                           NCALC = NEND + 1
   90                      NCALC = NCALC - 1
                           GO TO 120
                        END IF
  100                CONTINUE
                     N = NEND
                     EN = CONV(N+N) + (ALPHA+ALPHA)
c-------------------------------------------------------------------
c calculate special significance test for nbmx .gt. 2.
c-------------------------------------------------------------------
                     TEST = MAX(TEST,DSQRT(PLAST*ENSIG)*DSQRT(P+P))
                  END IF
c-------------------------------------------------------------------
c calculate p-sequence until significance test passed.
c-------------------------------------------------------------------
  110             N = N + 1
                     EN = EN + TWO
                     POLD = PLAST
                     PLAST = P
                     P = EN * PLAST/X + POLD
                  IF (P .LT. TEST) GO TO 110
c-------------------------------------------------------------------
c initialize the backward recursion and the normalization sum.
c-------------------------------------------------------------------
  120             N = N + 1
                  EN = EN + TWO
                  TEMPB = ZERO
                  TEMPA = ONE / P
                  EM = CONV(N) - ONE
                  EMPAL = EM + ALPHA
                  EMP2AL = (EM - ONE) + (ALPHA + ALPHA)
                  SUM = TEMPA * EMPAL * EMP2AL / EM
                  NEND = N - NB
                  IF (NEND .LT. 0) THEN
c-------------------------------------------------------------------
c n .lt. nb, so store b(n) and set higher orders to zero.
c-------------------------------------------------------------------
                        B(N) = TEMPA
                        NEND = -NEND
                        DO 130 L = 1, NEND
  130                      B(N+L) = ZERO
                     ELSE
                        IF (NEND .GT. 0) THEN
c-------------------------------------------------------------------
c recur backward via difference equation, calculating (but
c not storing) b(n), until n = nb.
c-------------------------------------------------------------------
                           DO 140 L = 1, NEND
                              N = N - 1
                              EN = EN - TWO
                              TEMPC = TEMPB
                              TEMPB = TEMPA
                              TEMPA = (EN*TEMPB) / X + TEMPC
                              EM = EM - ONE
                              EMP2AL = EMP2AL - ONE
                              IF (N .EQ. 1) GO TO 150
                              IF (N .EQ. 2) EMP2AL = ONE
                              EMPAL = EMPAL - ONE
                              SUM = (SUM + TEMPA*EMPAL) * EMP2AL / EM
  140                      CONTINUE
                        END IF
c-------------------------------------------------------------------
c store b(nb)
c-------------------------------------------------------------------
  150                   B(N) = TEMPA
                        IF (NB .LE. 1) THEN
                           SUM = (SUM + SUM) + TEMPA
                           GO TO 230
                        END IF
c-------------------------------------------------------------------
c calculate and store b(nb-1)
c-------------------------------------------------------------------
                        N = N - 1
                        EN = EN - TWO
                        B(N)  = (EN*TEMPA) / X + TEMPB
                        IF (N .EQ. 1) GO TO 220
                        EM = EM - ONE
                        EMP2AL = EMP2AL - ONE
                        IF (N .EQ. 2) EMP2AL = ONE
                        EMPAL = EMPAL - ONE
                        SUM = (SUM + B(N)*EMPAL) * EMP2AL / EM
                  END IF
                  NEND = N - 2
                  IF (NEND .GT. 0) THEN
c-------------------------------------------------------------------
c calculate via difference equation and store b(n), until n = 2.
c-------------------------------------------------------------------
                     DO 200 L = 1, NEND
                        N = N - 1
                        EN = EN - TWO
                        B(N) = (EN*B(N+1)) / X +B(N+2)
                        EM = EM - ONE
                        EMP2AL = EMP2AL - ONE
                        IF (N .EQ. 2) EMP2AL = ONE
                        EMPAL = EMPAL - ONE
                        SUM = (SUM + B(N)*EMPAL) * EMP2AL / EM
  200                CONTINUE
                  END IF
c-------------------------------------------------------------------
c calculate b(1)
c-------------------------------------------------------------------
                  B(1) = TWO*EMPAL*B(2) / X + B(3)
  220             SUM = (SUM + SUM) + B(1)
c-------------------------------------------------------------------
c normalize.  divide all b(n) by sum.
c-------------------------------------------------------------------
  230             IF (ALPHA .NE. ZERO)
     1               SUM = SUM * FUNC(ONE+ALPHA) * (X*HALF)**(-ALPHA)
                  IF (IZE .EQ. 1) SUM = SUM * DEXP(-X)
                  TEMPA = ENMTEN
                  IF (SUM .GT. ONE) TEMPA = TEMPA * SUM
                  DO 260 N = 1, NB
                     IF (B(N) .LT. TEMPA) B(N) = ZERO
                     B(N) = B(N) / SUM
  260             CONTINUE
                  RETURN
c-------------------------------------------------------------------
c two-term ascending series for small x.
c-------------------------------------------------------------------
               ELSE
                  TEMPA = ONE
                  EMPAL = ONE + ALPHA
                  HALFX = ZERO
                  IF (X .GT. ENMTEN) HALFX = HALF * X
                  IF (ALPHA .NE. ZERO) TEMPA = HALFX**ALPHA /FUNC(EMPAL)
                  IF (IZE .EQ. 2) TEMPA = TEMPA * DEXP(-X)
                  TEMPB = ZERO
                  IF ((X+ONE) .GT. ONE) TEMPB = HALFX * HALFX
                  B(1) = TEMPA + TEMPA*TEMPB / EMPAL
                  IF ((X .NE. ZERO) .AND. (B(1) .EQ. ZERO)) NCALC = 0
                  IF (NB .GT. 1) THEN
                     IF (X .EQ. ZERO) THEN
                           DO 310 N = 2, NB
                              B(N) = ZERO
  310                      CONTINUE
                        ELSE
c-------------------------------------------------------------------
c calculate higher-order functions.
c-------------------------------------------------------------------
                           TEMPC = HALFX
                           TOVER = (ENMTEN + ENMTEN) / X
                           IF (TEMPB .NE. ZERO) TOVER = ENMTEN / TEMPB
                           DO 340 N = 2, NB
                              TEMPA = TEMPA / EMPAL
                              EMPAL = EMPAL + ONE
                              TEMPA = TEMPA * TEMPC
                              IF (TEMPA .LE. TOVER*EMPAL) TEMPA = ZERO
                              B(N) = TEMPA + TEMPA*TEMPB / EMPAL
                              IF ((B(N) .EQ. ZERO) .AND. (NCALC .GT. N))
     1                             NCALC = N-1
  340                      CONTINUE
                     END IF
                  END IF
            END IF
         ELSE
            NCALC = MIN0(NB,0)-1
      END IF
      RETURN
c---------- last line of ribesl ----------
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      DOUBLE PRECISION FUNCTION DGAMMA(X)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c this routine calculates the gamma function for a real argument x.
c   computation is based on an algorithm outlined in reference 1.
c   the program uses rational functions that approximate the gamma
c   function to at least 20 significant decimal digits.  coefficients
c   for the approximation over the interval (1,2) are unpublished.
c   those for the approximation for x .ge. 12 are from reference 2.
c   the accuracy achieved depends on the arithmetic system, the
c   compiler, the intrinsic functions, and proper selection of the
c   machine-dependent constants.
c
c
c*******************************************************************
c*******************************************************************
c
c explanation of machine-dependent constants
c
c beta   - radix for the floating-point representation
c maxexp - the smallest positive power of beta that overflows
c xbig   - the largest argument for which gamma(x) is representable
c          in the machine, i.e., the solution to the equation
c                  gamma(xbig) = beta**maxexp
c xinf   - the largest machine representable floating-point number;
c          approximately beta**maxexp
c eps    - the smallest positive floating-point number such that
c          1.0+eps .gt. 1.0
c xminin - the smallest positive floating-point number such that
c          1/xminin is machine representable
c
c     approximate values for some important machines are:
c
c                            beta       maxexp        xbig
c
c cray-1         (s.p.)        2         8191        966.961
c cyber 180/855
c   under nos    (s.p.)        2         1070        177.803
c ieee (ibm/xt,
c   sun, etc.)   (s.p.)        2          128        35.040
c ieee (ibm/xt,
c   sun, etc.)   (d.p.)        2         1024        171.624
c ibm 3033       (d.p.)       16           63        57.574
c vax d-format   (d.p.)        2          127        34.844
c vax g-format   (d.p.)        2         1023        171.489
c
c                            xinf         eps        xminin
c
c cray-1         (s.p.)   5.45e+2465   7.11e-15    1.84e-2466
c cyber 180/855
c   under nos    (s.p.)   1.26e+322    3.55e-15    3.14e-294
c ieee (ibm/xt,
c   sun, etc.)   (s.p.)   3.40e+38     1.19e-7     1.18e-38
c ieee (ibm/xt,
c   sun, etc.)   (d.p.)   1.79d+308    2.22d-16    2.23d-308
c ibm 3033       (d.p.)   7.23d+75     2.22d-16    1.39d-76
c vax d-format   (d.p.)   1.70d+38     1.39d-17    5.88d-39
c vax g-format   (d.p.)   8.98d+307    1.11d-16    1.12d-308
c
c*******************************************************************
c*******************************************************************
c
c error returns
c
c  the program returns the value xinf for singularities or
c     when overflow would occur.  the computation is believed
c     to be free of underflow and overflow.
c
c
c  intrinsic functions required are:
c
c     int, dble, exp, log, real, sin
c
c
c references: "an overview of software development for special
c              functions", w. j. cody, lecture notes in mathematics,
c              506, numerical analysis dundee, 1975, g. a. watson
c              (ed.), springer verlag, berlin, 1976.
c
c              computer approximations, hart, et. al., wiley and
c              sons, new york, 1968.
c
c  latest modification: october 12, 1989
c
c  authors: w. j. cody and l. stoltz
c           applied mathematics division
c           argonne national laboratory
c           argonne, il 60439
c
c----------------------------------------------------------------------
      IMPLICIT REAL *8 (A-H,O-Z)
      INTEGER *4 I,N
      LOGICAL PARITY
      DOUBLE PRECISION
     1    C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE,
     2    TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
      DIMENSION C(7),P(8),Q(8)
c----------------------------------------------------------------------
c  mathematical constants
c----------------------------------------------------------------------
      DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/,
     1     SQRTPI/0.9189385332046727417803297D0/,
     2     PI/3.1415926535897932384626434D0/
c----------------------------------------------------------------------
c  machine dependent parameters
c----------------------------------------------------------------------
      DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/,
     1     XINF/1.79D308/
c----------------------------------------------------------------------
c  numerator and denominator coefficients for rational minimax
c     approximation over (1,2).
c----------------------------------------------------------------------
      DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1,
     1       -3.79804256470945635097577D+2,6.29331155312818442661052D+2,
     2       8.66966202790413211295064D+2,-3.14512729688483675254357D+4,
     3       -3.61444134186911729807069D+4,6.64561438202405440627855D+4/
      DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2,
     1      -1.01515636749021914166146D+3,-3.10777167157231109440444D+3,
     2        2.25381184209801510330112D+4,4.75584627752788110767815D+3,
     3      -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/
c----------------------------------------------------------------------
c  coefficients for minimax approximation over (12, inf).
c----------------------------------------------------------------------
      DATA C/-1.910444077728D-03,8.4171387781295D-04,
     1     -5.952379913043012D-04,7.93650793500350248D-04,
     2     -2.777777777777681622553D-03,8.333333333333333331554247D-02,
     3      5.7083835261D-03/
c----------------------------------------------------------------------
c  statement functions for conversion between integer and float
c----------------------------------------------------------------------
      CONV(I) = DBLE(I)
      PARITY = .FALSE.
      FACT = ONE
      N = 0
      Y = X
      IF (Y .LE. ZERO) THEN
c----------------------------------------------------------------------
c  argument is negative
c----------------------------------------------------------------------
            Y = -X
            Y1 = DINT(Y)
            RES = Y - Y1
            IF (RES .NE. ZERO) THEN
                  IF (Y1 .NE. DINT(Y1*HALF)*TWO) PARITY = .TRUE.
                  FACT = -PI / SIN(PI*RES)
                  Y = Y + ONE
               ELSE
                  RES = XINF
                  GO TO 900
            END IF
      END IF
c----------------------------------------------------------------------
c  argument is positive
c----------------------------------------------------------------------
      IF (Y .LT. EPS) THEN
c----------------------------------------------------------------------
c  argument .lt. eps
c----------------------------------------------------------------------
            IF (Y .GE. XMININ) THEN
                  RES = ONE / Y
               ELSE
                  RES = XINF
                  GO TO 900
            END IF
         ELSE IF (Y .LT. TWELVE) THEN
            Y1 = Y
            IF (Y .LT. ONE) THEN
c----------------------------------------------------------------------
c  0.0 .lt. argument .lt. 1.0
c----------------------------------------------------------------------
                  Z = Y
                  Y = Y + ONE
               ELSE
c----------------------------------------------------------------------
c  1.0 .lt. argument .lt. 12.0, reduce argument if necessary
c----------------------------------------------------------------------
                  N = INT(Y) - 1
                  Y = Y - CONV(N)
                  Z = Y - ONE
            END IF
c----------------------------------------------------------------------
c  evaluate approximation for 1.0 .lt. argument .lt. 2.0
c----------------------------------------------------------------------
            XNUM = ZERO
            XDEN = ONE
            DO 260 I = 1, 8
               XNUM = (XNUM + P(I)) * Z
               XDEN = XDEN * Z + Q(I)
  260       CONTINUE
            RES = XNUM / XDEN + ONE
            IF (Y1 .LT. Y) THEN
c----------------------------------------------------------------------
c  adjust result for case  0.0 .lt. argument .lt. 1.0
c----------------------------------------------------------------------
                  RES = RES / Y1
               ELSE IF (Y1 .GT. Y) THEN
c----------------------------------------------------------------------
c  adjust result for case  2.0 .lt. argument .lt. 12.0
c----------------------------------------------------------------------
                  DO 290 I = 1, N
                     RES = RES * Y
                     Y = Y + ONE
  290             CONTINUE
            END IF
         ELSE
c----------------------------------------------------------------------
c  evaluate for argument .ge. 12.0,
c----------------------------------------------------------------------
            IF (Y .LE. XBIG) THEN
                  YSQ = Y * Y
                  SUM = C(7)
                  DO 350 I = 1, 6
                     SUM = SUM / YSQ + C(I)
  350             CONTINUE
                  SUM = SUM/Y - Y + SQRTPI
                  SUM = SUM + (Y-HALF)*DLOG(Y)
                  RES = DEXP(SUM)
               ELSE
                  RES = XINF
                  GO TO 900
            END IF
      END IF
c----------------------------------------------------------------------
c  final adjustments and return
c----------------------------------------------------------------------
      IF (PARITY) RES = -RES
      IF (FACT .NE. ONE) RES = FACT / RES
  900 DGAMMA = RES
      RETURN
c ---------- last line of gamma ----------
      END
c
c**********************************************************************
      SUBROUTINE LGNDRGT1(SCAL,NMAX, X, Y)
c**********************************************************************
c
c-----this subroutine calculates the legendre function for x >1.
c     the scaled version. scal^n* p_n^m
c     scal should be less than 1.
c
c     input :
c       scal : the scaling factor.
c       nmax : the number of terms we want to compute.
c       x : the value at which we want to compute.
c
c     output :
c       y() : the scal^n* p_n^m.
c
c***********************************************************************
      IMPLICIT REAL *8 (A-H,O-Z)
      INTEGER *4 NMAX, M
      REAL *8  SCAL, X, Y(0:NMAX,0:NMAX), U, V, W
c
      V=SCAL*X
      W=SCAL*SCAL
      U=DSQRT(X*X-1)*SCAL
      Y(0,0)=1.0D0
      DO 10 M=0, NMAX
         IF (M.GT.0)  Y(M,M)=Y(M-1,M-1)*U*DBLE(2*M-1)
         IF (M.LT.NMAX)  Y(M+1,M)=DBLE(2*M+1)*V*Y(M,M)
         DO 20 N=M+2, NMAX
            Y(N,M)=((2.0D0*DBLE(N)-1.0D0)*V*Y(N-1,M)-DBLE(N+M-1)
     1         *W*Y(N-2,M)) / DBLE(N-M)
 20      CONTINUE
 10   CONTINUE
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE LGNDR(NMAX, X, Y)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    this subroutine computes the lengre polynomial expansion using
c    a recursive expansion.
c
c  on input:
c    nmax: the max number of terms in the expansion.
c    x: where we want to evaluate the expansion.
c
c  on output:
c    y: the function value at x.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NMAX
      REAL *8 X,Y(0:NMAX,0:NMAX)
c
c-----local variables.
c
      INTEGER *4 M,N
      REAL *8 U,DSQRT,DBLE
c
      U=-DSQRT(1.0D0-X*X)
      Y(0,0)=1.0D0
      DO 10 M=0, NMAX
        IF (M.GT.0)  Y(M,M)=Y(M-1,M-1)*U*DBLE(2*M-1)
        IF (M.LT.NMAX)  Y(M+1,M)=DBLE(2*M+1)*X*Y(M,M)
        DO 20 N=M+2, NMAX
          Y(N,M)=((2.0D0*DBLE(N)-1.0D0)*X*Y(N-1,M)-DBLE(N+M-1)
     1      *Y(N-2,M)) / DBLE(N-M)
 20     CONTINUE
 10   CONTINUE
c
      RETURN
      END
c
