C
C  ---------------------------------------
C  Compare scalar fnlib routines to vfnlib
C  ---------------------------------------
C
C  Ron Boisvert and Bonita Saunders
C  Computing and Applied Mathematics Laboratory
C  National Institute of Standards and Technology
C  Gaithersburg, MD 20899
C
C  boisvert@cam.nist.gov
C  saunders@cam.nist.gov
C
C  21 Feb 91
C
C----------------------------------------------------------------------
C
      INTEGER M
      PARAMETER ( M=2000 )
C
      REAL      XLO, XHI, X, F, FX, WORK
      INTEGER   I, IWORK, J
C
      DIMENSION XLO(3,8), XHI(3,8), X(M), F(M), FX(M), WORK(7*M),
     +          IWORK(M)
C
      SAVE XLO, XHI
C
      EXTERNAL  VI0, BESI0, VI1, BESI1, VJ0, BESJ0, VJ1, BESJ1,
     +          VK0, BESK0, VK1, BESK1, VY0, BESY0, VY1, BESY1
C
      DATA ((XLO(I,J),XHI(I,J),I=1,3),J=1,8) /
     +     0.0E0, 3.0E0 , -8.0E0, -3.0E0, 10.0E0, 50.0E0,
     +     0.0E0, 3.0E0 , -8.0E0, -3.0E0, 10.0E0, 50.0E0,
     +     0.0E0, 4.0E0 , -8.0E0, -4.0E0,  8.0E0, 50.0E0,
     +     0.0E0, 4.0E0 , -8.0E0, -4.0E0,  8.0E0, 50.0E0,
     +     0.0E0, 2.0E0 ,  2.0E0,  8.0E0, 10.0E0, 50.0E0,
     +     0.0E0, 2.0E0 ,  2.0E0,  8.0E0, 10.0E0, 50.0E0,
     +     0.0E0, 1.0E-3,  1.0E0,  4.0E0,  4.0E0, 50.0E0,
     +     0.0E0, 1.0E-3,  1.0E0,  4.0E0,  4.0E0, 50.0E0  /
C
C----------------------------------------------------------------------
C
C     ... TEST EACH USER-CALLABLE ROUTINE
C
      CALL TEST('I0',BESI0,VI0,XLO(1,1),XHI(1,1),M,X,F,FX,WORK,IWORK)
      CALL TEST('I1',BESI1,VI1,XLO(1,2),XHI(1,2),M,X,F,FX,WORK,IWORK)
      CALL TEST('J0',BESJ0,VJ0,XLO(1,3),XHI(1,3),M,X,F,FX,WORK,IWORK)
      CALL TEST('J1',BESJ1,VJ1,XLO(1,4),XHI(1,4),M,X,F,FX,WORK,IWORK)
      CALL TEST('K0',BESK0,VK0,XLO(1,5),XHI(1,5),M,X,F,FX,WORK,IWORK)
      CALL TEST('K1',BESK1,VK1,XLO(1,6),XHI(1,6),M,X,F,FX,WORK,IWORK)
      CALL TEST('Y0',BESY0,VY0,XLO(1,7),XHI(1,7),M,X,F,FX,WORK,IWORK)
      CALL TEST('Y1',BESY1,VY1,XLO(1,8),XHI(1,8),M,X,F,FX,WORK,IWORK)
C
      END

      SUBROUTINE TEST (NAME, SFUN, VSUB, XLOW, XHIGH, M, X, F, FX,
     +                 WORK, IWORK)
C
C----------------------------------------------------------------------
C
C  compare scalar function sfun to vectorized subroutine vsub for m
C  arguments distributed in various ways in the intervals
C  (xlo(i),xhi(i)), i=1,2,3.
C
C----------------------------------------------------------------------
C
C  ----------
C  PARAMETERS
C  ----------
C
      REAL    SFUN, XLOW, XHIGH, X, F, FX, WORK
      INTEGER IWORK, M
      CHARACTER*(*) NAME
C
      DIMENSION XLOW(3), XHIGH(3), X(M), F(M), FX(M), WORK(6*M),
     +          IWORK(M)
C
      EXTERNAL SFUN, VSUB
C
C  ---------------
C  LOCAL VARIABLES
C  ---------------
C
      LOGICAL DEBUG
      INTEGER NREP, NTRY
      PARAMETER ( DEBUG=.FALSE., NREP=15, NTRY=19 )
C
      INTEGER I, INFO, IREP, IOFF, IPC(NTRY), JPC(NTRY), K, KPC,
     +        MA, MB, MC, NFAIL, NREPV
      REAL    DIFF, DMAX, DX, EPS, R1MACH
      REAL SECOND, SU, SUMAX, SUMIN, SUAVE, T0, T1, TS, TSAVE, TV,
     +        TVAVE
C
      SAVE IPC, JPC
C
      DATA IPC /  0,  0,  0,  0,  0, 25, 25, 25, 25, 33, 50, 50, 50,
     +           75, 75,100,  1,  1, 98/
      DATA JPC /  0, 25, 50, 75,100,  0, 25, 50, 75, 33,  0, 25, 50,
     +            0, 25,  0,  1, 98,  1/
C
C----------------------------------------------------------------------
C
C  ntry = number of trials (each with different argument distribution)
C  nrep = number of times test repeated to get reliable timings
C
C  the distribution of arguments in trial k is
C
C     ipc(k)%            in (xlow(1),xhigh(1))
C     jpc(k)%            in (xlow(2),xhigh(2))
C     100-ipc(k)-jpc(k)% in (xlow(3),xhigh(3))
C
C  scalar results (fx) and vector results (f) are compared and accepted
C  if   abs((f - fx)/fx) < eps  for fx > 1  (relative difference)
C  or   abs(f - fx)      < eps  for fx < 1  (absolute difference)
C
C----------------------------------------------------------------------
C
      EPS = 5.0E0*R1MACH(4)

C     ... TRIGGER INITIALIZATIONS
C
      X(1) = 1.0E0
      CALL VSUB(1,X,F,WORK,IWORK,INFO)
      FX(1) = SFUN(X(1))
C
      PRINT *
      PRINT *
      PRINT *,' ----------'
      PRINT *,' TEST OF ',NAME
      PRINT *,' ----------'
      PRINT *
      PRINT *,' Vector Length = ',M
      PRINT *,' Repetitions   = ',NREP
      PRINT *,' Eps           = ',EPS
      PRINT *
      PRINT *,' Range A = (',XLOW(1),',',XHIGH(1),')'
      PRINT *,' Range B = (',XLOW(2),',',XHIGH(2),')'
      PRINT *,' Range C = (',XLOW(3),',',XHIGH(3),')'
      PRINT *
      PRINT *,' %A      = percent arguments in range A'
      PRINT *,' %B      = percent arguments in range B'
      PRINT *,' %C      = percent arguments in range C'
      PRINT *,' Scalar  = time for scalar code (FNLIB)'
      PRINT *,' Vector  = time for vector code (VFNLIB)'
      PRINT *,' SU      = speedup (Scalar/Vector)'
      PRINT *,' Nfail   = number of scalar-vector differences .gt. eps'
      PRINT *,' Dmax    = maximum scalar-vector difference'
      PRINT *
      PRINT *,' All times in seconds/argument.'
      PRINT *
      PRINT *,' Scalar-vector differences measured by'
      PRINT *,'     abs(s-v)/max(abs(s),1)'
      PRINT *,' where s,v are scalar,vector results.'
      PRINT *
      PRINT *,'  %A  %B  %C    Scalar    Vector   SU   Nfail   Dmax'
      PRINT *,' -------------------------------------------------------'
C
      SUMAX = 0.0
      SUMIN = 1000.0
      SUAVE = 0.0
      TSAVE = 0.0
      TVAVE = 0.0
C
C     ... PERFORM EACH TEST
C
      DO 100 K=1,NTRY
C
         MA = IPC(K)*M/100
         MB = JPC(K)*M/100
         KPC = 100 - IPC(K) - JPC(K)
         MC = M - MA - MB
C
C        ... SELECT X VECTOR
C
         IOFF = 0
C
         IF (MA .GT. 0) THEN
            DX = (XHIGH(1) - XLOW(1))/MA
            DO 10 I=1,MA
               X(I) = XLOW(1) + I*DX
   10       CONTINUE
            IOFF = IOFF + MA
         ENDIF
C
         IF (MB .GT. 0) THEN
            DX = (XHIGH(2) - XLOW(2))/MB
            DO 20 I=1,MB
               X(I+IOFF) = XLOW(2) + I*DX
   20       CONTINUE
            IOFF = IOFF + MB
         ENDIF
C
         IF (MC .GT. 0) THEN
            DX = (XHIGH(3) - XLOW(3))/MC
            DO 30 I=1,MC
               X(I+IOFF) = XLOW(3) + I*DX
   30       CONTINUE
         ENDIF
C
C        ... RUN VECTOR CODE
C
	 NREPV = 5*NREP
         T0 = SECOND()
         DO 40 IREP=1,NREPV
            CALL VSUB(M,X,F,WORK,IWORK,INFO)
   40    CONTINUE
         T1 = SECOND()
         TV = (T1 - T0)/NREPV
         IF (INFO .LT. 0) PRINT *,'Warning in ',NAME,' : Info = ',INFO
         IF (INFO .GT. 0) PRINT *,'Failure in ',NAME,' : Info = ',INFO,
     +                            '  Position = ',IWORK(1)
C
C        ... RUN SCALAR CODE
C
         T0 = SECOND()
         DO 60 IREP=1,NREP
            DO 50 I=1,M
               FX(I) = SFUN(X(I))
   50       CONTINUE
   60    CONTINUE
         T1 = SECOND()
         TS = (T1 - T0)/NREP
C
C        ... COMPUTE TIMINGS
C
         SU = TS/TV
         SUMIN = MIN(SU,SUMIN)
         SUMAX = MAX(SU,SUMAX)
         SUAVE = SUAVE + SU
         TSAVE = TSAVE + TS
         TVAVE = TVAVE + TV
C
C        ... EVALUATE RESULTS
C
         DMAX = 0.0E0
         NFAIL = 0
         DO 70 I=1,M
            DIFF = ABS(FX(I)-F(I))/MAX(ABS(FX(I)),1.0E0)
            DMAX = MAX(DIFF,DMAX)
            IF (DIFF .GT. EPS) THEN
               NFAIL = NFAIL + 1
               IF (DEBUG) PRINT *,'Failure ',I,X(I),F(I),FX(I),DIFF
            ENDIF
  70     CONTINUE
C
C         ... PRINT RESULTS
C
          PRINT 1000, IPC(K),JPC(K),KPC,TS/M,TV/M,SU,NFAIL,DMAX
C
  100 CONTINUE
C
      SUAVE = SUAVE/NTRY
      TSAVE = TSAVE/M/NTRY
      TVAVE = TVAVE/M/NTRY
C
      PRINT *
      PRINT 1002,'Ave Scalar   = ',TSAVE,' sec/argument'
      PRINT 1002,'Ave Vector   = ',TVAVE,' sec/argument'
      PRINT *
      PRINT 1001,'Max speed up = ',SUMAX
      PRINT 1001,'Min speed up = ',SUMIN
      PRINT *
      PRINT 1001,'Ave speed up = ',SUAVE
C
      RETURN
 1000 FORMAT(3(1X,I3),1P,2(3X,E7.1),0P,1X,F5.1,1X,I4,3X,1P,E7.1)
 1001 FORMAT(1X,A,F5.1,A)
 1002 FORMAT(1X,A,1P,E7.1,A)
      END
      function besi0 (x)
c april 1977 version.  w. fullerton, c3, los alamos scientific lab.
      dimension bi0cs(12)
c     external alog, besi0e, csevl, exp, inits, r1mach, sqrt
c
c series for bi0        on the interval  0.          to  9.00000d+00
c                                        with weighted error   2.46e-18
c                                         log weighted error  17.61
c                               significant figures required  17.90
c                                    decimal places required  18.15
c
      data bi0 cs( 1) /   -.0766054725 2839144951e0 /
      data bi0 cs( 2) /   1.9273379539 93808270e0 /
      data bi0 cs( 3) /    .2282644586 920301339e0 /
      data bi0 cs( 4) /    .0130489146 6707290428e0 /
      data bi0 cs( 5) /    .0004344270 9008164874e0 /
      data bi0 cs( 6) /    .0000094226 5768600193e0 /
      data bi0 cs( 7) /    .0000001434 0062895106e0 /
      data bi0 cs( 8) /    .0000000016 1384906966e0 /
      data bi0 cs( 9) /    .0000000000 1396650044e0 /
      data bi0 cs(10) /    .0000000000 0009579451e0 /
      data bi0 cs(11) /    .0000000000 0000053339e0 /
      data bi0 cs(12) /    .0000000000 0000000245e0 /
c
      data nti0, xsml, xmax / 0, 0., 0. /
c
      if (nti0.ne.0) go to 10
      nti0 = inits (bi0cs, 12, 0.1*r1mach(3))
      xsml = sqrt (4.0*r1mach(3))
      xmax = alog (r1mach(2))
c
 10   y = abs(x)
      if (y.gt.3.0) go to 20
c
      besi0 = 1.0
      if (y.gt.xsml) besi0 = 2.75 + csevl (y*y/4.5-1.0, bi0cs, nti0)
      return
c
 20   if (y.gt.xmax) call seteru (34hbesi0   abs(x) so big i0 overflows,
     1  34, 1, 2)
c
      besi0 = exp(y) * besi0e(x)
c
      return
      end
      function besi1 (x)
c oct 1983 version.  w. fullerton, c3, los alamos scientific lab.
      dimension bi1cs(11)
c     external alog, besi1e, csevl, exp, inits, r1mach, sqrt
c
c series for bi1        on the interval  0.          to  9.00000d+00
c                                        with weighted error   2.40e-17
c                                         log weighted error  16.62
c                               significant figures required  16.23
c                                    decimal places required  17.14
c
      data bi1 cs( 1) /   -.0019717132 61099859e0 /
      data bi1 cs( 2) /    .4073488766 7546481e0 /
      data bi1 cs( 3) /    .0348389942 99959456e0 /
      data bi1 cs( 4) /    .0015453945 56300123e0 /
      data bi1 cs( 5) /    .0000418885 21098377e0 /
      data bi1 cs( 6) /    .0000007649 02676483e0 /
      data bi1 cs( 7) /    .0000000100 42493924e0 /
      data bi1 cs( 8) /    .0000000000 99322077e0 /
      data bi1 cs( 9) /    .0000000000 00766380e0 /
      data bi1 cs(10) /    .0000000000 00004741e0 /
      data bi1 cs(11) /    .0000000000 00000024e0 /
c
c
      data nti1, xmin, xsml, xmax / 0, 3*0. /
c
      if (nti1.ne.0) go to 10
      nti1 = inits (bi1cs, 11, 0.1*r1mach(3))
      xmin = 2.0*r1mach(1)
      xsml = sqrt (8.0*r1mach(3))
      xmax = alog (r1mach(2))
c
 10   y = abs(x)
      if (y.gt.3.0) go to 20
c
      besi1 = 0.0
      if (y.eq.0.0) return
c
      if (y.le.xmin) call seteru (
     1  37hbesi1   abs(x) so small i1 underflows, 37, 1, 0)
      if (y.gt.xmin) besi1 = 0.5*x
      if (y.gt.xsml) besi1 = x * (.875 + csevl(y*y/4.5-1., bi1cs, nti1))
      return
c
 20   if (y.gt.xmax) call seteru (
     1  34hbesi1   abs(x) so big i1 overflows, 34, 2, 2)
c
      besi1 = exp(y) * besi1e(x)
c
      return
      end
      function besj0 (x)
c april 1977 version.  w. fullerton, c3, los alamos scientific lab.
      dimension bj0cs(13), bm0cs(21), bth0cs(24)
c     external cos, csevl, inits, r1mach, sqrt
c
c series for bj0        on the interval  0.          to  1.60000d+01
c                                        with weighted error   7.47e-18
c                                         log weighted error  17.13
c                               significant figures required  16.98
c                                    decimal places required  17.68
c
      data bj0 cs( 1) /    .1002541619 68939137e0 /
      data bj0 cs( 2) /   -.6652230077 64405132e0 /
      data bj0 cs( 3) /    .2489837034 98281314e0 /
      data bj0 cs( 4) /   -.0332527231 700357697e0 /
      data bj0 cs( 5) /    .0023114179 304694015e0 /
      data bj0 cs( 6) /   -.0000991127 741995080e0 /
      data bj0 cs( 7) /    .0000028916 708643998e0 /
      data bj0 cs( 8) /   -.0000000612 108586630e0 /
      data bj0 cs( 9) /    .0000000009 838650793e0 /
      data bj0 cs(10) /   -.0000000000 124235515e0 /
      data bj0 cs(11) /    .0000000000 001265433e0 /
      data bj0 cs(12) /   -.0000000000 000010619e0 /
      data bj0 cs(13) /    .0000000000 000000074e0 /
c
c series for bm0        on the interval  0.          to  6.25000d-02
c                                        with weighted error   4.98e-17
c                                         log weighted error  16.30
c                               significant figures required  14.97
c                                    decimal places required  16.96
c
      data bm0 cs( 1) /    .0928496163 7381644e0 /
      data bm0 cs( 2) /   -.0014298770 7403484e0 /
      data bm0 cs( 3) /    .0000283057 9271257e0 /
      data bm0 cs( 4) /   -.0000014330 0611424e0 /
      data bm0 cs( 5) /    .0000001202 8628046e0 /
      data bm0 cs( 6) /   -.0000000139 7113013e0 /
      data bm0 cs( 7) /    .0000000020 4076188e0 /
      data bm0 cs( 8) /   -.0000000003 5399669e0 /
      data bm0 cs( 9) /    .0000000000 7024759e0 /
      data bm0 cs(10) /   -.0000000000 1554107e0 /
      data bm0 cs(11) /    .0000000000 0376226e0 /
      data bm0 cs(12) /   -.0000000000 0098282e0 /
      data bm0 cs(13) /    .0000000000 0027408e0 /
      data bm0 cs(14) /   -.0000000000 0008091e0 /
      data bm0 cs(15) /    .0000000000 0002511e0 /
      data bm0 cs(16) /   -.0000000000 0000814e0 /
      data bm0 cs(17) /    .0000000000 0000275e0 /
      data bm0 cs(18) /   -.0000000000 0000096e0 /
      data bm0 cs(19) /    .0000000000 0000034e0 /
      data bm0 cs(20) /   -.0000000000 0000012e0 /
      data bm0 cs(21) /    .0000000000 0000004e0 /
c
c series for bth0       on the interval  0.          to  6.25000d-02
c                                        with weighted error   3.67e-17
c                                         log weighted error  16.44
c                               significant figures required  15.53
c                                    decimal places required  17.13
c
      data bth0cs( 1) /   -.2463916377 4300119e0 /
      data bth0cs( 2) /    .0017370983 07508963e0 /
      data bth0cs( 3) /   -.0000621836 33402968e0 /
      data bth0cs( 4) /    .0000043680 50165742e0 /
      data bth0cs( 5) /   -.0000004560 93019869e0 /
      data bth0cs( 6) /    .0000000621 97400101e0 /
      data bth0cs( 7) /   -.0000000103 00442889e0 /
      data bth0cs( 8) /    .0000000019 79526776e0 /
      data bth0cs( 9) /   -.0000000004 28198396e0 /
      data bth0cs(10) /    .0000000001 02035840e0 /
      data bth0cs(11) /   -.0000000000 26363898e0 /
      data bth0cs(12) /    .0000000000 07297935e0 /
      data bth0cs(13) /   -.0000000000 02144188e0 /
      data bth0cs(14) /    .0000000000 00663693e0 /
      data bth0cs(15) /   -.0000000000 00215126e0 /
      data bth0cs(16) /    .0000000000 00072659e0 /
      data bth0cs(17) /   -.0000000000 00025465e0 /
      data bth0cs(18) /    .0000000000 00009229e0 /
      data bth0cs(19) /   -.0000000000 00003448e0 /
      data bth0cs(20) /    .0000000000 00001325e0 /
      data bth0cs(21) /   -.0000000000 00000522e0 /
      data bth0cs(22) /    .0000000000 00000210e0 /
      data bth0cs(23) /   -.0000000000 00000087e0 /
      data bth0cs(24) /    .0000000000 00000036e0 /
c
      data pi4 / 0.7853981633 9744831e0 /
      data ntj0, ntm0, ntth0, xsml, xmax / 3*0, 2*0./
c
      if (ntj0.ne.0) go to 10
      ntj0 = inits (bj0cs, 13, 0.1*r1mach(3))
      ntm0 = inits (bm0cs, 21, 0.1*r1mach(3))
      ntth0 = inits (bth0cs, 24, 0.1*r1mach(3))
c
      xsml = sqrt (4.0*r1mach(3))
      xmax = 1.0/r1mach(4)
c
 10   y = abs(x)
      if (y.gt.4.0) go to 20
c
      besj0 = 1.0
      if (y.gt.xsml) besj0 = csevl (.125*y*y-1., bj0cs, ntj0)
      return
c
 20   if (y.gt.xmax) call seteru (
     1  42hbesj0   no precision because abs(x) is big, 42, 1, 2)
c
      z = 32.0/y**2 - 1.0
      ampl = (0.75 + csevl (z, bm0cs, ntm0)) / sqrt(y)
      theta = y - pi4 + csevl (z, bth0cs, ntth0) / y
      besj0 = ampl * cos (theta)
c
      return
      end
      function besj1 (x)
c sept 1983 edition.  w. fullerton, c3, los alamos scientific lab.
      dimension bj1cs(12), bm1cs(21), bth1cs(24)
c     external cos, csevl, inits, r1mach, sqrt
c
c series for bj1        on the interval  0.          to  1.60000d+01
c                                        with weighted error   4.48e-17
c                                         log weighted error  16.35
c                               significant figures required  15.77
c                                    decimal places required  16.89
c
      data bj1 cs( 1) /   -.1172614151 3332787e0 /
      data bj1 cs( 2) /   -.2536152183 0790640e0 /
      data bj1 cs( 3) /    .0501270809 84469569e0 /
      data bj1 cs( 4) /   -.0046315148 09625081e0 /
      data bj1 cs( 5) /    .0002479962 29415914e0 /
      data bj1 cs( 6) /   -.0000086789 48686278e0 /
      data bj1 cs( 7) /    .0000002142 93917143e0 /
      data bj1 cs( 8) /   -.0000000039 36093079e0 /
      data bj1 cs( 9) /    .0000000000 55911823e0 /
      data bj1 cs(10) /   -.0000000000 00632761e0 /
      data bj1 cs(11) /    .0000000000 00005840e0 /
      data bj1 cs(12) /   -.0000000000 00000044e0 /
c
c series for bm1        on the interval  0.          to  6.25000d-02
c                                        with weighted error   5.61e-17
c                                         log weighted error  16.25
c                               significant figures required  14.97
c                                    decimal places required  16.91
c
      data bm1 cs( 1) /    .1047362510 931285e0 /
      data bm1 cs( 2) /    .0044244389 3702345e0 /
      data bm1 cs( 3) /   -.0000566163 9504035e0 /
      data bm1 cs( 4) /    .0000023134 9417339e0 /
      data bm1 cs( 5) /   -.0000001737 7182007e0 /
      data bm1 cs( 6) /    .0000000189 3209930e0 /
      data bm1 cs( 7) /   -.0000000026 5416023e0 /
      data bm1 cs( 8) /    .0000000004 4740209e0 /
      data bm1 cs( 9) /   -.0000000000 8691795e0 /
      data bm1 cs(10) /    .0000000000 1891492e0 /
      data bm1 cs(11) /   -.0000000000 0451884e0 /
      data bm1 cs(12) /    .0000000000 0116765e0 /
      data bm1 cs(13) /   -.0000000000 0032265e0 /
      data bm1 cs(14) /    .0000000000 0009450e0 /
      data bm1 cs(15) /   -.0000000000 0002913e0 /
      data bm1 cs(16) /    .0000000000 0000939e0 /
      data bm1 cs(17) /   -.0000000000 0000315e0 /
      data bm1 cs(18) /    .0000000000 0000109e0 /
      data bm1 cs(19) /   -.0000000000 0000039e0 /
      data bm1 cs(20) /    .0000000000 0000014e0 /
      data bm1 cs(21) /   -.0000000000 0000005e0 /
c
c series for bth1       on the interval  0.          to  6.25000d-02
c                                        with weighted error   4.10e-17
c                                         log weighted error  16.39
c                               significant figures required  15.96
c                                    decimal places required  17.08
c
      data bth1cs( 1) /    .7406014102 6313850e0 /
      data bth1cs( 2) /   -.0045717556 59637690e0 /
      data bth1cs( 3) /    .0001198185 10964326e0 /
      data bth1cs( 4) /   -.0000069645 61891648e0 /
      data bth1cs( 5) /    .0000006554 95621447e0 /
      data bth1cs( 6) /   -.0000000840 66228945e0 /
      data bth1cs( 7) /    .0000000133 76886564e0 /
      data bth1cs( 8) /   -.0000000024 99565654e0 /
      data bth1cs( 9) /    .0000000005 29495100e0 /
      data bth1cs(10) /   -.0000000001 24135944e0 /
      data bth1cs(11) /    .0000000000 31656485e0 /
      data bth1cs(12) /   -.0000000000 08668640e0 /
      data bth1cs(13) /    .0000000000 02523758e0 /
      data bth1cs(14) /   -.0000000000 00775085e0 /
      data bth1cs(15) /    .0000000000 00249527e0 /
      data bth1cs(16) /   -.0000000000 00083773e0 /
      data bth1cs(17) /    .0000000000 00029205e0 /
      data bth1cs(18) /   -.0000000000 00010534e0 /
      data bth1cs(19) /    .0000000000 00003919e0 /
      data bth1cs(20) /   -.0000000000 00001500e0 /
      data bth1cs(21) /    .0000000000 00000589e0 /
      data bth1cs(22) /   -.0000000000 00000237e0 /
      data bth1cs(23) /    .0000000000 00000097e0 /
      data bth1cs(24) /   -.0000000000 00000040e0 /
c
c
      data pi4 / 0.7853981633 9744831e0 /
      data ntj1, ntm1, ntth1, xsml, xmin, xmax / 3*0, 3*0./
c
      if (ntj1.ne.0) go to 10
      ntj1 = inits (bj1cs, 12, 0.1*r1mach(3))
      ntm1 = inits (bm1cs, 21, 0.1*r1mach(3))
      ntth1 = inits (bth1cs, 24, 0.1*r1mach(3))
c
      xsml = sqrt (8.0*r1mach(3))
      xmin = 2.0*r1mach(1)
      xmax = 1.0/r1mach(4)
c
 10   y = abs(x)
      if (y.gt.4.0) go to 20
c
      besj1 = 0.0
      if (y.eq.0.0) return
      if (y.le.xmin) call seteru (
     1  37hbesj1   abs(x) so small j1 underflows, 37, 1, 0)
      if (y.gt.xmin) besj1 = 0.5*x
      if (y.gt.xsml) besj1 = x * (.25 + csevl(.125*y*y-1., bj1cs, ntj1))
      return
c
 20   if (y.gt.xmax) call seteru (
     1  42hbesj1   no precision because abs(x) is big, 42, 2, 2)
      z = 32.0/y**2 - 1.0
      ampl = (0.75 + csevl (z, bm1cs, ntm1)) / sqrt(y)
      theta = y - 3.0*pi4 + csevl (z, bth1cs, ntth1) / y
      besj1 = sign (ampl, x) * cos (theta)
c
      return
      end
      function besk0 (x)
c april 1977 version.  w. fullerton, c3, los alamos scientific lab.
      dimension bk0cs(11)
c     external alog, besi0, besk0e, csevl, exp, inits, r1mach, sqrt
c
c series for bk0        on the interval  0.          to  4.00000d+00
c                                        with weighted error   3.57e-19
c                                         log weighted error  18.45
c                               significant figures required  17.99
c                                    decimal places required  18.97
c
      data bk0 cs( 1) /   -.0353273932 3390276872e0 /
      data bk0 cs( 2) /    .3442898999 246284869e0 /
      data bk0 cs( 3) /    .0359799365 1536150163e0 /
      data bk0 cs( 4) /    .0012646154 1144692592e0 /
      data bk0 cs( 5) /    .0000228621 2103119451e0 /
      data bk0 cs( 6) /    .0000002534 7910790261e0 /
      data bk0 cs( 7) /    .0000000019 0451637722e0 /
      data bk0 cs( 8) /    .0000000000 1034969525e0 /
      data bk0 cs( 9) /    .0000000000 0004259816e0 /
      data bk0 cs(10) /    .0000000000 0000013744e0 /
      data bk0 cs(11) /    .0000000000 0000000035e0 /
c
      data ntk0, xsml, xmax / 0, 0., 0. /
c
      if (ntk0.ne.0) go to 10
      ntk0 = inits (bk0cs, 11, 0.1*r1mach(3))
      xsml = sqrt (4.0*r1mach(3))
      xmax = -alog(r1mach(1))
      xmax = xmax - 0.5*xmax*alog(xmax)/(xmax+0.5) - 0.01
c
 10   if (x.le.0.) call seteru (29hbesk0   x is zero or negative, 29,
     1  2, 2)
      if (x.gt.2.) go to 20
c
      y = 0.
      if (x.gt.xsml) y = x*x
      besk0 = -alog(0.5*x)*besi0(x) - .25 + csevl (.5*y-1., bk0cs, ntk0)
      return
c
 20   besk0 = 0.
      if (x.gt.xmax) call seteru (30hbesk0   x so big k0 underflows, 30,
     1  1, 0)
      if (x.gt.xmax) return
c
      besk0 = exp(-x) * besk0e(x)
c
      return
      end
      function besk1 (x)
c july 1980 version.  w. fullerton, c3, los alamos scientific lab.
      dimension bk1cs(11)
c     external alog, besi1, besk1e, csevl, exp, inits, r1mach, sqrt
c
c series for bk1        on the interval  0.          to  4.00000d+00
c                                        with weighted error   7.02e-18
c                                         log weighted error  17.15
c                               significant figures required  16.73
c                                    decimal places required  17.67
c
      data bk1 cs( 1) /    .0253002273 389477705e0 /
      data bk1 cs( 2) /   -.3531559607 76544876e0 /
      data bk1 cs( 3) /   -.1226111808 22657148e0 /
      data bk1 cs( 4) /   -.0069757238 596398643e0 /
      data bk1 cs( 5) /   -.0001730288 957513052e0 /
      data bk1 cs( 6) /   -.0000024334 061415659e0 /
      data bk1 cs( 7) /   -.0000000221 338763073e0 /
      data bk1 cs( 8) /   -.0000000001 411488392e0 /
      data bk1 cs( 9) /   -.0000000000 006666901e0 /
      data bk1 cs(10) /   -.0000000000 000024274e0 /
      data bk1 cs(11) /   -.0000000000 000000070e0 /
c
      data ntk1, xmin, xsml, xmax /0, 0., 0., 0. /
c
      if (ntk1.ne.0) go to 10
      ntk1 = inits (bk1cs, 11, 0.1*r1mach(3))
      xmin = exp (amax1(alog(r1mach(1)), -alog(r1mach(2))) + .01)
      xsml = sqrt (4.0*r1mach(3))
      xmax = -alog(r1mach(1))
      xmax = xmax - 0.5*xmax*alog(xmax)/(xmax+0.5) - 0.01
c
 10   if (x.le.0.) call seteru (29hbesk1   x is zero or negative, 29,
     1  2, 2)
      if (x.gt.2.0) go to 20
c
      if (x.lt.xmin) call seteru (31hbesk1   x so small k1 overflows,
     1  31, 3, 2)
      y = 0.
      if (x.gt.xsml) y = x*x
      besk1 = alog(0.5*x)*besi1(x) +
     1  (0.75 + csevl (.5*y-1., bk1cs, ntk1))/x
      return
c
 20   besk1 = 0.
      if (x.gt.xmax) call seteru (30hbesk1   x so big k1 underflows,
     1  30, 1, 0)
      if (x.gt.xmax) return
c
      besk1 = exp(-x) * besk1e(x)
c
      return
      end
      function besy0 (x)
c august 1980 version.  w. fullerton, c3, los alamos scientific lab.
      dimension by0cs(13), bm0cs(21), bth0cs(24)
c     external alog, besj0, csevl, inits, r1mach, sin, sqrt
c
c series for by0        on the interval  0.          to  1.60000d+01
c                                        with weighted error   1.20e-17
c                                         log weighted error  16.92
c                               significant figures required  16.15
c                                    decimal places required  17.48
c
      data by0 cs( 1) /   -.0112778393 92865573e0 /
      data by0 cs( 2) /   -.1283452375 6042035e0 /
      data by0 cs( 3) /   -.1043788479 9794249e0 /
      data by0 cs( 4) /    .0236627491 83969695e0 /
      data by0 cs( 5) /   -.0020903916 47700486e0 /
      data by0 cs( 6) /    .0001039754 53939057e0 /
      data by0 cs( 7) /   -.0000033697 47162423e0 /
      data by0 cs( 8) /    .0000000772 93842676e0 /
      data by0 cs( 9) /   -.0000000013 24976772e0 /
      data by0 cs(10) /    .0000000000 17648232e0 /
      data by0 cs(11) /   -.0000000000 00188105e0 /
      data by0 cs(12) /    .0000000000 00001641e0 /
      data by0 cs(13) /   -.0000000000 00000011e0 /
c
c series for bm0        on the interval  0.          to  6.25000d-02
c                                        with weighted error   4.98e-17
c                                         log weighted error  16.30
c                               significant figures required  14.97
c                                    decimal places required  16.96
c
      data bm0 cs( 1) /    .0928496163 7381644e0 /
      data bm0 cs( 2) /   -.0014298770 7403484e0 /
      data bm0 cs( 3) /    .0000283057 9271257e0 /
      data bm0 cs( 4) /   -.0000014330 0611424e0 /
      data bm0 cs( 5) /    .0000001202 8628046e0 /
      data bm0 cs( 6) /   -.0000000139 7113013e0 /
      data bm0 cs( 7) /    .0000000020 4076188e0 /
      data bm0 cs( 8) /   -.0000000003 5399669e0 /
      data bm0 cs( 9) /    .0000000000 7024759e0 /
      data bm0 cs(10) /   -.0000000000 1554107e0 /
      data bm0 cs(11) /    .0000000000 0376226e0 /
      data bm0 cs(12) /   -.0000000000 0098282e0 /
      data bm0 cs(13) /    .0000000000 0027408e0 /
      data bm0 cs(14) /   -.0000000000 0008091e0 /
      data bm0 cs(15) /    .0000000000 0002511e0 /
      data bm0 cs(16) /   -.0000000000 0000814e0 /
      data bm0 cs(17) /    .0000000000 0000275e0 /
      data bm0 cs(18) /   -.0000000000 0000096e0 /
      data bm0 cs(19) /    .0000000000 0000034e0 /
      data bm0 cs(20) /   -.0000000000 0000012e0 /
      data bm0 cs(21) /    .0000000000 0000004e0 /
c
c series for bth0       on the interval  0.          to  6.25000d-02
c                                        with weighted error   3.67e-17
c                                         log weighted error  16.44
c                               significant figures required  15.53
c                                    decimal places required  17.13
c
      data bth0cs( 1) /   -.2463916377 4300119e0 /
      data bth0cs( 2) /    .0017370983 07508963e0 /
      data bth0cs( 3) /   -.0000621836 33402968e0 /
      data bth0cs( 4) /    .0000043680 50165742e0 /
      data bth0cs( 5) /   -.0000004560 93019869e0 /
      data bth0cs( 6) /    .0000000621 97400101e0 /
      data bth0cs( 7) /   -.0000000103 00442889e0 /
      data bth0cs( 8) /    .0000000019 79526776e0 /
      data bth0cs( 9) /   -.0000000004 28198396e0 /
      data bth0cs(10) /    .0000000001 02035840e0 /
      data bth0cs(11) /   -.0000000000 26363898e0 /
      data bth0cs(12) /    .0000000000 07297935e0 /
      data bth0cs(13) /   -.0000000000 02144188e0 /
      data bth0cs(14) /    .0000000000 00663693e0 /
      data bth0cs(15) /   -.0000000000 00215126e0 /
      data bth0cs(16) /    .0000000000 00072659e0 /
      data bth0cs(17) /   -.0000000000 00025465e0 /
      data bth0cs(18) /    .0000000000 00009229e0 /
      data bth0cs(19) /   -.0000000000 00003448e0 /
      data bth0cs(20) /    .0000000000 00001325e0 /
      data bth0cs(21) /   -.0000000000 00000522e0 /
      data bth0cs(22) /    .0000000000 00000210e0 /
      data bth0cs(23) /   -.0000000000 00000087e0 /
      data bth0cs(24) /    .0000000000 00000036e0 /
c
c twodpi = 2.0/pi
      data twodpi / 0.6366197723 6758134e0 /
      data pi4 / 0.7853981633 9744831e0 /
      data alnhaf / -0.6931471805 59945309e0 /
      data nty0, ntm0, ntth0, xsml, xmax / 3*0, 2*0./
c
      if (nty0.ne.0) go to 10
      nty0 = inits (by0cs, 13, 0.1*r1mach(3))
      ntm0 = inits (bm0cs, 21, 0.1*r1mach(3))
      ntth0 = inits (bth0cs, 24, 0.1*r1mach(3))
c
      xsml = sqrt (4.0*r1mach(3))
      xmax = 1.0/r1mach(4)
c
 10   if (x.le.0.) call seteru (29hbesy0   x is zero or negative, 29,
     1  1, 2)
      if (x.gt.4.0) go to 20
c
      y = 0.
      if (x.gt.xsml) y = x*x
      besy0 = twodpi*(alnhaf+alog(x))*besj0(x) + .375 +
     1  csevl (.125*y-1., by0cs, nty0)
      return
c
 20   if (x.gt.xmax) call seteru (
     1  37hbesy0   no precision because x is big, 37, 2, 2)
c
      z = 32.0/x**2 - 1.0
      ampl = (0.75 + csevl (z, bm0cs, ntm0)) / sqrt(x)
      theta = x - pi4 + csevl (z, bth0cs, ntth0) / x
      besy0 = ampl * sin (theta)
c
      return
      end
      function besy1 (x)
c april 1977 version.  w. fullerton, c3, los alamos scientific lab.
      dimension by1cs(14), bm1cs(21), bth1cs(24)
c     external alog, besj1, csevl, exp, inits, r1mach, sin, sqrt
c
c series for by1        on the interval  0.          to  1.60000d+01
c                                        with weighted error   1.87e-18
c                                         log weighted error  17.73
c                               significant figures required  17.83
c                                    decimal places required  18.30
c
      data by1 cs( 1) /    .0320804710 0611908629e0 /
      data by1 cs( 2) /   1.2627078974 33500450e0 /
      data by1 cs( 3) /    .0064999618 9992317500e0 /
      data by1 cs( 4) /   -.0893616452 8860504117e0 /
      data by1 cs( 5) /    .0132508812 2175709545e0 /
      data by1 cs( 6) /   -.0008979059 1196483523e0 /
      data by1 cs( 7) /    .0000364736 1487958306e0 /
      data by1 cs( 8) /   -.0000010013 7438166600e0 /
      data by1 cs( 9) /    .0000000199 4539657390e0 /
      data by1 cs(10) /   -.0000000003 0230656018e0 /
      data by1 cs(11) /    .0000000000 0360987815e0 /
      data by1 cs(12) /   -.0000000000 0003487488e0 /
      data by1 cs(13) /    .0000000000 0000027838e0 /
      data by1 cs(14) /   -.0000000000 0000000186e0 /
c
c series for bm1        on the interval  0.          to  6.25000d-02
c                                        with weighted error   5.61e-17
c                                         log weighted error  16.25
c                               significant figures required  14.97
c                                    decimal places required  16.91
c
      data bm1 cs( 1) /    .1047362510 931285e0 /
      data bm1 cs( 2) /    .0044244389 3702345e0 /
      data bm1 cs( 3) /   -.0000566163 9504035e0 /
      data bm1 cs( 4) /    .0000023134 9417339e0 /
      data bm1 cs( 5) /   -.0000001737 7182007e0 /
      data bm1 cs( 6) /    .0000000189 3209930e0 /
      data bm1 cs( 7) /   -.0000000026 5416023e0 /
      data bm1 cs( 8) /    .0000000004 4740209e0 /
      data bm1 cs( 9) /   -.0000000000 8691795e0 /
      data bm1 cs(10) /    .0000000000 1891492e0 /
      data bm1 cs(11) /   -.0000000000 0451884e0 /
      data bm1 cs(12) /    .0000000000 0116765e0 /
      data bm1 cs(13) /   -.0000000000 0032265e0 /
      data bm1 cs(14) /    .0000000000 0009450e0 /
      data bm1 cs(15) /   -.0000000000 0002913e0 /
      data bm1 cs(16) /    .0000000000 0000939e0 /
      data bm1 cs(17) /   -.0000000000 0000315e0 /
      data bm1 cs(18) /    .0000000000 0000109e0 /
      data bm1 cs(19) /   -.0000000000 0000039e0 /
      data bm1 cs(20) /    .0000000000 0000014e0 /
      data bm1 cs(21) /   -.0000000000 0000005e0 /
c
c series for bth1       on the interval  0.          to  6.25000d-02
c                                        with weighted error   4.10e-17
c                                         log weighted error  16.39
c                               significant figures required  15.96
c                                    decimal places required  17.08
c
      data bth1cs( 1) /    .7406014102 6313850e0 /
      data bth1cs( 2) /   -.0045717556 59637690e0 /
      data bth1cs( 3) /    .0001198185 10964326e0 /
      data bth1cs( 4) /   -.0000069645 61891648e0 /
      data bth1cs( 5) /    .0000006554 95621447e0 /
      data bth1cs( 6) /   -.0000000840 66228945e0 /
      data bth1cs( 7) /    .0000000133 76886564e0 /
      data bth1cs( 8) /   -.0000000024 99565654e0 /
      data bth1cs( 9) /    .0000000005 29495100e0 /
      data bth1cs(10) /   -.0000000001 24135944e0 /
      data bth1cs(11) /    .0000000000 31656485e0 /
      data bth1cs(12) /   -.0000000000 08668640e0 /
      data bth1cs(13) /    .0000000000 02523758e0 /
      data bth1cs(14) /   -.0000000000 00775085e0 /
      data bth1cs(15) /    .0000000000 00249527e0 /
      data bth1cs(16) /   -.0000000000 00083773e0 /
      data bth1cs(17) /    .0000000000 00029205e0 /
      data bth1cs(18) /   -.0000000000 00010534e0 /
      data bth1cs(19) /    .0000000000 00003919e0 /
      data bth1cs(20) /   -.0000000000 00001500e0 /
      data bth1cs(21) /    .0000000000 00000589e0 /
      data bth1cs(22) /   -.0000000000 00000237e0 /
      data bth1cs(23) /    .0000000000 00000097e0 /
      data bth1cs(24) /   -.0000000000 00000040e0 /
c
c twodpi = 2.0/pi
      data twodpi / 0.6366197723 6758134e0 /
      data pi4 / 0.7853981633 9744831e0 /
      data nty1, ntm1, ntth1, xmin, xsml, xmax / 3*0, 3*0./
c
      if (nty1.ne.0) go to 10
      nty1 = inits (by1cs, 14, 0.1*r1mach(3))
      ntm1 = inits (bm1cs, 21, 0.1*r1mach(3))
      ntth1 = inits (bth1cs, 24, 0.1*r1mach(3))
c
      xmin = 1.571*exp ( amax1(alog(r1mach(1)), -alog(r1mach(2)))+.01)
      xsml = sqrt (4.0*r1mach(3))
      xmax = 1.0/r1mach(4)
c
 10   if (x.le.0.) call seteru (29hbesy1   x is zero or negative, 29,
     1  1, 2)
      if (x.gt.4.0) go to 20
c
      if (x.lt.xmin) call seteru (31hbesy1   x so small y1 overflows,
     1  31, 3, 2)
      y = 0.
      if (x.gt.xsml) y = x*x
      besy1 = twodpi*alog(0.5*x)*besj1(x) +
     1  (0.5 + csevl (.125*y-1., by1cs, nty1))/x
      return
c
 20   if (x.gt.xmax) call seteru (
     1  37hbesy1   no precision because x is big, 37, 2, 2)
c
      z = 32.0/x**2 - 1.0
      ampl = (0.75 + csevl (z, bm1cs, ntm1)) / sqrt(x)
      theta = x - 3.0*pi4 + csevl (z, bth1cs, ntth1) / x
      besy1 = ampl * sin (theta)
c
      return
      end
      function inits (os, nos, eta)
c april 1977 version.  w. fullerton, c3, los alamos scientific lab.
c
c initialize the orthogonal series so that inits is the number of terms
c needed to insure the error is no larger than eta.  ordinarily, eta
c will be chosen to be one-tenth machine precision.
c
c             input arguments --
c os     array of nos coefficients in an orthogonal series.
c nos    number of coefficients in os.
c eta    requested accuracy of series.
c
      dimension os(nos)
c
      if (nos.lt.1) call seteru (
     1  35hinits   number of coefficients lt 1, 35, 2, 2)
c
      err = 0.
      do 10 ii=1,nos
        i = nos + 1 - ii
        err = err + abs(os(i))
        if (err.gt.eta) go to 20
 10   continue
c
 20   if (i.eq.nos) call seteru (28hinits   eta may be too small, 28,
     1  1, 2)
      inits = i
c
      return
      end
      function csevl (x, cs, n)
c april 1977 version.  w. fullerton, c3, los alamos scientific lab.
c
c evaluate the n-term chebyshev series cs at x.  adapted from
c r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973).  also see fox
c and parker, chebyshev polys in numerical analysis, oxford press, p.56.
c
c             input arguments --
c x      value at which the series is to be evaluated.
c cs     array of n terms of a chebyshev series.  in eval-
c        uating cs, only half the first coef is summed.
c n      number of terms in array cs.
c
      dimension cs(1)
c
      if (n.lt.1) call seteru (28hcsevl   number of terms le 0, 28, 2,2)
      if (n.gt.1000) call seteru (31hcsevl   number of terms gt 1000,
     1  31, 3, 2)
      if (x.lt.(-1.1) .or. x.gt.1.1) call seteru (
     1  25hcsevl   x outside (-1,+1), 25, 1, 1)
c
      b1 = 0.
      b0 = 0.
      twox = 2.*x
      do 10 i=1,n
        b2 = b1
        b1 = b0
        ni = n + 1 - i
        b0 = twox*b1 - b2 + cs(ni)
 10   continue
c
      csevl = 0.5 * (b0-b2)
c
      return
      end
      subroutine seteru (messg, nmessg, nerr, iopt)
      common /cseter/ iunflo
      integer messg(1)
      data iunflo / 0 /
c
      if (iopt.ne.0) call seterr (messg, nmessg, nerr, iopt)
      if (iopt.ne.0) return
c
      if (iunflo.le.0) return
      call seterr (messg, nmessg, nerr, 1)
c
      return
      end
      function besi0e (x)
c july 1977 version.  w. fullerton, c3, los alamos scientific lab.
      dimension bi0cs(12), ai0cs(21), ai02cs(22)
c     external csevl, exp, inits, r1mach, sqrt
c
c series for bi0        on the interval  0.          to  9.00000d+00
c                                        with weighted error   2.46e-18
c                                         log weighted error  17.61
c                               significant figures required  17.90
c                                    decimal places required  18.15
c
      data bi0 cs( 1) /   -.0766054725 2839144951e0 /
      data bi0 cs( 2) /   1.9273379539 93808270e0 /
      data bi0 cs( 3) /    .2282644586 920301339e0 /
      data bi0 cs( 4) /    .0130489146 6707290428e0 /
      data bi0 cs( 5) /    .0004344270 9008164874e0 /
      data bi0 cs( 6) /    .0000094226 5768600193e0 /
      data bi0 cs( 7) /    .0000001434 0062895106e0 /
      data bi0 cs( 8) /    .0000000016 1384906966e0 /
      data bi0 cs( 9) /    .0000000000 1396650044e0 /
      data bi0 cs(10) /    .0000000000 0009579451e0 /
      data bi0 cs(11) /    .0000000000 0000053339e0 /
      data bi0 cs(12) /    .0000000000 0000000245e0 /
c
c series for ai0        on the interval  1.25000d-01 to  3.33333d-01
c                                        with weighted error   7.87e-17
c                                         log weighted error  16.10
c                               significant figures required  14.69
c                                    decimal places required  16.76
c
      data ai0 cs( 1) /    .0757599449 4023796e0 /
      data ai0 cs( 2) /    .0075913808 1082334e0 /
      data ai0 cs( 3) /    .0004153131 3389237e0 /
      data ai0 cs( 4) /    .0000107007 6463439e0 /
      data ai0 cs( 5) /   -.0000079011 7997921e0 /
      data ai0 cs( 6) /   -.0000007826 1435014e0 /
      data ai0 cs( 7) /    .0000002783 8499429e0 /
      data ai0 cs( 8) /    .0000000082 5247260e0 /
      data ai0 cs( 9) /   -.0000000120 4463945e0 /
      data ai0 cs(10) /    .0000000015 5964859e0 /
      data ai0 cs(11) /    .0000000002 2925563e0 /
      data ai0 cs(12) /   -.0000000001 1916228e0 /
      data ai0 cs(13) /    .0000000000 1757854e0 /
      data ai0 cs(14) /    .0000000000 0112822e0 /
      data ai0 cs(15) /   -.0000000000 0114684e0 /
      data ai0 cs(16) /    .0000000000 0027155e0 /
      data ai0 cs(17) /   -.0000000000 0002415e0 /
      data ai0 cs(18) /   -.0000000000 0000608e0 /
      data ai0 cs(19) /    .0000000000 0000314e0 /
      data ai0 cs(20) /   -.0000000000 0000071e0 /
      data ai0 cs(21) /    .0000000000 0000007e0 /
c
c series for ai02       on the interval  0.          to  1.25000d-01
c                                        with weighted error   3.79e-17
c                                         log weighted error  16.42
c                               significant figures required  14.86
c                                    decimal places required  17.09
c
      data ai02cs( 1) /    .0544904110 1410882e0 /
      data ai02cs( 2) /    .0033691164 7825569e0 /
      data ai02cs( 3) /    .0000688975 8346918e0 /
      data ai02cs( 4) /    .0000028913 7052082e0 /
      data ai02cs( 5) /    .0000002048 9185893e0 /
      data ai02cs( 6) /    .0000000226 6668991e0 /
      data ai02cs( 7) /    .0000000033 9623203e0 /
      data ai02cs( 8) /    .0000000004 9406022e0 /
      data ai02cs( 9) /    .0000000000 1188914e0 /
      data ai02cs(10) /   -.0000000000 3149915e0 /
      data ai02cs(11) /   -.0000000000 1321580e0 /
      data ai02cs(12) /   -.0000000000 0179419e0 /
      data ai02cs(13) /    .0000000000 0071801e0 /
      data ai02cs(14) /    .0000000000 0038529e0 /
      data ai02cs(15) /    .0000000000 0001539e0 /
      data ai02cs(16) /   -.0000000000 0004151e0 /
      data ai02cs(17) /   -.0000000000 0000954e0 /
      data ai02cs(18) /    .0000000000 0000382e0 /
      data ai02cs(19) /    .0000000000 0000176e0 /
      data ai02cs(20) /   -.0000000000 0000034e0 /
      data ai02cs(21) /   -.0000000000 0000027e0 /
      data ai02cs(22) /    .0000000000 0000003e0 /
c
      data nti0, ntai0, ntai02, xsml / 3*0, 0. /
c
      if (nti0.ne.0) go to 10
      nti0 = inits (bi0cs, 12, 0.1*r1mach(3))
      ntai0 = inits (ai0cs, 21, 0.1*r1mach(3))
      ntai02 = inits (ai02cs, 22, 0.1*r1mach(3))
      xsml = sqrt (4.0*r1mach(3))
c
 10   y = abs(x)
      if (y.gt.3.0) go to 20
c
      besi0e = 1.0
      if (y.gt.xsml) besi0e = exp(-y) * ( 2.75 +
     1  csevl (y*y/4.5-1.0, bi0cs, nti0) )
      return
c
 20   if (y.le.8.) besi0e = (.375 + csevl ((48./y-11.)/5., ai0cs, ntai0)
     1  ) / sqrt(y)
      if (y.gt.8.) besi0e = (.375 + csevl (16./y-1., ai02cs, ntai02))
     1  / sqrt(y)
c
      return
      end
      function besi1e (x)
c oct 1983 version.  w. fullerton, c3, los alamos scientific lab.
      dimension bi1cs(11), ai1cs(21), ai12cs(22)
c     external csevl, exp, inits, r1mach, sqrt
c
c series for bi1        on the interval  0.          to  9.00000d+00
c                                        with weighted error   2.40e-17
c                                         log weighted error  16.62
c                               significant figures required  16.23
c                                    decimal places required  17.14
c
      data bi1 cs( 1) /   -.0019717132 61099859e0 /
      data bi1 cs( 2) /    .4073488766 7546481e0 /
      data bi1 cs( 3) /    .0348389942 99959456e0 /
      data bi1 cs( 4) /    .0015453945 56300123e0 /
      data bi1 cs( 5) /    .0000418885 21098377e0 /
      data bi1 cs( 6) /    .0000007649 02676483e0 /
      data bi1 cs( 7) /    .0000000100 42493924e0 /
      data bi1 cs( 8) /    .0000000000 99322077e0 /
      data bi1 cs( 9) /    .0000000000 00766380e0 /
      data bi1 cs(10) /    .0000000000 00004741e0 /
      data bi1 cs(11) /    .0000000000 00000024e0 /
c
c series for ai1        on the interval  1.25000d-01 to  3.33333d-01
c                                        with weighted error   6.98e-17
c                                         log weighted error  16.16
c                               significant figures required  14.53
c                                    decimal places required  16.82
c
      data ai1 cs( 1) /   -.0284674418 1881479e0 /
      data ai1 cs( 2) /   -.0192295323 1443221e0 /
      data ai1 cs( 3) /   -.0006115185 8579437e0 /
      data ai1 cs( 4) /   -.0000206997 1253350e0 /
      data ai1 cs( 5) /    .0000085856 1914581e0 /
      data ai1 cs( 6) /    .0000010494 9824671e0 /
      data ai1 cs( 7) /   -.0000002918 3389184e0 /
      data ai1 cs( 8) /   -.0000000155 9378146e0 /
      data ai1 cs( 9) /    .0000000131 8012367e0 /
      data ai1 cs(10) /   -.0000000014 4842341e0 /
      data ai1 cs(11) /   -.0000000002 9085122e0 /
      data ai1 cs(12) /    .0000000001 2663889e0 /
      data ai1 cs(13) /   -.0000000000 1664947e0 /
      data ai1 cs(14) /   -.0000000000 0166665e0 /
      data ai1 cs(15) /    .0000000000 0124260e0 /
      data ai1 cs(16) /   -.0000000000 0027315e0 /
      data ai1 cs(17) /    .0000000000 0002023e0 /
      data ai1 cs(18) /    .0000000000 0000730e0 /
      data ai1 cs(19) /   -.0000000000 0000333e0 /
      data ai1 cs(20) /    .0000000000 0000071e0 /
      data ai1 cs(21) /   -.0000000000 0000006e0 /
c
c series for ai12       on the interval  0.          to  1.25000d-01
c                                        with weighted error   3.55e-17
c                                         log weighted error  16.45
c                               significant figures required  14.69
c                                    decimal places required  17.12
c
      data ai12cs( 1) /    .0285762350 1828014e0 /
      data ai12cs( 2) /   -.0097610974 9136147e0 /
      data ai12cs( 3) /   -.0001105889 3876263e0 /
      data ai12cs( 4) /   -.0000038825 6480887e0 /
      data ai12cs( 5) /   -.0000002512 2362377e0 /
      data ai12cs( 6) /   -.0000000263 1468847e0 /
      data ai12cs( 7) /   -.0000000038 3538039e0 /
      data ai12cs( 8) /   -.0000000005 5897433e0 /
      data ai12cs( 9) /   -.0000000000 1897495e0 /
      data ai12cs(10) /    .0000000000 3252602e0 /
      data ai12cs(11) /    .0000000000 1412580e0 /
      data ai12cs(12) /    .0000000000 0203564e0 /
      data ai12cs(13) /   -.0000000000 0071985e0 /
      data ai12cs(14) /   -.0000000000 0040836e0 /
      data ai12cs(15) /   -.0000000000 0002101e0 /
      data ai12cs(16) /    .0000000000 0004273e0 /
      data ai12cs(17) /    .0000000000 0001041e0 /
      data ai12cs(18) /   -.0000000000 0000382e0 /
      data ai12cs(19) /   -.0000000000 0000186e0 /
      data ai12cs(20) /    .0000000000 0000033e0 /
      data ai12cs(21) /    .0000000000 0000028e0 /
      data ai12cs(22) /   -.0000000000 0000003e0 /
c
      data nti1, ntai1, ntai12, xmin, xsml / 3*0, 2*0. /
c
      if (nti1.ne.0) go to 10
      nti1 = inits (bi1cs, 11, 0.1*r1mach(3))
      ntai1 = inits (ai1cs, 21, 0.1*r1mach(3))
      ntai12 = inits (ai12cs, 22, 0.1*r1mach(3))
c
      xmin = 2.0*r1mach(1)
      xsml = sqrt (8.0*r1mach(3))
c
 10   y = abs(x)
      if (y.gt.3.0) go to 20
c
      besi1e = 0.0
      if (y.eq.0.0) return
c
      if (y.le.xmin) call seteru (
     1  37hbesi1e  abs(x) so small i1 underflows, 37, 1, 0)
      besi1e = 0.0
      if (y.gt.xmin) besi1e = 0.5*x
      if (y.gt.xsml) besi1e = x * (.875 + csevl(y*y/4.5-1., bi1cs,nti1))
      besi1e = exp(-y) * besi1e
      return
c
 20   if (y.le.8.) besi1e = (.375 + csevl ((48./y-11.)/5., ai1cs, ntai1)
     1  ) / sqrt(y)
      if (y.gt.8.) besi1e = (.375 + csevl (16./y-1.0, ai12cs, ntai12))
     1  / sqrt(y)
      besi1e = sign (besi1e, x)
c
      return
      end
      function besk0e (x)
c july 1980 version.  w. fullerton, c3, los alamos scientific lab.
      dimension bk0cs(11), ak0cs(17), ak02cs(14)
c     external alog, besi0, csevl, exp, inits, r1mach, sqrt
c
c series for bk0        on the interval  0.          to  4.00000d+00
c                                        with weighted error   3.57e-19
c                                         log weighted error  18.45
c                               significant figures required  17.99
c                                    decimal places required  18.97
c
      data bk0 cs( 1) /   -.0353273932 3390276872e0 /
      data bk0 cs( 2) /    .3442898999 246284869e0 /
      data bk0 cs( 3) /    .0359799365 1536150163e0 /
      data bk0 cs( 4) /    .0012646154 1144692592e0 /
      data bk0 cs( 5) /    .0000228621 2103119451e0 /
      data bk0 cs( 6) /    .0000002534 7910790261e0 /
      data bk0 cs( 7) /    .0000000019 0451637722e0 /
      data bk0 cs( 8) /    .0000000000 1034969525e0 /
      data bk0 cs( 9) /    .0000000000 0004259816e0 /
      data bk0 cs(10) /    .0000000000 0000013744e0 /
      data bk0 cs(11) /    .0000000000 0000000035e0 /
c
c series for ak0        on the interval  1.25000d-01 to  5.00000d-01
c                                        with weighted error   5.34e-17
c                                         log weighted error  16.27
c                               significant figures required  14.92
c                                    decimal places required  16.89
c
      data ak0 cs( 1) /   -.0764394790 3327941e0 /
      data ak0 cs( 2) /   -.0223565260 5699819e0 /
      data ak0 cs( 3) /    .0007734181 1546938e0 /
      data ak0 cs( 4) /   -.0000428100 6688886e0 /
      data ak0 cs( 5) /    .0000030817 0017386e0 /
      data ak0 cs( 6) /   -.0000002639 3672220e0 /
      data ak0 cs( 7) /    .0000000256 3713036e0 /
      data ak0 cs( 8) /   -.0000000027 4270554e0 /
      data ak0 cs( 9) /    .0000000003 1694296e0 /
      data ak0 cs(10) /   -.0000000000 3902353e0 /
      data ak0 cs(11) /    .0000000000 0506804e0 /
      data ak0 cs(12) /   -.0000000000 0068895e0 /
      data ak0 cs(13) /    .0000000000 0009744e0 /
      data ak0 cs(14) /   -.0000000000 0001427e0 /
      data ak0 cs(15) /    .0000000000 0000215e0 /
      data ak0 cs(16) /   -.0000000000 0000033e0 /
      data ak0 cs(17) /    .0000000000 0000005e0 /
c
c series for ak02       on the interval  0.          to  1.25000d-01
c                                        with weighted error   2.34e-17
c                                         log weighted error  16.63
c                               significant figures required  14.67
c                                    decimal places required  17.20
c
      data ak02cs( 1) /   -.0120186982 6307592e0 /
      data ak02cs( 2) /   -.0091748526 9102569e0 /
      data ak02cs( 3) /    .0001444550 9317750e0 /
      data ak02cs( 4) /   -.0000040136 1417543e0 /
      data ak02cs( 5) /    .0000001567 8318108e0 /
      data ak02cs( 6) /   -.0000000077 7011043e0 /
      data ak02cs( 7) /    .0000000004 6111825e0 /
      data ak02cs( 8) /   -.0000000000 3158592e0 /
      data ak02cs( 9) /    .0000000000 0243501e0 /
      data ak02cs(10) /   -.0000000000 0020743e0 /
      data ak02cs(11) /    .0000000000 0001925e0 /
      data ak02cs(12) /   -.0000000000 0000192e0 /
      data ak02cs(13) /    .0000000000 0000020e0 /
      data ak02cs(14) /   -.0000000000 0000002e0 /
c
      data ntk0, ntak0, ntak02, xsml / 3*0, 0. /
c
      if (ntk0.ne.0) go to 10
      ntk0 = inits (bk0cs, 11, 0.1*r1mach(3))
      ntak0 = inits (ak0cs, 17, 0.1*r1mach(3))
      ntak02 = inits (ak02cs, 14, 0.1*r1mach(3))
      xsml = sqrt (4.0*r1mach(3))
c
 10   if (x.le.0.) call seteru (29hbesk0e  x is zero or negative, 29,
     1  1, 2)
      if (x.gt.2.) go to 20
c
      y = 0.
      if (x.gt.xsml) y = x*x
      besk0e = exp(x) * (-alog(0.5*x)*besi0(x)
     1  - .25 + csevl (.5*y-1., bk0cs, ntk0) )
      return
c
 20   if (x.le.8.) besk0e = (1.25 + csevl ((16./x-5.)/3., ak0cs, ntak0))
     1  / sqrt(x)
      if (x.gt.8.) besk0e = (1.25 + csevl (16./x-1., ak02cs, ntak02))
     1  / sqrt(x)
c
      return
      end
      function besk1e (x)
c july 1980 version.  w. fullerton, c3, los alamos scientific lab.
      dimension bk1cs(11), ak1cs(17), ak12cs(14)
c     external alog, besi1, csevl, exp, inits, r1mach, sqrt
c
c series for bk1        on the interval  0.          to  4.00000d+00
c                                        with weighted error   7.02e-18
c                                         log weighted error  17.15
c                               significant figures required  16.73
c                                    decimal places required  17.67
c
      data bk1 cs( 1) /    .0253002273 389477705e0 /
      data bk1 cs( 2) /   -.3531559607 76544876e0 /
      data bk1 cs( 3) /   -.1226111808 22657148e0 /
      data bk1 cs( 4) /   -.0069757238 596398643e0 /
      data bk1 cs( 5) /   -.0001730288 957513052e0 /
      data bk1 cs( 6) /   -.0000024334 061415659e0 /
      data bk1 cs( 7) /   -.0000000221 338763073e0 /
      data bk1 cs( 8) /   -.0000000001 411488392e0 /
      data bk1 cs( 9) /   -.0000000000 006666901e0 /
      data bk1 cs(10) /   -.0000000000 000024274e0 /
      data bk1 cs(11) /   -.0000000000 000000070e0 /
c
c series for ak1        on the interval  1.25000d-01 to  5.00000d-01
c                                        with weighted error   6.06e-17
c                                         log weighted error  16.22
c                               significant figures required  15.41
c                                    decimal places required  16.83
c
      data ak1 cs( 1) /    .2744313406 973883e0 /
      data ak1 cs( 2) /    .0757198995 3199368e0 /
      data ak1 cs( 3) /   -.0014410515 5647540e0 /
      data ak1 cs( 4) /    .0000665011 6955125e0 /
      data ak1 cs( 5) /   -.0000043699 8470952e0 /
      data ak1 cs( 6) /    .0000003540 2774997e0 /
      data ak1 cs( 7) /   -.0000000331 1163779e0 /
      data ak1 cs( 8) /    .0000000034 4597758e0 /
      data ak1 cs( 9) /   -.0000000003 8989323e0 /
      data ak1 cs(10) /    .0000000000 4720819e0 /
      data ak1 cs(11) /   -.0000000000 0604783e0 /
      data ak1 cs(12) /    .0000000000 0081284e0 /
      data ak1 cs(13) /   -.0000000000 0011386e0 /
      data ak1 cs(14) /    .0000000000 0001654e0 /
      data ak1 cs(15) /   -.0000000000 0000248e0 /
      data ak1 cs(16) /    .0000000000 0000038e0 /
      data ak1 cs(17) /   -.0000000000 0000006e0 /
c
c series for ak12       on the interval  0.          to  1.25000d-01
c                                        with weighted error   2.58e-17
c                                         log weighted error  16.59
c                               significant figures required  15.22
c                                    decimal places required  17.16
c
      data ak12cs( 1) /    .0637930834 3739001e0 /
      data ak12cs( 2) /    .0283288781 3049721e0 /
      data ak12cs( 3) /   -.0002475370 6739052e0 /
      data ak12cs( 4) /    .0000057719 7245160e0 /
      data ak12cs( 5) /   -.0000002068 9392195e0 /
      data ak12cs( 6) /    .0000000097 3998344e0 /
      data ak12cs( 7) /   -.0000000005 5853361e0 /
      data ak12cs( 8) /    .0000000000 3732996e0 /
      data ak12cs( 9) /   -.0000000000 0282505e0 /
      data ak12cs(10) /    .0000000000 0023720e0 /
      data ak12cs(11) /   -.0000000000 0002176e0 /
      data ak12cs(12) /    .0000000000 0000215e0 /
      data ak12cs(13) /   -.0000000000 0000022e0 /
      data ak12cs(14) /    .0000000000 0000002e0 /
c
      data ntk1, ntak1, ntak12, xmin, xsml / 3*0, 2*0. /
c
      if (ntk1.ne.0) go to 10
      ntk1 = inits (bk1cs, 11, 0.1*r1mach(3))
      ntak1 = inits (ak1cs, 17, 0.1*r1mach(3))
      ntak12 = inits (ak12cs, 14, 0.1*r1mach(3))
c
      xmin = exp (amax1(alog(r1mach(1)), -alog(r1mach(2))) + .01)
      xsml = sqrt (4.0*r1mach(3))
c
 10   if (x.le.0.) call seteru (29hbesk1e  x is zero or negative, 29,
     1  1, 2)
      if (x.gt.2.0) go to 20
c
      if (x.lt.xmin) call seteru (31hbesk1e  x so small k1 overflows,
     1  31, 2, 2)
      y = 0.
      if (x.gt.xsml) y = x*x
      besk1e = exp(x) * (alog(0.5*x)*besi1(x) +
     1  (0.75 + csevl (.5*y-1., bk1cs, ntk1))/x )
      return
c
 20   if (x.le.8.) besk1e = (1.25 + csevl ((16./x-5.)/3., ak1cs, ntak1))
     1  / sqrt(x)
      if (x.gt.8.) besk1e = (1.25 + csevl (16./x-1., ak12cs, ntak12))
     1  / sqrt(x)
c
      return
      end
      integer function i8save(isw,ivalue,set)
c
c  if (isw = 1) i8save returns the current error number and
c               sets it to ivalue if set = .true. .
c
c  if (isw = 2) i8save returns the current recovery switch and
c               sets it to ivalue if set = .true. .
c
      logical set
c
      integer iparam(2)
c  iparam(1) is the error number and iparam(2) is the recovery switch.
c
c  start execution error free and with recovery turned off.
c
      data iparam(1) /0/,  iparam(2) /2/
c
      i8save=iparam(isw)
      if (set) iparam(isw)=ivalue
c
      return
c
      end
      subroutine eprint
c
c  this subroutine prints the last error message, if any.
c
      integer messg(1)
c
      call e9rint(messg,1,1,.false.)
      return
c
      end
      subroutine s88fmt( n, w, ifmt )
c
c  s88fmt  replaces ifmt(1), ... , ifmt(n) with
c  the characters corresponding to the n least significant
c  digits of w.
c
      integer n,w,ifmt(n)
c
      integer nt,wt,digits(10)
c
      data digits( 1) / 1h0 /
      data digits( 2) / 1h1 /
      data digits( 3) / 1h2 /
      data digits( 4) / 1h3 /
      data digits( 5) / 1h4 /
      data digits( 6) / 1h5 /
      data digits( 7) / 1h6 /
      data digits( 8) / 1h7 /
      data digits( 9) / 1h8 /
      data digits(10) / 1h9 /
c
      nt = n
      wt = w
c
 10   if (nt .le. 0) return
        idigit = mod( wt, 10 )
        ifmt(nt) = digits(idigit+1)
        wt = wt/10
        nt = nt - 1
        go to 10
c
      end
      subroutine seterr (messg, nmessg, nerr, iopt)
c
c  this version modified by w. fullerton to dump if iopt = 1 and
c  not recovering.
c  seterr sets lerror = nerr, optionally prints the message and dumps
c  according to the following rules...
c
c    if iopt = 1 and recovering      - just remember the error.
c    if iopt = 1 and not recovering  - print, dump and stop.
c    if iopt = 2                     - print, dump and stop.
c
c  input
c
c    messg  - the error message.
c    nmessg - the length of the message, in characters.
c    nerr   - the error number. must have nerr non-zero.
c    iopt   - the option. must have iopt=1 or 2.
c
c  error states -
c
c    1 - message length not positive.
c    2 - cannot have nerr=0.
c    3 - an unrecovered error followed by another error.
c    4 - bad value for iopt.
c
c  only the first 72 characters of the message are printed.
c
c  the error handler calls a subroutine named fdump to produce a
c  symbolic dump. to complete the package, a dummy version of fdump
c  is supplied, but it should be replaced by a locally written version
c  which at least gives a trace-back.
c
      integer messg(1)
      external i1mach, i8save
c
c  the unit for error messages.
c
      iwunit=i1mach(4)
c
      if (nmessg.ge.1) go to 10
c
c  a message of non-positive length is fatal.
c
        write(iwunit,9000)
 9000   format(52h1error    1 in seterr - message length not positive.)
        go to 60
c
c  nw is the number of words the message occupies.
c
 10   nw=(min0(nmessg,72)-1)/i1mach(6)+1
c
      if (nerr.ne.0) go to 20
c
c  cannot turn the error state off using seterr.
c
        write(iwunit,9001)
 9001   format(42h1error    2 in seterr - cannot have nerr=0//
     1         34h the current error message follows///)
        call e9rint(messg,nw,nerr,.true.)
        itemp=i8save(1,1,.true.)
        go to 50
c
c  set lerror and test for a previous unrecovered error.
c
 20   if (i8save(1,nerr,.true.).eq.0) go to 30
c
        write(iwunit,9002)
 9002   format(23h1error    3 in seterr -,
     1         48h an unrecovered error followed by another error.//
     2         48h the previous and current error messages follow.///)
        call eprint
        call e9rint(messg,nw,nerr,.true.)
        go to 50
c
c  save this message in case it is not recovered from properly.
c
 30   call e9rint(messg,nw,nerr,.true.)
c
      if (iopt.eq.1 .or. iopt.eq.2) go to 40
c
c  must have iopt = 1 or 2.
c
        write(iwunit,9003)
 9003   format(42h1error    4 in seterr - bad value for iopt//
     1         34h the current error message follows///)
        go to 50
c
c  test for recovery.
c
 40   if (iopt.eq.2) go to 50
c
      if (i8save(2,0,.false.).eq.1) return
c
c     call eprint
c     stop
c
 50   call eprint
 60   call fdump
      stop
c
      end
      subroutine e9rint(messg,nw,nerr,save)
c
c  this routine stores the current error message or prints the old one,
c  if any, depending on whether or not save = .true. .
c
      integer messg(nw)
      logical save
      external i1mach, i8save
c
c  messgp stores at least the first 72 characters of the previous
c  message. its length is machine dependent and must be at least
c
c       1 + 71/(the number of characters stored per integer word).
c
      integer messgp(36),fmt(14),ccplus
c
c  start with no previous message.
c
      data messgp(1)/1h1/, nwp/0/, nerrp/0/
c
c  set up the format for printing the error message.
c  the format is simply (a1,14x,72axx) where xx=i1mach(6) is the
c  number of characters stored per integer word.
c
      data ccplus  / 1h+ /
c
      data fmt( 1) / 1h( /
      data fmt( 2) / 1ha /
      data fmt( 3) / 1h1 /
      data fmt( 4) / 1h, /
      data fmt( 5) / 1h1 /
      data fmt( 6) / 1h4 /
      data fmt( 7) / 1hx /
      data fmt( 8) / 1h, /
      data fmt( 9) / 1h7 /
      data fmt(10) / 1h2 /
      data fmt(11) / 1ha /
      data fmt(12) / 1hx /
      data fmt(13) / 1hx /
      data fmt(14) / 1h) /
c
      if (.not.save) go to 20
c
c  save the message.
c
        nwp=nw
        nerrp=nerr
        do 10 i=1,nw
 10     messgp(i)=messg(i)
c
        go to 30
c
 20   if (i8save(1,0,.false.).eq.0) go to 30
c
c  print the message.
c
        iwunit=i1mach(4)
        write(iwunit,9000) nerrp
 9000   format(7h error ,i4,4h in )
c
        call s88fmt(2,i1mach(6),fmt(12))
        write(iwunit,fmt) ccplus,(messgp(i),i=1,nwp)
c
 30   return
c
      end

      subroutine fdump
      return
      end
      REAL FUNCTION SECOND ()
C
C  OBTAIN ELAPSED CPU TIME SINCE START OF RUN.
C
C  ON CRAY -- THIS ROUTINE NOT NEEDED.
C
C  ON CONVEX -- MORE ACCURATE TIMINGS ARE OBTAINED BY CALLING
C               THE VECLIB ROUTINE CPUTIME AS FOLLOWS
C
C     SECOND = CPUTIME(0.0E0)
C
C  ON UNIX -- (SUN, FOR EXAMPLE) THE FOLLOWING MAY BE USED
C
      REAL TARRAY(2), ETIME, TOTAL
      TOTAL = ETIME(TARRAY)
      SECOND = TARRAY(1) + TARRAY(2)
C
      RETURN
      END