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=1000 )
C
      DOUBLE PRECISION   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  DVI0, DBESI0, DVI1, DBESI1, DVJ0, DBESJ0, DVJ1, DBESJ1,
     +          DVK0, DBESK0, DVK1, DBESK1, DVY0, DBESY0, DVY1, DBESY1
C
      DATA ((XLO(I,J),XHI(I,J),I=1,3),J=1,8) /
     +     0.0D0, 3.0D0 , -8.0D0, -3.0D0, 10.0D0, 50.0D0,
     +     0.0D0, 3.0D0 , -8.0D0, -3.0D0, 10.0D0, 50.0D0,
     +     0.0D0, 4.0D0 , -8.0D0, -4.0D0,  8.0D0, 50.0D0,
     +     0.0D0, 4.0D0 , -8.0D0, -4.0D0,  8.0D0, 50.0D0,
     +     0.0D0, 2.0D0 ,  2.0D0,  8.0D0, 10.0D0, 50.0D0,
     +     0.0D0, 2.0D0 ,  2.0D0,  8.0D0, 10.0D0, 50.0D0,
     +     0.0D0, 1.0D-3,  1.0D0,  4.0D0,  4.0D0, 50.0D0,
     +     0.0D0, 1.0D-3,  1.0D0,  4.0D0,  4.0D0, 50.0D0  /
C
C----------------------------------------------------------------------
C
C     ... TEST EACH USER-CALLABLE ROUTINE
C
      CALL TEST('I0',DBESI0,DVI0,XLO(1,1),XHI(1,1),M,X,F,FX,WORK,IWORK)
      CALL TEST('I1',DBESI1,DVI1,XLO(1,2),XHI(1,2),M,X,F,FX,WORK,IWORK)
      CALL TEST('J0',DBESJ0,DVJ0,XLO(1,3),XHI(1,3),M,X,F,FX,WORK,IWORK)
      CALL TEST('J1',DBESJ1,DVJ1,XLO(1,4),XHI(1,4),M,X,F,FX,WORK,IWORK)
      CALL TEST('K0',DBESK0,DVK0,XLO(1,5),XHI(1,5),M,X,F,FX,WORK,IWORK)
      CALL TEST('K1',DBESK1,DVK1,XLO(1,6),XHI(1,6),M,X,F,FX,WORK,IWORK)
      CALL TEST('Y0',DBESY0,DVY0,XLO(1,7),XHI(1,7),M,X,F,FX,WORK,IWORK)
      CALL TEST('Y1',DBESY1,DVY1,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
      DOUBLE PRECISION 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=5, NTRY=19 )
C
      INTEGER I, INFO, IREP, IOFF, IPC(NTRY), JPC(NTRY), K, KPC,
     +        MA, MB, MC, NFAIL, NREPV
      DOUBLE PRECISION DIFF, DMAX, DX, EPS, D1MACH
      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.0D0*D1MACH(4)

C     ... TRIGGER INITIALIZATIONS
C
      X(1) = 1.0D0
      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 differences .gt. eps'
      PRINT *,' Dmax    = maximum difference'
      PRINT *
      PRINT *,' Times given in seconds/argument.'
      PRINT *
      PRINT *,' Differences are abs(s-v)/max(abs(s),1))'
      PRINT *,' where s = FNLIB result and v = VFNLIB result.'
      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.0D0
         NFAIL = 0
         DO 70 I=1,M
            DIFF = ABS(FX(I)-F(I))/MAX(ABS(FX(I)),1.0D0)
            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
      double precision function dbesi0 (x)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, bi0cs(18), xmax, xsml, y, d1mach,
     1  dcsevl, dbsi0e, dexp, dlog, dsqrt
c     external d1mach, dbsi0e, dcsevl, dexp, dlog, dsqrt, initds
c
c series for bi0        on the interval  0.          to  9.00000e+00
c                                        with weighted error   9.51e-34
c                                         log weighted error  33.02
c                               significant figures required  33.31
c                                    decimal places required  33.65
c
      data bi0 cs(  1) / -.7660547252 8391449510 8189497624 3285 d-1   /
      data bi0 cs(  2) / +.1927337953 9938082699 5240875088 1196 d+1   /
      data bi0 cs(  3) / +.2282644586 9203013389 3702929233 0415 d+0   /
      data bi0 cs(  4) / +.1304891466 7072904280 7933421069 1888 d-1   /
      data bi0 cs(  5) / +.4344270900 8164874513 7868268102 6107 d-3   /
      data bi0 cs(  6) / +.9422657686 0019346639 2317174411 8766 d-5   /
      data bi0 cs(  7) / +.1434006289 5106910799 6209187817 9957 d-6   /
      data bi0 cs(  8) / +.1613849069 6617490699 1541971999 4611 d-8   /
      data bi0 cs(  9) / +.1396650044 5356696994 9509270814 2522 d-10  /
      data bi0 cs( 10) / +.9579451725 5054453446 2752317189 3333 d-13  /
      data bi0 cs( 11) / +.5333981859 8625021310 1510774400 0000 d-15  /
      data bi0 cs( 12) / +.2458716088 4374707746 9678591999 9999 d-17  /
      data bi0 cs( 13) / +.9535680890 2487700269 4434133333 3333 d-20  /
      data bi0 cs( 14) / +.3154382039 7214273367 8933333333 3333 d-22  /
      data bi0 cs( 15) / +.9004564101 0946374314 6666666666 6666 d-25  /
      data bi0 cs( 16) / +.2240647369 1236700160 0000000000 0000 d-27  /
      data bi0 cs( 17) / +.4903034603 2428373333 3333333333 3333 d-30  /
      data bi0 cs( 18) / +.9508172606 1226666666 6666666666 6666 d-33  /
c
      data nti0, xsml, xmax / 0, 2*0.d0 /
c
      if (nti0.ne.0) go to 10
      nti0 = initds (bi0cs, 18, 0.1*sngl(d1mach(3)))
      xsml = dsqrt (8.0d0*d1mach(3))
      xmax = dlog (d1mach(2))
c
 10   y = dabs(x)
      if (y.gt.3.0d0) go to 20
c
      dbesi0 = 1.0d0
      if (y.gt.xsml) dbesi0 = 2.75d0 + dcsevl (y*y/4.5d0-1.d0, bi0cs,
     1  nti0)
      return
c
 20   if (y.gt.xmax) call seteru (
     1  35hdbesi0  dabs(x) so big i0 overflows, 35, 2, 2)
c
      dbesi0 = dexp(y) * dbsi0e(x)
c
      return
      end
      double precision function dbesi1 (x)
c oct 1983 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, bi1cs(17), xmax, xmin, xsml, y, d1mach,
     1  dcsevl, dbsi1e, dexp, dlog, dsqrt
c     external d1mach, dbsi1e, dcsevl, dexp, dlog, dsqrt, initds
c
c series for bi1        on the interval  0.          to  9.00000e+00
c                                        with weighted error   1.44e-32
c                                         log weighted error  31.84
c                               significant figures required  31.45
c                                    decimal places required  32.46
c
      data bi1 cs(  1) / -.1971713261 0998597316 1385032181 49 d-2     /
      data bi1 cs(  2) / +.4073488766 7546480608 1553936520 14 d+0     /
      data bi1 cs(  3) / +.3483899429 9959455866 2450377837 87 d-1     /
      data bi1 cs(  4) / +.1545394556 3001236038 5984010584 89 d-2     /
      data bi1 cs(  5) / +.4188852109 8377784129 4588320041 20 d-4     /
      data bi1 cs(  6) / +.7649026764 8362114741 9597039660 69 d-6     /
      data bi1 cs(  7) / +.1004249392 4741178689 1798080372 38 d-7     /
      data bi1 cs(  8) / +.9932207791 9238106481 3712980548 63 d-10    /
      data bi1 cs(  9) / +.7663801791 8447637275 2001716813 49 d-12    /
      data bi1 cs( 10) / +.4741418923 8167394980 3880919481 60 d-14    /
      data bi1 cs( 11) / +.2404114404 0745181799 8631720320 00 d-16    /
      data bi1 cs( 12) / +.1017150500 7093713649 1211007999 99 d-18    /
      data bi1 cs( 13) / +.3645093565 7866949458 4917333333 33 d-21    /
      data bi1 cs( 14) / +.1120574950 2562039344 8106666666 66 d-23    /
      data bi1 cs( 15) / +.2987544193 4468088832 0000000000 00 d-26    /
      data bi1 cs( 16) / +.6973231093 9194709333 3333333333 33 d-29    /
      data bi1 cs( 17) / +.1436794822 0620800000 0000000000 00 d-31    /
c
      data nti1, xmin, xsml, xmax / 0, 3*0.d0 /
c
      if (nti1.ne.0) go to 10
      nti1 = initds (bi1cs, 17, 0.1*sngl(d1mach(3)))
      xmin = 2.0d0*d1mach(1)
      xsml = dsqrt (8.0d0*d1mach(3))
      xmax = dlog (d1mach(2))
c
 10   y = dabs(x)
      if (y.gt.3.0d0) go to 20
c
      dbesi1 = 0.0d0
      if (y.eq.0.0d0) return
c
      if (y.le.xmin) call seteru (
     1  38hdbesi1  dabs(x) so small i1 underflows, 38, 1, 0)
      if (y.gt.xmin) dbesi1 = 0.5d0*x
      if (y.gt.xsml) dbesi1 = x*(0.875d0 + dcsevl (y*y/4.5d0-1.d0,
     1  bi1cs, nti1))
      return
c
 20   if (y.gt.xmax) call seteru (
     1  35hdbesi1  dabs(x) so big i1 overflows, 35, 2, 2)
c
      dbesi1 = dexp(y) * dbsi1e(x)
c
      return
      end
      double precision function dbesj0 (x)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, bj0cs(19), ampl, theta, xsml, y,
     1  d1mach, dcsevl, dcos, dsqrt
c     external d1mach, dcos, dcsevl, dsqrt, initds
c
c series for bj0        on the interval  0.          to  1.60000e+01
c                                        with weighted error   4.39e-32
c                                         log weighted error  31.36
c                               significant figures required  31.21
c                                    decimal places required  32.00
c
      data bj0 cs(  1) / +.1002541619 6893913701 0731272640 74 d+0     /
      data bj0 cs(  2) / -.6652230077 6440513177 6787578311 24 d+0     /
      data bj0 cs(  3) / +.2489837034 9828131370 4604687266 80 d+0     /
      data bj0 cs(  4) / -.3325272317 0035769653 8843415038 54 d-1     /
      data bj0 cs(  5) / +.2311417930 4694015462 9049241177 29 d-2     /
      data bj0 cs(  6) / -.9911277419 9508092339 0485193365 49 d-4     /
      data bj0 cs(  7) / +.2891670864 3998808884 7339037470 78 d-5     /
      data bj0 cs(  8) / -.6121085866 3032635057 8184074815 16 d-7     /
      data bj0 cs(  9) / +.9838650793 8567841324 7687486364 15 d-9     /
      data bj0 cs( 10) / -.1242355159 7301765145 5158970068 36 d-10    /
      data bj0 cs( 11) / +.1265433630 2559045797 9158272103 63 d-12    /
      data bj0 cs( 12) / -.1061945649 5287244546 9148175129 59 d-14    /
      data bj0 cs( 13) / +.7470621075 8024567437 0989155840 00 d-17    /
      data bj0 cs( 14) / -.4469703227 4412780547 6270079999 99 d-19    /
      data bj0 cs( 15) / +.2302428158 4337436200 5230933333 33 d-21    /
      data bj0 cs( 16) / -.1031914479 4166698148 5226666666 66 d-23    /
      data bj0 cs( 17) / +.4060817827 4873322700 8000000000 00 d-26    /
      data bj0 cs( 18) / -.1414383600 5240913919 9999999999 99 d-28    /
      data bj0 cs( 19) / +.4391090549 6698880000 0000000000 00 d-31    /
c
      data ntj0, xsml / 0, 0.d0 /
c
      if (ntj0.ne.0) go to 10
      ntj0 = initds (bj0cs, 19, 0.1*sngl(d1mach(3)))
      xsml = dsqrt (4.0d0*d1mach(3))
c
 10   y = dabs(x)
      if (y.gt.4.0d0) go to 20
c
      dbesj0 = 1.0d0
      if (y.gt.xsml) dbesj0 = dcsevl (.125d0*y*y-1.d0, bj0cs, ntj0)
      return
c
 20   call d9b0mp (y, ampl, theta)
      dbesj0 = ampl * dcos(theta)
c
      return
      end
      double precision function dbesj1 (x)
c sept 1983 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, bj1cs(19), ampl, theta, xsml, xmin, y,
     1  d1mach, dcsevl, dcos, dsqrt
c     external d1mach, dcos, dcsevl, dsqrt, initds
c
c series for bj1        on the interval  0.          to  1.60000e+01
c                                        with weighted error   1.16e-33
c                                         log weighted error  32.93
c                               significant figures required  32.36
c                                    decimal places required  33.57
c
      data bj1 cs(  1) / -.1172614151 3332786560 6240574524 003 d+0    /
      data bj1 cs(  2) / -.2536152183 0790639562 3030884554 698 d+0    /
      data bj1 cs(  3) / +.5012708098 4469568505 3656363203 743 d-1    /
      data bj1 cs(  4) / -.4631514809 6250819184 2619728789 772 d-2    /
      data bj1 cs(  5) / +.2479962294 1591402453 9124064592 364 d-3    /
      data bj1 cs(  6) / -.8678948686 2788258452 1246435176 416 d-5    /
      data bj1 cs(  7) / +.2142939171 4379369150 2766250991 292 d-6    /
      data bj1 cs(  8) / -.3936093079 1831797922 9322764073 061 d-8    /
      data bj1 cs(  9) / +.5591182317 9468800401 8248059864 032 d-10   /
      data bj1 cs( 10) / -.6327616404 6613930247 7695274014 880 d-12   /
      data bj1 cs( 11) / +.5840991610 8572470032 6945563268 266 d-14   /
      data bj1 cs( 12) / -.4482533818 7012581903 9135059199 999 d-16   /
      data bj1 cs( 13) / +.2905384492 6250246630 6018688000 000 d-18   /
      data bj1 cs( 14) / -.1611732197 8414416541 2118186666 666 d-20   /
      data bj1 cs( 15) / +.7739478819 3927463729 8346666666 666 d-23   /
      data bj1 cs( 16) / -.3248693782 1119984114 3466666666 666 d-25   /
      data bj1 cs( 17) / +.1202237677 2274102272 0000000000 000 d-27   /
      data bj1 cs( 18) / -.3952012212 6513493333 3333333333 333 d-30   /
      data bj1 cs( 19) / +.1161678082 2664533333 3333333333 333 d-32   /
c
      data ntj1, xsml, xmin / 0, 2*0.d0 /
c
      if (ntj1.ne.0) go to 10
      ntj1 = initds (bj1cs, 19, 0.1*sngl(d1mach(3)))
c
      xsml = dsqrt (4.0d0*d1mach(3))
      xmin = 2.0d0*d1mach(1)
c
 10   y = dabs(x)
      if (y.gt.4.0d0) go to 20
c
      dbesj1 = 0.0d0
      if (y.eq.0.0d0) return
c
      if (y.le.xmin) call seteru (
     1  38hdbesj1  dabs(x) so small j1 underflows, 38, 1, 0)
      if (y.gt.xmin) dbesj1 = 0.5d0*x
      if (y.gt.xsml) dbesj1 = x*(.25d0 + dcsevl (.125d0*y*y-1.d0,
     1  bj1cs, ntj1) )
      return
c
 20   call d9b1mp (y, ampl, theta)
C**Next line changed from*****************
C     dbesj1 = ampl * dcos(theta)
C  to
      dbesj1 = sign(ampl,x) * dcos(theta)
C**by Ron Boisvert, 19 Mar 91 ************
c
      return
      end
      double precision function dbesk0 (x)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, bk0cs(16), xmax, xsml, y,
     1  d1mach, dcsevl, dbesi0, dbsk0e, dexp, dlog, dsqrt
c     external d1mach, dbesi0, dbsk0e, dcsevl, dexp, dlog, dsqrt, initds
c
c series for bk0        on the interval  0.          to  4.00000e+00
c                                        with weighted error   3.08e-33
c                                         log weighted error  32.51
c                               significant figures required  32.05
c                                    decimal places required  33.11
c
      data bk0 cs(  1) / -.3532739323 3902768720 1140060063 153 d-1    /
      data bk0 cs(  2) / +.3442898999 2462848688 6344927529 213 d+0    /
      data bk0 cs(  3) / +.3597993651 5361501626 5721303687 231 d-1    /
      data bk0 cs(  4) / +.1264615411 4469259233 8479508673 447 d-2    /
      data bk0 cs(  5) / +.2286212103 1194517860 8269830297 585 d-4    /
      data bk0 cs(  6) / +.2534791079 0261494573 0790013428 354 d-6    /
      data bk0 cs(  7) / +.1904516377 2202088589 7214059381 366 d-8    /
      data bk0 cs(  8) / +.1034969525 7633624585 1008317853 089 d-10   /
      data bk0 cs(  9) / +.4259816142 7910825765 2445327170 133 d-13   /
      data bk0 cs( 10) / +.1374465435 8807508969 4238325440 000 d-15   /
      data bk0 cs( 11) / +.3570896528 5083735909 9688597333 333 d-18   /
      data bk0 cs( 12) / +.7631643660 1164373766 7498666666 666 d-21   /
      data bk0 cs( 13) / +.1365424988 4407818590 8053333333 333 d-23   /
      data bk0 cs( 14) / +.2075275266 9066680831 9999999999 999 d-26   /
      data bk0 cs( 15) / +.2712814218 0729856000 0000000000 000 d-29   /
      data bk0 cs( 16) / +.3082593887 9146666666 6666666666 666 d-32   /
c
      data ntk0, xsml, xmax / 0, 2*0.d0 /
c
      if (ntk0.ne.0) go to 10
      ntk0 = initds (bk0cs, 16, 0.1*sngl(d1mach(3)))
      xsml = dsqrt (4.0d0*d1mach(3))
      xmax = -dlog(d1mach(1))
      xmax = xmax - 0.5d0*xmax*dlog(xmax)/(xmax+0.5d0)
c
 10   if (x.le.0.d0) call seteru (29hdbesk0  x is zero or negative,
     1  29, 2, 2)
      if (x.gt.2.0d0) go to 20
c
      y = 0.d0
      if (x.gt.xsml) y = x*x
      dbesk0 = -dlog(0.5d0*x)*dbesi0(x) - 0.25d0 + dcsevl (.5d0*y-1.d0,
     1  bk0cs, ntk0)
      return
c
 20   dbesk0 = 0.d0
      if (x.gt.xmax) call seteru (30hdbesk0  x so big k0 underflows,
     1  30, 1, 0)
      if (x.gt.xmax) return
c
      dbesk0 = dexp(-x) * dbsk0e(x)
c
      return
      end
      double precision function dbesk1 (x)
c july 1980 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, bk1cs(16), xmax, xmin, xsml, y,
     1  d1mach, dcsevl, dbesi1, dbsk1e, dexp, dlog, dsqrt
c     external d1mach, dbesi1, dbsk1e, dcsevl, dexp, dlog, dsqrt, initds
c
c series for bk1        on the interval  0.          to  4.00000e+00
c                                        with weighted error   9.16e-32
c                                         log weighted error  31.04
c                               significant figures required  30.61
c                                    decimal places required  31.64
c
      data bk1 cs(  1) / +.2530022733 8947770532 5311208685 33 d-1     /
      data bk1 cs(  2) / -.3531559607 7654487566 7238316918 01 d+0     /
      data bk1 cs(  3) / -.1226111808 2265714823 4790679300 42 d+0     /
      data bk1 cs(  4) / -.6975723859 6398643501 8129202960 83 d-2     /
      data bk1 cs(  5) / -.1730288957 5130520630 1765073689 79 d-3     /
      data bk1 cs(  6) / -.2433406141 5659682349 6007350301 64 d-5     /
      data bk1 cs(  7) / -.2213387630 7347258558 3152525451 26 d-7     /
      data bk1 cs(  8) / -.1411488392 6335277610 9583302126 08 d-9     /
      data bk1 cs(  9) / -.6666901694 1993290060 8537512643 73 d-12    /
      data bk1 cs( 10) / -.2427449850 5193659339 2631968648 53 d-14    /
      data bk1 cs( 11) / -.7023863479 3862875971 7837971200 00 d-17    /
      data bk1 cs( 12) / -.1654327515 5100994675 4910293333 33 d-19    /
      data bk1 cs( 13) / -.3233834745 9944491991 8933333333 33 d-22    /
      data bk1 cs( 14) / -.5331275052 9265274999 4666666666 66 d-25    /
      data bk1 cs( 15) / -.7513040716 2157226666 6666666666 66 d-28    /
      data bk1 cs( 16) / -.9155085717 6541866666 6666666666 66 d-31    /
c
      data ntk1, xmin, xsml, xmax / 0, 3*0.d0 /
c
      if (ntk1.ne.0) go to 10
      ntk1 = initds (bk1cs, 16, 0.1*sngl(d1mach(3)))
      xmin = dexp (dmax1(dlog(d1mach(1)), -dlog(d1mach(2))) + 0.01d0)
      xsml = dsqrt (4.0d0*d1mach(3))
      xmax = -dlog(d1mach(1))
      xmax = xmax - 0.5d0*xmax*dlog(xmax)/(xmax+0.5d0) - 0.01d0
c
 10   if (x.le.0.d0) call seteru (29hdbesk1  x is zero or negative,
     1  29, 2, 2)
      if (x.gt.2.0d0) go to 20
c
      if (x.lt.xmin) call seteru (31hdbesk1  x so small k1 overflows,
     1  31, 3, 2)
      y = 0.d0
      if (x.gt.xsml) y = x*x
      dbesk1 = dlog(0.5d0*x)*dbesi1(x) + (0.75d0 + dcsevl (.5d0*y-1.d0,
     1  bk1cs, ntk1))/x
      return
c
 20   dbesk1 = 0.d0
      if (x.gt.xmax) call seteru (30hdbesk1  x so big k1 underflows,
     1  30, 1, 0)
      if (x.gt.xmax) return
c
      dbesk1 = dexp(-x) * dbsk1e(x)
c
      return
      end
      double precision function dbesy0 (x)
c august 1980 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, by0cs(19), ampl, theta, twodpi, xsml,
     1  y, alnhaf, d1mach, dcsevl, dbesj0, dlog, dsin, dsqrt
c     external d1mach, dbesj0, dcsevl, dlog, dsin, dsqrt, initds
c
c series for by0        on the interval  0.          to  1.60000e+01
c                                        with weighted error   8.14e-32
c                                         log weighted error  31.09
c                               significant figures required  30.31
c                                    decimal places required  31.73
c
      data by0 cs(  1) / -.1127783939 2865573217 9398054602 8 d-1      /
      data by0 cs(  2) / -.1283452375 6042034604 8088453183 8 d+0      /
      data by0 cs(  3) / -.1043788479 9794249365 8176227661 8 d+0      /
      data by0 cs(  4) / +.2366274918 3969695409 2415926461 3 d-1      /
      data by0 cs(  5) / -.2090391647 7004862391 9622395034 2 d-2      /
      data by0 cs(  6) / +.1039754539 3905725209 9924657638 1 d-3      /
      data by0 cs(  7) / -.3369747162 4239720967 1877534503 7 d-5      /
      data by0 cs(  8) / +.7729384267 6706671585 2136721637 1 d-7      /
      data by0 cs(  9) / -.1324976772 6642595914 4347606896 4 d-8      /
      data by0 cs( 10) / +.1764823261 5404527921 0038936315 8 d-10     /
      data by0 cs( 11) / -.1881055071 5801962006 0282301206 9 d-12     /
      data by0 cs( 12) / +.1641865485 3661495027 9223718574 9 d-14     /
      data by0 cs( 13) / -.1195659438 6046060857 4599100672 0 d-16     /
      data by0 cs( 14) / +.7377296297 4401858424 9411242666 6 d-19     /
      data by0 cs( 15) / -.3906843476 7104373307 4090666666 6 d-21     /
      data by0 cs( 16) / +.1795503664 4361579498 2912000000 0 d-23     /
      data by0 cs( 17) / -.7229627125 4480104789 3333333333 3 d-26     /
      data by0 cs( 18) / +.2571727931 6351685973 3333333333 3 d-28     /
      data by0 cs( 19) / -.8141268814 1636949333 3333333333 3 d-31     /
c
      data twodpi / 0.6366197723 6758134307 5535053490 057 d0 /
      data alnhaf /-0.6931471805 5994530941 7232121458 18d0 /
      data nty0, xsml / 0, 0.d0 /
c
      if (nty0.ne.0) go to 10
      nty0 = initds (by0cs, 19, 0.1*sngl(d1mach(3)))
      xsml = dsqrt (4.0d0*d1mach(3))
c
 10   if (x.le.0.d0) call seteru (29hdbesy0  x is zero or negative, 29,
     1  1, 2)
      if (x.gt.4.0d0) go to 20
c
      y = 0.d0
      if (x.gt.xsml) y = x*x
      dbesy0 = twodpi*(alnhaf+dlog(x))*dbesj0(x) + .375d0
     1  + dcsevl (.125d0*y-1.d0, by0cs, nty0)
      return
c
 20   call d9b0mp (x, ampl, theta)
      dbesy0 = ampl * dsin(theta)
      return
c
      end
      double precision function dbesy1 (x)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, by1cs(20), ampl, theta, twodpi, xmin, xsml,
     1  y, d1mach, dcsevl, dbesj1, dexp, dlog, dsin, dsqrt
c     external d1mach, dbesj1, dcsevl, dexp, dlog, dsin, dsqrt, initds
c
c series for by1        on the interval  0.          to  1.60000e+01
c                                        with weighted error   8.65e-33
c                                         log weighted error  32.06
c                               significant figures required  32.17
c                                    decimal places required  32.71
c
      data by1 cs(  1) / +.3208047100 6119086293 2352018628 015 d-1    /
      data by1 cs(  2) / +.1262707897 4335004495 3431725999 727 d+1    /
      data by1 cs(  3) / +.6499961899 9231750009 7490637314 144 d-2    /
      data by1 cs(  4) / -.8936164528 8605041165 3144160009 712 d-1    /
      data by1 cs(  5) / +.1325088122 1757095451 2375510370 043 d-1    /
      data by1 cs(  6) / -.8979059119 6483523775 3039508298 105 d-3    /
      data by1 cs(  7) / +.3647361487 9583067824 2287368165 349 d-4    /
      data by1 cs(  8) / -.1001374381 6660005554 9075523845 295 d-5    /
      data by1 cs(  9) / +.1994539657 3901739703 1159372421 243 d-7    /
      data by1 cs( 10) / -.3023065601 8033816728 4799332520 743 d-9    /
      data by1 cs( 11) / +.3609878156 9478119611 6252914242 474 d-11   /
      data by1 cs( 12) / -.3487488297 2875824241 4552947409 066 d-13   /
      data by1 cs( 13) / +.2783878971 5591766581 3507698517 333 d-15   /
      data by1 cs( 14) / -.1867870968 6194876876 6825352533 333 d-17   /
      data by1 cs( 15) / +.1068531533 9116825975 7070336000 000 d-19   /
      data by1 cs( 16) / -.5274721956 6844822894 3872000000 000 d-22   /
      data by1 cs( 17) / +.2270199403 1556641437 0133333333 333 d-24   /
      data by1 cs( 18) / -.8595390353 9452310869 3333333333 333 d-27   /
      data by1 cs( 19) / +.2885404379 8337945600 0000000000 000 d-29   /
      data by1 cs( 20) / -.8647541138 9371733333 3333333333 333 d-32   /
c
      data twodpi / 0.6366197723 6758134307 5535053490 057 d0 /
      data nty1, xmin, xsml / 0, 2*0.d0 /
c
      if (nty1.ne.0) go to 10
      nty1 = initds (by1cs, 20, 0.1*sngl(d1mach(3)))
c
      xmin = 1.571d0 * dexp (dmax1(dlog(d1mach(1)), -dlog(d1mach(2))) +
     1  0.01d0)
      xsml = dsqrt (4.0d0*d1mach(3))
c
 10   if (x.le.0.d0) call seteru (29hdbesy1  x is zero or negative, 29,
     1  1, 2)
      if (x.gt.4.0d0) go to 20
c
      if (x.lt.xmin) call seteru (31hdbesy1  x so small y1 overflows,
     1  31, 3, 2)
      y = 0.d0
      if (x.gt.xsml) y = x*x
      dbesy1 = twodpi * dlog(0.5d0*x)*dbesj1(x) + (0.5d0 +
     1  dcsevl (.125d0*y-1.d0, by1cs, nty1))/x
      return
c
 20   call d9b1mp (x, ampl, theta)
      dbesy1 = ampl * dsin(theta)
      return
c
      end
      function initds (dos, nos, eta)
c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
c
c initialize the double precision orthogonal series dos so that initds
c is the number of terms needed to insure the error is no larger than
c eta.  ordinarily eta will be chosen to be one-tenth machine precision.
c
c             input arguments --
c dos    dble prec array of nos coefficients in an orthogonal series.
c nos    number of coefficients in dos.
c eta    requested accuracy of series.
c
      double precision dos(nos)
c
      if (nos.lt.1) call seteru (
     1  35hinitds  number of coefficients lt 1, 35, 2, 2)
c
      err = 0.
      do 10 ii=1,nos
        i = nos + 1 - ii
        err = err + abs(sngl(dos(i)))
        if (err.gt.eta) go to 20
 10   continue
c
 20   if (i.eq.nos) call seteru (28hinitds  eta may be too small, 28,
     1  1, 2)
      initds = i
c
      return
      end
      double precision function dcsevl (x, a, n)
c
c evaluate the n-term chebyshev series a at x.  adapted from
c r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973).
c
c             input arguments --
c x      dble prec value at which the series is to be evaluated.
c a      dble prec array of n terms of a chebyshev series.  in eval-
c        uating a, only half the first coef is summed.
c n      number of terms in array a.
c
      double precision a(n), x, twox, b0, b1, b2
c
      if (n.lt.1) call seteru (28hdcsevl  number of terms le 0, 28, 2,2)
      if (n.gt.1000) call seteru (31hdcsevl  number of terms gt 1000,
     1  31, 3, 2)
      if (x.lt.(-1.1d0) .or. x.gt.1.1d0) call seteru (
     1  25hdcsevl  x outside (-1,+1), 25, 1, 1)
c
      twox = 2.0d0*x
      b1 = 0.d0
      b0 = 0.d0
      do 10 i=1,n
        b2 = b1
        b1 = b0
        ni = n - i + 1
        b0 = twox*b1 - b2 + a(ni)
 10   continue
c
      dcsevl = 0.5d0 * (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
      double precision function dbsi0e (x)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, bi0cs(18), ai0cs(46), ai02cs(69),
     1  xsml, y, d1mach, dcsevl, dexp, dsqrt
c     external d1mach, dcsevl, dexp, dsqrt, initds
c
c series for bi0        on the interval  0.          to  9.00000e+00
c                                        with weighted error   9.51e-34
c                                         log weighted error  33.02
c                               significant figures required  33.31
c                                    decimal places required  33.65
c
      data bi0 cs(  1) / -.7660547252 8391449510 8189497624 3285 d-1   /
      data bi0 cs(  2) / +.1927337953 9938082699 5240875088 1196 d+1   /
      data bi0 cs(  3) / +.2282644586 9203013389 3702929233 0415 d+0   /
      data bi0 cs(  4) / +.1304891466 7072904280 7933421069 1888 d-1   /
      data bi0 cs(  5) / +.4344270900 8164874513 7868268102 6107 d-3   /
      data bi0 cs(  6) / +.9422657686 0019346639 2317174411 8766 d-5   /
      data bi0 cs(  7) / +.1434006289 5106910799 6209187817 9957 d-6   /
      data bi0 cs(  8) / +.1613849069 6617490699 1541971999 4611 d-8   /
      data bi0 cs(  9) / +.1396650044 5356696994 9509270814 2522 d-10  /
      data bi0 cs( 10) / +.9579451725 5054453446 2752317189 3333 d-13  /
      data bi0 cs( 11) / +.5333981859 8625021310 1510774400 0000 d-15  /
      data bi0 cs( 12) / +.2458716088 4374707746 9678591999 9999 d-17  /
      data bi0 cs( 13) / +.9535680890 2487700269 4434133333 3333 d-20  /
      data bi0 cs( 14) / +.3154382039 7214273367 8933333333 3333 d-22  /
      data bi0 cs( 15) / +.9004564101 0946374314 6666666666 6666 d-25  /
      data bi0 cs( 16) / +.2240647369 1236700160 0000000000 0000 d-27  /
      data bi0 cs( 17) / +.4903034603 2428373333 3333333333 3333 d-30  /
      data bi0 cs( 18) / +.9508172606 1226666666 6666666666 6666 d-33  /
c
c series for ai0        on the interval  1.25000e-01 to  3.33333e-01
c                                        with weighted error   2.74e-32
c                                         log weighted error  31.56
c                               significant figures required  30.15
c                                    decimal places required  32.39
c
      data ai0 cs(  1) / +.7575994494 0237959427 2987203743 8 d-1      /
      data ai0 cs(  2) / +.7591380810 8233455072 9297873320 4 d-2      /
      data ai0 cs(  3) / +.4153131338 9237505018 6319749138 2 d-3      /
      data ai0 cs(  4) / +.1070076463 4390730735 8242970217 0 d-4      /
      data ai0 cs(  5) / -.7901179979 2128946607 5031948573 0 d-5      /
      data ai0 cs(  6) / -.7826143501 4387522697 8898980690 9 d-6      /
      data ai0 cs(  7) / +.2783849942 9488708063 8118538985 7 d-6      /
      data ai0 cs(  8) / +.8252472600 6120271919 6682913319 8 d-8      /
      data ai0 cs(  9) / -.1204463945 5201991790 5496089110 3 d-7      /
      data ai0 cs( 10) / +.1559648598 5060764436 1228752792 8 d-8      /
      data ai0 cs( 11) / +.2292556367 1033165434 7725480285 7 d-9      /
      data ai0 cs( 12) / -.1191622884 2790646036 7777423447 8 d-9      /
      data ai0 cs( 13) / +.1757854916 0324098302 1833124774 3 d-10     /
      data ai0 cs( 14) / +.1128224463 2189005171 4441135682 4 d-11     /
      data ai0 cs( 15) / -.1146848625 9272988777 2963387698 2 d-11     /
      data ai0 cs( 16) / +.2715592054 8036628726 4365192160 6 d-12     /
      data ai0 cs( 17) / -.2415874666 5626878384 4247572028 1 d-13     /
      data ai0 cs( 18) / -.6084469888 2551250646 0609963922 4 d-14     /
      data ai0 cs( 19) / +.3145705077 1754772937 0836026730 3 d-14     /
      data ai0 cs( 20) / -.7172212924 8711877179 6217505917 6 d-15     /
      data ai0 cs( 21) / +.7874493403 4541033960 8390960332 7 d-16     /
      data ai0 cs( 22) / +.1004802753 0094624023 4524457183 9 d-16     /
      data ai0 cs( 23) / -.7566895365 3505348534 2843588881 0 d-17     /
      data ai0 cs( 24) / +.2150380106 8761198878 1205128784 5 d-17     /
      data ai0 cs( 25) / -.3754858341 8308744291 5158445260 8 d-18     /
      data ai0 cs( 26) / +.2354065842 2269925769 0075710532 2 d-19     /
      data ai0 cs( 27) / +.1114667612 0479285302 2637335511 0 d-19     /
      data ai0 cs( 28) / -.5398891884 3969903786 9677932270 9 d-20     /
      data ai0 cs( 29) / +.1439598792 2407526770 4285840452 2 d-20     /
      data ai0 cs( 30) / -.2591916360 1110934064 6081840196 2 d-21     /
      data ai0 cs( 31) / +.2238133183 9985839074 3409229824 0 d-22     /
      data ai0 cs( 32) / +.5250672575 3647711727 7221683199 9 d-23     /
      data ai0 cs( 33) / -.3249904138 5332307841 7343228586 6 d-23     /
      data ai0 cs( 34) / +.9924214103 2050379278 5728471040 0 d-24     /
      data ai0 cs( 35) / -.2164992254 2446695231 4655429973 3 d-24     /
      data ai0 cs( 36) / +.3233609471 9435940839 7333299199 9 d-25     /
      data ai0 cs( 37) / -.1184620207 3967424898 2473386666 6 d-26     /
      data ai0 cs( 38) / -.1281671853 9504986505 4833868799 9 d-26     /
      data ai0 cs( 39) / +.5827015182 2793905116 0556885333 3 d-27     /
      data ai0 cs( 40) / -.1668222326 0261097193 6450150399 9 d-27     /
      data ai0 cs( 41) / +.3625309510 5415699757 0068480000 0 d-28     /
      data ai0 cs( 42) / -.5733627999 0557135899 4595839999 9 d-29     /
      data ai0 cs( 43) / +.3736796722 0630982296 4258133333 3 d-30     /
      data ai0 cs( 44) / +.1602073983 1568519633 6551253333 3 d-30     /
      data ai0 cs( 45) / -.8700424864 0572298845 2249599999 9 d-31     /
      data ai0 cs( 46) / +.2741320937 9374811456 0341333333 3 d-31     /
c
c series for ai02       on the interval  0.          to  1.25000e-01
c                                        with weighted error   1.97e-32
c                                         log weighted error  31.71
c                               significant figures required  30.15
c                                    decimal places required  32.63
c
      data ai02cs(  1) / +.5449041101 4108831607 8960962268 0 d-1      /
      data ai02cs(  2) / +.3369116478 2556940898 9785662979 9 d-2      /
      data ai02cs(  3) / +.6889758346 9168239842 6263914301 1 d-4      /
      data ai02cs(  4) / +.2891370520 8347564829 6692402323 2 d-5      /
      data ai02cs(  5) / +.2048918589 4690637418 2760534093 1 d-6      /
      data ai02cs(  6) / +.2266668990 4981780645 9327743136 1 d-7      /
      data ai02cs(  7) / +.3396232025 7083863451 5084396952 3 d-8      /
      data ai02cs(  8) / +.4940602388 2249695891 0482449783 5 d-9      /
      data ai02cs(  9) / +.1188914710 7846438342 4084525196 3 d-10     /
      data ai02cs( 10) / -.3149916527 9632413645 3864862961 9 d-10     /
      data ai02cs( 11) / -.1321581184 0447713118 7540739926 7 d-10     /
      data ai02cs( 12) / -.1794178531 5068061177 7943574026 9 d-11     /
      data ai02cs( 13) / +.7180124451 3836662336 7106429346 9 d-12     /
      data ai02cs( 14) / +.3852778382 7421427011 4089801777 6 d-12     /
      data ai02cs( 15) / +.1540086217 5214098269 1325823339 7 d-13     /
      data ai02cs( 16) / -.4150569347 2872220866 2689972015 6 d-13     /
      data ai02cs( 17) / -.9554846698 8283076487 0214494312 5 d-14     /
      data ai02cs( 18) / +.3811680669 3526224207 4605535511 8 d-14     /
      data ai02cs( 19) / +.1772560133 0565263836 0493266675 8 d-14     /
      data ai02cs( 20) / -.3425485619 6772191346 1924790328 2 d-15     /
      data ai02cs( 21) / -.2827623980 5165834849 4205593759 4 d-15     /
      data ai02cs( 22) / +.3461222867 6974610930 9706250813 4 d-16     /
      data ai02cs( 23) / +.4465621420 2967599990 1042054284 3 d-16     /
      data ai02cs( 24) / -.4830504485 9441820712 5525403795 4 d-17     /
      data ai02cs( 25) / -.7233180487 8747539545 6227240924 5 d-17     /
      data ai02cs( 26) / +.9921475412 1736985988 8046093981 0 d-18     /
      data ai02cs( 27) / +.1193650890 8459820855 0439949924 2 d-17     /
      data ai02cs( 28) / -.2488709837 1508072357 2054491660 2 d-18     /
      data ai02cs( 29) / -.1938426454 1609059289 8469781132 6 d-18     /
      data ai02cs( 30) / +.6444656697 3734438687 8301949394 9 d-19     /
      data ai02cs( 31) / +.2886051596 2892243264 8171383073 4 d-19     /
      data ai02cs( 32) / -.1601954907 1749718070 6167156200 7 d-19     /
      data ai02cs( 33) / -.3270815010 5923147208 9193567485 9 d-20     /
      data ai02cs( 34) / +.3686932283 8264091811 4600723939 3 d-20     /
      data ai02cs( 35) / +.1268297648 0309501530 1359529710 9 d-22     /
      data ai02cs( 36) / -.7549825019 3772739076 9636664410 1 d-21     /
      data ai02cs( 37) / +.1502133571 3778353496 3712789053 4 d-21     /
      data ai02cs( 38) / +.1265195883 5096485349 3208799248 3 d-21     /
      data ai02cs( 39) / -.6100998370 0836807086 2940891600 2 d-22     /
      data ai02cs( 40) / -.1268809629 2601282643 6872095924 2 d-22     /
      data ai02cs( 41) / +.1661016099 8907414578 4038487490 5 d-22     /
      data ai02cs( 42) / -.1585194335 7658855793 7970504881 4 d-23     /
      data ai02cs( 43) / -.3302645405 9682178009 5381766755 6 d-23     /
      data ai02cs( 44) / +.1313580902 8392397817 4039623117 4 d-23     /
      data ai02cs( 45) / +.3689040246 6711567933 1425637280 4 d-24     /
      data ai02cs( 46) / -.4210141910 4616891492 1978247249 9 d-24     /
      data ai02cs( 47) / +.4791954591 0828657806 3171401373 0 d-25     /
      data ai02cs( 48) / +.8459470390 2218217952 9971707412 4 d-25     /
      data ai02cs( 49) / -.4039800940 8728324931 4607937181 0 d-25     /
      data ai02cs( 50) / -.6434714653 6504313473 0100850469 5 d-26     /
      data ai02cs( 51) / +.1225743398 8756659903 4464736990 5 d-25     /
      data ai02cs( 52) / -.2934391316 0257089231 9879821175 4 d-26     /
      data ai02cs( 53) / -.1961311309 1949829262 0371205728 9 d-26     /
      data ai02cs( 54) / +.1503520374 8221934241 6229900309 8 d-26     /
      data ai02cs( 55) / -.9588720515 7448265520 3386388206 9 d-28     /
      data ai02cs( 56) / -.3483339380 8170454863 9441108511 4 d-27     /
      data ai02cs( 57) / +.1690903610 2630436730 6244960725 6 d-27     /
      data ai02cs( 58) / +.1982866538 7356030438 9400115718 8 d-28     /
      data ai02cs( 59) / -.5317498081 4918162145 7583002528 4 d-28     /
      data ai02cs( 60) / +.1803306629 8883929462 3501450390 1 d-28     /
      data ai02cs( 61) / +.6213093341 4548931758 8405311242 2 d-29     /
      data ai02cs( 62) / -.7692189292 7721618632 0072806673 0 d-29     /
      data ai02cs( 63) / +.1858252826 1117025426 2556016596 3 d-29     /
      data ai02cs( 64) / +.1237585142 2813957248 9927154554 1 d-29     /
      data ai02cs( 65) / -.1102259120 4092238032 1779478779 2 d-29     /
      data ai02cs( 66) / +.1886287118 0397044900 7787447943 1 d-30     /
      data ai02cs( 67) / +.2160196872 2436589131 4903141406 0 d-30     /
      data ai02cs( 68) / -.1605454124 9197432005 8446594965 5 d-30     /
      data ai02cs( 69) / +.1965352984 5942906039 3884807331 8 d-31     /
c
      data nti0, ntai0, ntai02, xsml / 3*0, 0.d0 /
c
      if (nti0.ne.0) go to 10
      eta = 0.1*sngl(d1mach(3))
      nti0 = initds (bi0cs, 18, eta)
      ntai0 = initds (ai0cs, 46, eta)
      ntai02 = initds (ai02cs, 69, eta)
      xsml = dsqrt (8.0d0*d1mach(3))
c
 10   y = dabs(x)
      if (y.gt.3.0d0) go to 20
c
      dbsi0e = 1.0d0
      if (y.gt.xsml) dbsi0e = dexp(-y) * (2.75d0 +
     1  dcsevl (y*y/4.5d0-1.d0, bi0cs, nti0) )
      return
c
 20   if (y.le.8.d0) dbsi0e = (0.375d0 + dcsevl ((48.d0/y-11.d0)/5.d0,
     1  ai0cs, ntai0))/dsqrt(y)
      if (y.gt.8.d0) dbsi0e = (0.375d0 + dcsevl (16.d0/y-1.d0, ai02cs,
     1  ntai02))/dsqrt(y)
c
      return
      end
      double precision function dbsi1e (x)
c oct 1983 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, bi1cs(17), ai1cs(46), ai12cs(69), xmin,
     1  xsml, y, d1mach, dcsevl, dexp, dsqrt
c     external d1mach, dcsevl, dexp, dsqrt, initds
c
c series for bi1        on the interval  0.          to  9.00000e+00
c                                        with weighted error   1.44e-32
c                                         log weighted error  31.84
c                               significant figures required  31.45
c                                    decimal places required  32.46
c
      data bi1 cs(  1) / -.1971713261 0998597316 1385032181 49 d-2     /
      data bi1 cs(  2) / +.4073488766 7546480608 1553936520 14 d+0     /
      data bi1 cs(  3) / +.3483899429 9959455866 2450377837 87 d-1     /
      data bi1 cs(  4) / +.1545394556 3001236038 5984010584 89 d-2     /
      data bi1 cs(  5) / +.4188852109 8377784129 4588320041 20 d-4     /
      data bi1 cs(  6) / +.7649026764 8362114741 9597039660 69 d-6     /
      data bi1 cs(  7) / +.1004249392 4741178689 1798080372 38 d-7     /
      data bi1 cs(  8) / +.9932207791 9238106481 3712980548 63 d-10    /
      data bi1 cs(  9) / +.7663801791 8447637275 2001716813 49 d-12    /
      data bi1 cs( 10) / +.4741418923 8167394980 3880919481 60 d-14    /
      data bi1 cs( 11) / +.2404114404 0745181799 8631720320 00 d-16    /
      data bi1 cs( 12) / +.1017150500 7093713649 1211007999 99 d-18    /
      data bi1 cs( 13) / +.3645093565 7866949458 4917333333 33 d-21    /
      data bi1 cs( 14) / +.1120574950 2562039344 8106666666 66 d-23    /
      data bi1 cs( 15) / +.2987544193 4468088832 0000000000 00 d-26    /
      data bi1 cs( 16) / +.6973231093 9194709333 3333333333 33 d-29    /
      data bi1 cs( 17) / +.1436794822 0620800000 0000000000 00 d-31    /
c
c series for ai1        on the interval  1.25000e-01 to  3.33333e-01
c                                        with weighted error   2.81e-32
c                                         log weighted error  31.55
c                               significant figures required  29.93
c                                    decimal places required  32.38
c
      data ai1 cs(  1) / -.2846744181 8814786741 0037246830 7 d-1      /
      data ai1 cs(  2) / -.1922953231 4432206510 4444877497 9 d-1      /
      data ai1 cs(  3) / -.6115185857 9437889822 5624991778 5 d-3      /
      data ai1 cs(  4) / -.2069971253 3502277088 8282377797 9 d-4      /
      data ai1 cs(  5) / +.8585619145 8107255655 3694467313 8 d-5      /
      data ai1 cs(  6) / +.1049498246 7115908625 1745399786 0 d-5      /
      data ai1 cs(  7) / -.2918338918 4479022020 9343232669 7 d-6      /
      data ai1 cs(  8) / -.1559378146 6317390001 6068096907 7 d-7      /
      data ai1 cs(  9) / +.1318012367 1449447055 2530287390 9 d-7      /
      data ai1 cs( 10) / -.1448423418 1830783176 3913446781 5 d-8      /
      data ai1 cs( 11) / -.2908512243 9931420948 2504099301 0 d-9      /
      data ai1 cs( 12) / +.1266388917 8753823873 1115969040 3 d-9      /
      data ai1 cs( 13) / -.1664947772 9192206706 2417839858 0 d-10     /
      data ai1 cs( 14) / -.1666653644 6094329760 9593715499 9 d-11     /
      data ai1 cs( 15) / +.1242602414 2907682652 3216847201 7 d-11     /
      data ai1 cs( 16) / -.2731549379 6724323972 5146142863 3 d-12     /
      data ai1 cs( 17) / +.2023947881 6458037807 0026268898 1 d-13     /
      data ai1 cs( 18) / +.7307950018 1168836361 9869812612 3 d-14     /
      data ai1 cs( 19) / -.3332905634 4046749438 1377861713 3 d-14     /
      data ai1 cs( 20) / +.7175346558 5129537435 4225466567 0 d-15     /
      data ai1 cs( 21) / -.6982530324 7962563558 5062922365 6 d-16     /
      data ai1 cs( 22) / -.1299944201 5627607600 6044608058 7 d-16     /
      data ai1 cs( 23) / +.8120942864 2427988920 5467834286 0 d-17     /
      data ai1 cs( 24) / -.2194016207 4107368981 5626664378 3 d-17     /
      data ai1 cs( 25) / +.3630516170 0296548482 7986093233 4 d-18     /
      data ai1 cs( 26) / -.1695139772 4391041663 0686679039 9 d-19     /
      data ai1 cs( 27) / -.1288184829 8979078071 1688253822 2 d-19     /
      data ai1 cs( 28) / +.5694428604 9670527801 0999107310 9 d-20     /
      data ai1 cs( 29) / -.1459597009 0904800565 4550990028 7 d-20     /
      data ai1 cs( 30) / +.2514546010 6757173140 8469133448 5 d-21     /
      data ai1 cs( 31) / -.1844758883 1391248181 6040002901 3 d-22     /
      data ai1 cs( 32) / -.6339760596 2279486419 2860979199 9 d-23     /
      data ai1 cs( 33) / +.3461441102 0310111111 0814662656 0 d-23     /
      data ai1 cs( 34) / -.1017062335 3713935475 9654102357 3 d-23     /
      data ai1 cs( 35) / +.2149877147 0904314459 6250077866 6 d-24     /
      data ai1 cs( 36) / -.3045252425 2386764017 4620617386 6 d-25     /
      data ai1 cs( 37) / +.5238082144 7212859821 7763498666 6 d-27     /
      data ai1 cs( 38) / +.1443583107 0893824464 1678950399 9 d-26     /
      data ai1 cs( 39) / -.6121302074 8900427332 0067071999 9 d-27     /
      data ai1 cs( 40) / +.1700011117 4678184183 4918980266 6 d-27     /
      data ai1 cs( 41) / -.3596589107 9842441585 3521578666 6 d-28     /
      data ai1 cs( 42) / +.5448178578 9484185766 5051306666 6 d-29     /
      data ai1 cs( 43) / -.2731831789 6890849891 6256426666 6 d-30     /
      data ai1 cs( 44) / -.1858905021 7086007157 7190399999 9 d-30     /
      data ai1 cs( 45) / +.9212682974 5139334411 2776533333 3 d-31     /
      data ai1 cs( 46) / -.2813835155 6535611063 7083306666 6 d-31     /
c
c series for ai12       on the interval  0.          to  1.25000e-01
c                                        with weighted error   1.83e-32
c                                         log weighted error  31.74
c                               significant figures required  29.97
c                                    decimal places required  32.66
c
      data ai12cs(  1) / +.2857623501 8280120474 4984594846 9 d-1      /
      data ai12cs(  2) / -.9761097491 3614684077 6516445730 2 d-2      /
      data ai12cs(  3) / -.1105889387 6262371629 1256921277 5 d-3      /
      data ai12cs(  4) / -.3882564808 8776903934 5654477627 4 d-5      /
      data ai12cs(  5) / -.2512236237 8702089252 9452002212 1 d-6      /
      data ai12cs(  6) / -.2631468846 8895195068 3705236523 2 d-7      /
      data ai12cs(  7) / -.3835380385 9642370220 4500678796 8 d-8      /
      data ai12cs(  8) / -.5589743462 1965838068 6811252222 9 d-9      /
      data ai12cs(  9) / -.1897495812 3505412344 9892503323 8 d-10     /
      data ai12cs( 10) / +.3252603583 0154882385 5508067994 9 d-10     /
      data ai12cs( 11) / +.1412580743 6613781331 6336633284 6 d-10     /
      data ai12cs( 12) / +.2035628544 1470895072 2452613684 0 d-11     /
      data ai12cs( 13) / -.7198551776 2459085120 9258989044 6 d-12     /
      data ai12cs( 14) / -.4083551111 0921973182 2849963969 1 d-12     /
      data ai12cs( 15) / -.2101541842 7726643130 1984572746 2 d-13     /
      data ai12cs( 16) / +.4272440016 7119513542 9778833699 7 d-13     /
      data ai12cs( 17) / +.1042027698 4128802764 1741449994 8 d-13     /
      data ai12cs( 18) / -.3814403072 4370078047 6707253539 6 d-14     /
      data ai12cs( 19) / -.1880354775 5107824485 1273453396 3 d-14     /
      data ai12cs( 20) / +.3308202310 9209282827 3190335240 5 d-15     /
      data ai12cs( 21) / +.2962628997 6459501390 6854654205 2 d-15     /
      data ai12cs( 22) / -.3209525921 9934239587 7837353288 7 d-16     /
      data ai12cs( 23) / -.4650305368 4893583255 7128281897 9 d-16     /
      data ai12cs( 24) / +.4414348323 0717079499 4611375964 1 d-17     /
      data ai12cs( 25) / +.7517296310 8421048054 2545808029 5 d-17     /
      data ai12cs( 26) / -.9314178867 3268833756 8484784515 7 d-18     /
      data ai12cs( 27) / -.1242193275 1948909561 1678448869 7 d-17     /
      data ai12cs( 28) / +.2414276719 4548484690 0515390217 6 d-18     /
      data ai12cs( 29) / +.2026944384 0532851789 7192286069 2 d-18     /
      data ai12cs( 30) / -.6394267188 2690977870 4391988681 1 d-19     /
      data ai12cs( 31) / -.3049812452 3730958960 8488450357 1 d-19     /
      data ai12cs( 32) / +.1612841851 6514802251 3462230769 1 d-19     /
      data ai12cs( 33) / +.3560913964 3099250545 1027090462 0 d-20     /
      data ai12cs( 34) / -.3752017947 9364390796 6682800324 6 d-20     /
      data ai12cs( 35) / -.5787037427 0747993459 5198231074 1 d-22     /
      data ai12cs( 36) / +.7759997511 6481619619 8236963209 2 d-21     /
      data ai12cs( 37) / -.1452790897 2022333940 6445987408 5 d-21     /
      data ai12cs( 38) / -.1318225286 7390367021 2192275337 4 d-21     /
      data ai12cs( 39) / +.6116654862 9030707018 7999133171 7 d-22     /
      data ai12cs( 40) / +.1376279762 4271264277 3024338363 4 d-22     /
      data ai12cs( 41) / -.1690837689 9593478849 1983938230 6 d-22     /
      data ai12cs( 42) / +.1430596088 5954331539 8720108538 5 d-23     /
      data ai12cs( 43) / +.3409557828 0905940204 0536772990 2 d-23     /
      data ai12cs( 44) / -.1309457666 2707602278 4573872642 4 d-23     /
      data ai12cs( 45) / -.3940706411 2402574360 9352141755 7 d-24     /
      data ai12cs( 46) / +.4277137426 9808765808 0616679735 2 d-24     /
      data ai12cs( 47) / -.4424634830 9826068819 0028312302 9 d-25     /
      data ai12cs( 48) / -.8734113196 2307149721 1530978874 7 d-25     /
      data ai12cs( 49) / +.4045401335 6835333921 4340414242 8 d-25     /
      data ai12cs( 50) / +.7067100658 0946894656 5160771780 6 d-26     /
      data ai12cs( 51) / -.1249463344 5651052230 0286451860 5 d-25     /
      data ai12cs( 52) / +.2867392244 4034370329 7948339142 6 d-26     /
      data ai12cs( 53) / +.2044292892 5042926702 8177957421 0 d-26     /
      data ai12cs( 54) / -.1518636633 8204625683 7134680291 1 d-26     /
      data ai12cs( 55) / +.8110181098 1875758861 3227910703 7 d-28     /
      data ai12cs( 56) / +.3580379354 7735860911 2717370327 0 d-27     /
      data ai12cs( 57) / -.1692929018 9279025095 9305717544 8 d-27     /
      data ai12cs( 58) / -.2222902499 7024276390 6775852777 4 d-28     /
      data ai12cs( 59) / +.5424535127 1459696550 4860040112 8 d-28     /
      data ai12cs( 60) / -.1787068401 5780186887 6491299330 4 d-28     /
      data ai12cs( 61) / -.6565479068 7228149388 2392943788 0 d-29     /
      data ai12cs( 62) / +.7807013165 0611452809 2206770683 9 d-29     /
      data ai12cs( 63) / -.1816595260 6689797173 7933315222 1 d-29     /
      data ai12cs( 64) / -.1287704952 6600848203 7687559895 9 d-29     /
      data ai12cs( 65) / +.1114548172 9881645474 1370927369 4 d-29     /
      data ai12cs( 66) / -.1808343145 0393369391 5936887668 7 d-30     /
      data ai12cs( 67) / -.2231677718 2037719522 3244822893 9 d-30     /
      data ai12cs( 68) / +.1619029596 0803415106 1790980361 4 d-30     /
      data ai12cs( 69) / -.1834079908 8049414139 0130843921 0 d-31     /
c
      data nti1, ntai1, ntai12, xmin, xsml / 3*0, 2*0.d0 /
c
      if (nti1.ne.0) go to 10
      eta = 0.1*sngl(d1mach(3))
      nti1 = initds (bi1cs, 17, eta)
      ntai1 = initds (ai1cs, 46, eta)
      ntai12 = initds (ai12cs, 69, eta)
c
      xmin = 2.0d0*d1mach(1)
      xsml = dsqrt (8.0d0*d1mach(3))
c
 10   y = dabs(x)
      if (y.gt.3.0d0) go to 20
c
      dbsi1e = 0.0d0
      if (y.eq.0.0d0) return
c
      if (y.le.xmin) call seteru (
     1  38hdbsi1e  dabs(x) so small i1 underflows, 37, 1, 0)
      if (y.gt.xmin) dbsi1e = 0.5d0*x
      if (y.gt.xsml) dbsi1e = x*(0.875d0 + dcsevl (y*y/4.5d0-1.d0,
     1  bi1cs, nti1) )
      dbsi1e = dexp(-y) * dbsi1e
      return
c
 20   if (y.le.8.d0) dbsi1e = (0.375d0 + dcsevl ((48.d0/y-11.d0)/5.d0,
     1  ai1cs, ntai1))/dsqrt(y)
      if (y.gt.8.d0) dbsi1e = (0.375d0 + dcsevl (16.d0/y-1.d0, ai12cs,
     1  ntai12))/dsqrt(y)
      dbsi1e = dsign (dbsi1e, x)
c
      return
      end
      subroutine d9b0mp (x, ampl, theta)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
c
c evaluate the modulus and phase for the bessel j0 and y0 functions.
c
      double precision x, ampl, theta, bm0cs(37), bt02cs(39),
     1  bm02cs(40), bth0cs(44), xmax, pi4, z, d1mach, dcsevl,
     2  dsqrt
c     external d1mach, dcsevl, dsqrt, initds
c
c series for bm0        on the interval  1.56250e-02 to  6.25000e-02
c                                        with weighted error   4.40e-32
c                                         log weighted error  31.36
c                               significant figures required  30.02
c                                    decimal places required  32.14
c
      data bm0 cs(  1) / +.9211656246 8277427125 7376773018 2 d-1      /
      data bm0 cs(  2) / -.1050590997 2719051024 8071637175 5 d-2      /
      data bm0 cs(  3) / +.1470159840 7687597540 5639285095 2 d-4      /
      data bm0 cs(  4) / -.5058557606 0385542233 4792932770 2 d-6      /
      data bm0 cs(  5) / +.2787254538 6324441766 3035613788 1 d-7      /
      data bm0 cs(  6) / -.2062363611 7809148026 1884101897 3 d-8      /
      data bm0 cs(  7) / +.1870214313 1388796751 3817259626 1 d-9      /
      data bm0 cs(  8) / -.1969330971 1356362002 4173077782 5 d-10     /
      data bm0 cs(  9) / +.2325973793 9992754440 1250881805 2 d-11     /
      data bm0 cs( 10) / -.3009520344 9382502728 5122473448 2 d-12     /
      data bm0 cs( 11) / +.4194521333 8506691814 7120676864 6 d-13     /
      data bm0 cs( 12) / -.6219449312 1884458259 7326742956 4 d-14     /
      data bm0 cs( 13) / +.9718260411 3360684696 0176588526 9 d-15     /
      data bm0 cs( 14) / -.1588478585 7010752073 6663596693 7 d-15     /
      data bm0 cs( 15) / +.2700072193 6713088900 8621732445 8 d-16     /
      data bm0 cs( 16) / -.4750092365 2340089924 7750478677 3 d-17     /
      data bm0 cs( 17) / +.8615128162 6043708731 9170374656 0 d-18     /
      data bm0 cs( 18) / -.1605608686 9561448157 4560270335 9 d-18     /
      data bm0 cs( 19) / +.3066513987 3144829751 8853980159 9 d-19     /
      data bm0 cs( 20) / -.5987764223 1939564306 9650561706 6 d-20     /
      data bm0 cs( 21) / +.1192971253 7482483064 8906984106 6 d-20     /
      data bm0 cs( 22) / -.2420969142 0448054894 8468258133 3 d-21     /
      data bm0 cs( 23) / +.4996751760 5106164533 7100287999 9 d-22     /
      data bm0 cs( 24) / -.1047493639 3511585100 9504051199 9 d-22     /
      data bm0 cs( 25) / +.2227786843 7974681010 4818346666 6 d-23     /
      data bm0 cs( 26) / -.4801813239 3981628623 7054293333 3 d-24     /
      data bm0 cs( 27) / +.1047962723 4709599564 7699626666 6 d-24     /
      data bm0 cs( 28) / -.2313858165 6786153251 0126080000 0 d-25     /
      data bm0 cs( 29) / +.5164823088 4626742116 3519999999 9 d-26     /
      data bm0 cs( 30) / -.1164691191 8500653895 2540159999 9 d-26     /
      data bm0 cs( 31) / +.2651788486 0433192829 5833600000 0 d-27     /
      data bm0 cs( 32) / -.6092559503 8257284976 9130666666 6 d-28     /
      data bm0 cs( 33) / +.1411804686 1442593080 3882666666 6 d-28     /
      data bm0 cs( 34) / -.3298094961 2317372457 5061333333 3 d-29     /
      data bm0 cs( 35) / +.7763931143 0740650317 1413333333 3 d-30     /
      data bm0 cs( 36) / -.1841031343 6614584784 2133333333 3 d-30     /
      data bm0 cs( 37) / +.4395880138 5943107371 0079999999 9 d-31     /
c
c series for bth0       on the interval  0.          to  1.56250e-02
c                                        with weighted error   2.66e-32
c                                         log weighted error  31.57
c                               significant figures required  30.67
c                                    decimal places required  32.40
c
      data bth0cs(  1) / -.2490178086 2128936717 7097937899 67 d+0     /
      data bth0cs(  2) / +.4855029960 9623749241 0486155354 85 d-3     /
      data bth0cs(  3) / -.5451183734 5017204950 6562735635 05 d-5     /
      data bth0cs(  4) / +.1355867305 9405964054 3774459299 03 d-6     /
      data bth0cs(  5) / -.5569139890 2227626227 5832184149 20 d-8     /
      data bth0cs(  6) / +.3260903182 4994335304 0042057194 68 d-9     /
      data bth0cs(  7) / -.2491880786 2461341125 2379038779 93 d-10    /
      data bth0cs(  8) / +.2344937742 0882520554 3524135648 91 d-11    /
      data bth0cs(  9) / -.2609653444 4310387762 1775747661 36 d-12    /
      data bth0cs( 10) / +.3335314042 0097395105 8699550149 23 d-13    /
      data bth0cs( 11) / -.4789000044 0572684646 7507705574 09 d-14    /
      data bth0cs( 12) / +.7595617843 6192215972 6425685452 48 d-15    /
      data bth0cs( 13) / -.1313155601 6891440382 7733974876 33 d-15    /
      data bth0cs( 14) / +.2448361834 5240857495 4268207383 55 d-16    /
      data bth0cs( 15) / -.4880572981 0618777683 2567619183 31 d-17    /
      data bth0cs( 16) / +.1032728502 9786316149 2237563612 04 d-17    /
      data bth0cs( 17) / -.2305763381 5057217157 0047445270 25 d-18    /
      data bth0cs( 18) / +.5404444300 1892693993 0171084837 65 d-19    /
      data bth0cs( 19) / -.1324069519 4366572724 1550328823 85 d-19    /
      data bth0cs( 20) / +.3378079562 1371970203 4247921247 22 d-20    /
      data bth0cs( 21) / -.8945762915 7111779003 0269262922 99 d-21    /
      data bth0cs( 22) / +.2451990688 9219317090 8999086514 05 d-21    /
      data bth0cs( 23) / -.6938842287 6866318680 1399331576 57 d-22    /
      data bth0cs( 24) / +.2022827871 4890138392 9463033377 91 d-22    /
      data bth0cs( 25) / -.6062850000 2335483105 7941953717 64 d-23    /
      data bth0cs( 26) / +.1864974896 4037635381 8237883962 70 d-23    /
      data bth0cs( 27) / -.5878373238 4849894560 2450365308 67 d-24    /
      data bth0cs( 28) / +.1895859144 7999563485 5311795035 13 d-24    /
      data bth0cs( 29) / -.6248197937 2258858959 2916207285 65 d-25    /
      data bth0cs( 30) / +.2101790168 4551024686 6386335290 74 d-25    /
      data bth0cs( 31) / -.7208430093 5209253690 8139339924 46 d-26    /
      data bth0cs( 32) / +.2518136389 2474240867 1564059767 46 d-26    /
      data bth0cs( 33) / -.8951804225 8785778806 1439459536 43 d-27    /
      data bth0cs( 34) / +.3235723747 9762298533 2562358685 87 d-27    /
      data bth0cs( 35) / -.1188301051 9855353657 0471441137 96 d-27    /
      data bth0cs( 36) / +.4430628690 7358104820 5792319417 31 d-28    /
      data bth0cs( 37) / -.1676100964 8834829495 7920101356 81 d-28    /
      data bth0cs( 38) / +.6429294692 1207466972 5323939660 88 d-29    /
      data bth0cs( 39) / -.2499226116 6978652421 2072136827 63 d-29    /
      data bth0cs( 40) / +.9839979429 9521955672 8282603553 18 d-30    /
      data bth0cs( 41) / -.3922037524 2408016397 9891316261 58 d-30    /
      data bth0cs( 42) / +.1581810703 0056522138 5906188456 92 d-30    /
      data bth0cs( 43) / -.6452550614 4890715944 3440983654 26 d-31    /
      data bth0cs( 44) / +.2661111136 9199356137 1770183463 67 d-31    /
c
c series for bm02       on the interval  0.          to  1.56250e-02
c                                        with weighted error   4.72e-32
c                                         log weighted error  31.33
c                               significant figures required  30.00
c                                    decimal places required  32.13
c
      data bm02cs(  1) / +.9500415145 2283813693 3086133556 0 d-1      /
      data bm02cs(  2) / -.3801864682 3656709917 4808156685 1 d-3      /
      data bm02cs(  3) / +.2258339301 0314811929 5182992722 4 d-5      /
      data bm02cs(  4) / -.3895725802 3722287647 3062141260 5 d-7      /
      data bm02cs(  5) / +.1246886416 5120816979 3099052972 5 d-8      /
      data bm02cs(  6) / -.6065949022 1025037798 0383505838 7 d-10     /
      data bm02cs(  7) / +.4008461651 4217469910 1527597104 5 d-11     /
      data bm02cs(  8) / -.3350998183 3980942184 6729879457 4 d-12     /
      data bm02cs(  9) / +.3377119716 5174173670 6326434199 6 d-13     /
      data bm02cs( 10) / -.3964585901 6350127005 6935629582 3 d-14     /
      data bm02cs( 11) / +.5286111503 8838572173 8793974473 5 d-15     /
      data bm02cs( 12) / -.7852519083 4508523136 5464024349 3 d-16     /
      data bm02cs( 13) / +.1280300573 3866822010 1163407344 9 d-16     /
      data bm02cs( 14) / -.2263996296 3914297762 8709924488 4 d-17     /
      data bm02cs( 15) / +.4300496929 6567903886 4641029047 7 d-18     /
      data bm02cs( 16) / -.8705749805 1325870797 4753545145 5 d-19     /
      data bm02cs( 17) / +.1865862713 9620951411 8144277205 0 d-19     /
      data bm02cs( 18) / -.4210482486 0930654573 4508697230 1 d-20     /
      data bm02cs( 19) / +.9956676964 2284009915 8162741784 2 d-21     /
      data bm02cs( 20) / -.2457357442 8053133596 0592147854 7 d-21     /
      data bm02cs( 21) / +.6307692160 7620315680 8735370705 9 d-22     /
      data bm02cs( 22) / -.1678773691 4407401426 9333117238 8 d-22     /
      data bm02cs( 23) / +.4620259064 6739044337 7087813608 7 d-23     /
      data bm02cs( 24) / -.1311782266 8603087322 3769340249 6 d-23     /
      data bm02cs( 25) / +.3834087564 1163028277 4792244027 6 d-24     /
      data bm02cs( 26) / -.1151459324 0777412710 7261329357 6 d-24     /
      data bm02cs( 27) / +.3547210007 5233385230 7697134521 3 d-25     /
      data bm02cs( 28) / -.1119218385 8150046462 6435594217 6 d-25     /
      data bm02cs( 29) / +.3611879427 6298378316 9840499425 7 d-26     /
      data bm02cs( 30) / -.1190687765 9133331500 9264176246 3 d-26     /
      data bm02cs( 31) / +.4005094059 4039681318 0247644953 6 d-27     /
      data bm02cs( 32) / -.1373169422 4522123905 9519391601 7 d-27     /
      data bm02cs( 33) / +.4794199088 7425315859 9649152643 7 d-28     /
      data bm02cs( 34) / -.1702965627 6241095840 0699447645 2 d-28     /
      data bm02cs( 35) / +.6149512428 9363300715 0357516132 4 d-29     /
      data bm02cs( 36) / -.2255766896 5818283499 4430023724 2 d-29     /
      data bm02cs( 37) / +.8399707509 2942994860 6165835320 0 d-30     /
      data bm02cs( 38) / -.3172997595 5626023555 6742393615 2 d-30     /
      data bm02cs( 39) / +.1215205298 8812985545 8333302651 4 d-30     /
      data bm02cs( 40) / -.4715852749 7544386930 1321056804 5 d-31     /
c
c series for bt02       on the interval  1.56250e-02 to  6.25000e-02
c                                        with weighted error   2.99e-32
c                                         log weighted error  31.52
c                               significant figures required  30.61
c                                    decimal places required  32.32
c
      data bt02cs(  1) / -.2454829521 3424597462 0504672493 24 d+0     /
      data bt02cs(  2) / +.1254412103 9084615780 7853317782 99 d-2     /
      data bt02cs(  3) / -.3125395041 4871522854 9734467095 71 d-4     /
      data bt02cs(  4) / +.1470977824 9940831164 4534269693 14 d-5     /
      data bt02cs(  5) / -.9954348893 7950033643 4688503511 58 d-7     /
      data bt02cs(  6) / +.8549316673 3203041247 5787113977 51 d-8     /
      data bt02cs(  7) / -.8698975952 6554334557 9855121791 92 d-9     /
      data bt02cs(  8) / +.1005209953 3559791084 5401010821 53 d-9     /
      data bt02cs(  9) / -.1282823060 1708892903 4836236855 44 d-10    /
      data bt02cs( 10) / +.1773170078 1805131705 6557504510 23 d-11    /
      data bt02cs( 11) / -.2617457456 9485577488 6362841809 25 d-12    /
      data bt02cs( 12) / +.4082835138 9972059621 9664812211 03 d-13    /
      data bt02cs( 13) / -.6675166823 9742720054 6067495542 61 d-14    /
      data bt02cs( 14) / +.1136576139 3071629448 3924695499 51 d-14    /
      data bt02cs( 15) / -.2005118962 0647160250 5592664121 17 d-15    /
      data bt02cs( 16) / +.3649797879 4766269635 7205914641 06 d-16    /
      data bt02cs( 17) / -.6830963756 4582303169 3558437888 00 d-17    /
      data bt02cs( 18) / +.1310758314 5670756620 0571042679 46 d-17    /
      data bt02cs( 19) / -.2572336310 1850607778 7571306495 99 d-18    /
      data bt02cs( 20) / +.5152165744 1863959925 2677809493 33 d-19    /
      data bt02cs( 21) / -.1051301756 3758802637 9407414613 33 d-19    /
      data bt02cs( 22) / +.2182038199 1194813847 3010845013 33 d-20    /
      data bt02cs( 23) / -.4600470121 0362160577 2259054933 33 d-21    /
      data bt02cs( 24) / +.9840700692 5466818520 9536511999 99 d-22    /
      data bt02cs( 25) / -.2133403803 5728375844 7359863466 66 d-22    /
      data bt02cs( 26) / +.4683103642 3973365296 0662869333 33 d-23    /
      data bt02cs( 27) / -.1040021369 1985747236 5133823999 99 d-23    /
      data bt02cs( 28) / +.2334910567 7301510051 7777408000 00 d-24    /
      data bt02cs( 29) / -.5295682532 3318615788 0497493333 33 d-25    /
      data bt02cs( 30) / +.1212634195 2959756829 1962879999 99 d-25    /
      data bt02cs( 31) / -.2801889708 2289428760 2756266666 66 d-26    /
      data bt02cs( 32) / +.6529267898 7012873342 5937066666 66 d-27    /
      data bt02cs( 33) / -.1533798006 1873346427 8357333333 33 d-27    /
      data bt02cs( 34) / +.3630588430 6364536682 3594666666 66 d-28    /
      data bt02cs( 35) / -.8656075571 3629122479 1722666666 66 d-29    /
      data bt02cs( 36) / +.2077990997 2536284571 2383999999 99 d-29    /
      data bt02cs( 37) / -.5021117022 1417221674 3253333333 33 d-30    /
      data bt02cs( 38) / +.1220836027 9441714184 1919999999 99 d-30    /
      data bt02cs( 39) / -.2986005626 7039913454 2506666666 66 d-31    /
c
      data pi4 / 0.7853981633 9744830961 5660845819 876 d0 /
      data nbm0, nbt02, nbm02, nbth0, xmax / 4*0, 0.d0 /
c
      if (nbm0.ne.0) go to 10
      eta = 0.1*sngl(d1mach(3))
      nbm0 = initds (bm0cs, 37, eta)
      nbt02 = initds (bt02cs, 39, eta)
      nbm02 = initds (bm02cs, 40, eta)
      nbth0 = initds (bth0cs, 44, eta)
c
      xmax = 1.0d0/d1mach(4)
c
 10   if (x.lt.4.d0) call seteru (22hd9b0mp  x must be ge 4, 22, 1, 2)
c
      if (x.gt.8.d0) go to 20
      z = (128.d0/(x*x) - 5.d0)/3.d0
      ampl = (.75d0 + dcsevl (z, bm0cs, nbm0))/dsqrt(x)
      theta = x - pi4 + dcsevl (z, bt02cs, nbt02)/x
      return
c
 20   if (x.gt.xmax) call seteru (
     1  37hd9b0mp  no precision because x is big, 37, 2, 2)
c
      z = 128.d0/(x*x) - 1.d0
      ampl = (.75d0 + dcsevl (z, bm02cs, nbm02))/dsqrt(x)
      theta = x - pi4 + dcsevl (z, bth0cs, nbth0)/x
      return
c
      end
      subroutine d9b1mp (x, ampl, theta)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
c
c evaluate the modulus and phase for the bessel j1 and y1 functions.
c
      double precision x, ampl, theta, bm1cs(37), bt12cs(39),
     1  bm12cs(40), bth1cs(44), xmax, pi4, z, d1mach, dcsevl,
     2  dsqrt
c     external d1mach, dcsevl, dsqrt, initds
c
c series for bm1        on the interval  1.56250e-02 to  6.25000e-02
c                                        with weighted error   4.91e-32
c                                         log weighted error  31.31
c                               significant figures required  30.04
c                                    decimal places required  32.09
c
      data bm1 cs(  1) / +.1069845452 6180630149 6998530853 8 d+0      /
      data bm1 cs(  2) / +.3274915039 7159649007 2905514344 5 d-2      /
      data bm1 cs(  3) / -.2987783266 8316985920 3044577793 8 d-4      /
      data bm1 cs(  4) / +.8331237177 9919745313 9322266902 3 d-6      /
      data bm1 cs(  5) / -.4112665690 3020073048 9638172549 8 d-7      /
      data bm1 cs(  6) / +.2855344228 7892152207 1975766316 1 d-8      /
      data bm1 cs(  7) / -.2485408305 4156238780 6002659605 5 d-9      /
      data bm1 cs(  8) / +.2543393338 0725824427 4248439717 4 d-10     /
      data bm1 cs(  9) / -.2941045772 8229675234 8975082790 9 d-11     /
      data bm1 cs( 10) / +.3743392025 4939033092 6505615362 6 d-12     /
      data bm1 cs( 11) / -.5149118293 8211672187 2054824352 7 d-13     /
      data bm1 cs( 12) / +.7552535949 8651439080 3404076419 9 d-14     /
      data bm1 cs( 13) / -.1169409706 8288464441 6629062246 4 d-14     /
      data bm1 cs( 14) / +.1896562449 4347915717 2182460506 0 d-15     /
      data bm1 cs( 15) / -.3201955368 6932864206 6477531639 4 d-16     /
      data bm1 cs( 16) / +.5599548399 3162041144 8416990549 3 d-17     /
      data bm1 cs( 17) / -.1010215894 7304324431 1939044454 4 d-17     /
      data bm1 cs( 18) / +.1873844985 7275629833 0204271957 3 d-18     /
      data bm1 cs( 19) / -.3563537470 3285802192 7430143999 9 d-19     /
      data bm1 cs( 20) / +.6931283819 9712383304 2276351999 9 d-20     /
      data bm1 cs( 21) / -.1376059453 4065001522 5140893013 3 d-20     /
      data bm1 cs( 22) / +.2783430784 1070802205 9977932799 9 d-21     /
      data bm1 cs( 23) / -.5727595364 3205616893 4866943999 9 d-22     /
      data bm1 cs( 24) / +.1197361445 9188926725 3575679999 9 d-22     /
      data bm1 cs( 25) / -.2539928509 8918719766 4144042666 6 d-23     /
      data bm1 cs( 26) / +.5461378289 6572959730 6961919999 9 d-24     /
      data bm1 cs( 27) / -.1189211341 7733202889 8628949333 3 d-24     /
      data bm1 cs( 28) / +.2620150977 3400815949 5782400000 0 d-25     /
      data bm1 cs( 29) / -.5836810774 2556859019 2093866666 6 d-26     /
      data bm1 cs( 30) / +.1313743500 0805957734 2361599999 9 d-26     /
      data bm1 cs( 31) / -.2985814622 5103803553 3277866666 6 d-27     /
      data bm1 cs( 32) / +.6848390471 3346049376 2559999999 9 d-28     /
      data bm1 cs( 33) / -.1584401568 2224767211 9296000000 0 d-28     /
      data bm1 cs( 34) / +.3695641006 5709380543 0101333333 3 d-29     /
      data bm1 cs( 35) / -.8687115921 1446682430 1226666666 6 d-30     /
      data bm1 cs( 36) / +.2057080846 1587634629 2906666666 6 d-30     /
      data bm1 cs( 37) / -.4905225761 1162255185 2373333333 3 d-31     /
c
c series for bt12       on the interval  1.56250e-02 to  6.25000e-02
c                                        with weighted error   3.33e-32
c                                         log weighted error  31.48
c                               significant figures required  31.05
c                                    decimal places required  32.27
c
      data bt12cs(  1) / +.7382386012 8742974662 6208397927 64 d+0     /
      data bt12cs(  2) / -.3336111317 4483906384 4701476811 89 d-2     /
      data bt12cs(  3) / +.6146345488 8046964698 5148994201 86 d-4     /
      data bt12cs(  4) / -.2402458516 1602374264 9776354695 68 d-5     /
      data bt12cs(  5) / +.1466355557 7509746153 2105919972 04 d-6     /
      data bt12cs(  6) / -.1184191730 5589180567 0051475049 83 d-7     /
      data bt12cs(  7) / +.1157419896 3919197052 1254663030 55 d-8     /
      data bt12cs(  8) / -.1300116112 9439187449 3660077945 71 d-9     /
      data bt12cs(  9) / +.1624539114 1361731937 7421662736 67 d-10    /
      data bt12cs( 10) / -.2208963682 1403188752 1554417701 28 d-11    /
      data bt12cs( 11) / +.3218030425 8553177090 4743586537 78 d-12    /
      data bt12cs( 12) / -.4965314793 2768480785 5520211353 81 d-13    /
      data bt12cs( 13) / +.8043890043 2847825985 5588826393 17 d-14    /
      data bt12cs( 14) / -.1358912131 0161291384 6947126822 82 d-14    /
      data bt12cs( 15) / +.2381050439 7147214869 6765296059 73 d-15    /
      data bt12cs( 16) / -.4308146636 3849106724 4712414207 99 d-16    /
      data bt12cs( 17) / +.8020254403 2771002434 9935125504 00 d-17    /
      data bt12cs( 18) / -.1531631064 2462311864 2300274687 99 d-17    /
      data bt12cs( 19) / +.2992860635 2715568924 0730405546 66 d-18    /
      data bt12cs( 20) / -.5970996465 8085443393 8156366506 66 d-19    /
      data bt12cs( 21) / +.1214028966 9415185024 1608526506 66 d-19    /
      data bt12cs( 22) / -.2511511469 6612948901 0069777066 66 d-20    /
      data bt12cs( 23) / +.5279056717 0328744850 7383807999 99 d-21    /
      data bt12cs( 24) / -.1126050922 7550498324 3611613866 66 d-21    /
      data bt12cs( 25) / +.2434827735 9576326659 6634624000 00 d-22    /
      data bt12cs( 26) / -.5331726123 6931800130 0384426666 66 d-23    /
      data bt12cs( 27) / +.1181361505 9707121039 2059903999 99 d-23    /
      data bt12cs( 28) / -.2646536828 3353523514 8567893333 33 d-24    /
      data bt12cs( 29) / +.5990339404 1361503945 5778133333 33 d-25    /
      data bt12cs( 30) / -.1369085463 0829503109 1363839999 99 d-25    /
      data bt12cs( 31) / +.3157679015 4380228326 4136533333 33 d-26    /
      data bt12cs( 32) / -.7345791508 2084356491 4005333333 33 d-27    /
      data bt12cs( 33) / +.1722808148 0722747930 7059200000 00 d-27    /
      data bt12cs( 34) / -.4071690796 1286507941 0688000000 00 d-28    /
      data bt12cs( 35) / +.9693474513 6779622700 3733333333 33 d-29    /
      data bt12cs( 36) / -.2323763633 7765716765 3546666666 66 d-29    /
      data bt12cs( 37) / +.5607451067 3522029406 8906666666 66 d-30    /
      data bt12cs( 38) / -.1361646539 1539005860 5226666666 66 d-30    /
      data bt12cs( 39) / +.3326310923 3894654388 9066666666 66 d-31    /
c
c series for bm12       on the interval  0.          to  1.56250e-02
c                                        with weighted error   5.01e-32
c                                         log weighted error  31.30
c                               significant figures required  29.99
c                                    decimal places required  32.10
c
      data bm12cs(  1) / +.9807979156 2330500272 7209354693 7 d-1      /
      data bm12cs(  2) / +.1150961189 5046853061 7548348460 2 d-2      /
      data bm12cs(  3) / -.4312482164 3382054098 8935809773 2 d-5      /
      data bm12cs(  4) / +.5951839610 0888163078 1302980183 2 d-7      /
      data bm12cs(  5) / -.1704844019 8269098574 0070158647 8 d-8      /
      data bm12cs(  6) / +.7798265413 6111095086 5817382740 1 d-10     /
      data bm12cs(  7) / -.4958986126 7664158094 9175495186 5 d-11     /
      data bm12cs(  8) / +.4038432416 4211415168 3820226514 4 d-12     /
      data bm12cs(  9) / -.3993046163 7251754457 6548384664 5 d-13     /
      data bm12cs( 10) / +.4619886183 1189664943 1334243277 5 d-14     /
      data bm12cs( 11) / -.6089208019 0953833013 4547261933 3 d-15     /
      data bm12cs( 12) / +.8960930916 4338764821 5704804124 9 d-16     /
      data bm12cs( 13) / -.1449629423 9420231229 1651891892 5 d-16     /
      data bm12cs( 14) / +.2546463158 5377760561 6514964806 8 d-17     /
      data bm12cs( 15) / -.4809472874 6478364442 5926371862 0 d-18     /
      data bm12cs( 16) / +.9687684668 2925990490 8727583912 4 d-19     /
      data bm12cs( 17) / -.2067213372 2779660232 4503811755 1 d-19     /
      data bm12cs( 18) / +.4646651559 1503847318 0276780959 0 d-20     /
      data bm12cs( 19) / -.1094966128 8483341382 4135132833 9 d-20     /
      data bm12cs( 20) / +.2693892797 2886828609 0570761278 5 d-21     /
      data bm12cs( 21) / -.6894992910 9303744778 1897002685 7 d-22     /
      data bm12cs( 22) / +.1830268262 7520629098 9066855474 0 d-22     /
      data bm12cs( 23) / -.5025064246 3519164281 5611355322 4 d-23     /
      data bm12cs( 24) / +.1423545194 4548060396 3169363419 4 d-23     /
      data bm12cs( 25) / -.4152191203 6164503880 6888676980 1 d-24     /
      data bm12cs( 26) / +.1244609201 5039793258 8233007654 7 d-24     /
      data bm12cs( 27) / -.3827336370 5693042994 3191866128 6 d-25     /
      data bm12cs( 28) / +.1205591357 8156175353 7472398183 5 d-25     /
      data bm12cs( 29) / -.3884536246 3764880764 3185936112 4 d-26     /
      data bm12cs( 30) / +.1278689528 7204097219 0489528346 1 d-26     /
      data bm12cs( 31) / -.4295146689 4479462720 6193691591 2 d-27     /
      data bm12cs( 32) / +.1470689117 8290708864 5680270798 3 d-27     /
      data bm12cs( 33) / -.5128315665 1060731281 8037401779 6 d-28     /
      data bm12cs( 34) / +.1819509585 4711693854 8143737328 6 d-28     /
      data bm12cs( 35) / -.6563031314 8419808676 1863505037 3 d-29     /
      data bm12cs( 36) / +.2404898976 9199606531 9891487583 4 d-29     /
      data bm12cs( 37) / -.8945966744 6906124732 3495824297 9 d-30     /
      data bm12cs( 38) / +.3376085160 6572310266 3714897824 0 d-30     /
      data bm12cs( 39) / -.1291791454 6206563609 1309991696 6 d-30     /
      data bm12cs( 40) / +.5008634462 9588105206 8495150125 4 d-31     /
c
c series for bth1       on the interval  0.          to  1.56250e-02
c                                        with weighted error   2.82e-32
c                                         log weighted error  31.55
c                               significant figures required  31.12
c                                    decimal places required  32.37
c
      data bth1cs(  1) / +.7474995720 3587276055 4434839696 95 d+0     /
      data bth1cs(  2) / -.1240077714 4651711252 5457775413 84 d-2     /
      data bth1cs(  3) / +.9925244240 4424527376 6414976895 92 d-5     /
      data bth1cs(  4) / -.2030369073 7159711052 4193753756 08 d-6     /
      data bth1cs(  5) / +.7535961770 5690885712 1840175836 29 d-8     /
      data bth1cs(  6) / -.4166161271 5343550107 6300238562 28 d-9     /
      data bth1cs(  7) / +.3070161807 0834890481 2451020912 16 d-10    /
      data bth1cs(  8) / -.2817849963 7605213992 3240088839 24 d-11    /
      data bth1cs(  9) / +.3079069673 9040295476 0281468216 47 d-12    /
      data bth1cs( 10) / -.3880330026 2803434112 7873475547 81 d-13    /
      data bth1cs( 11) / +.5509603960 8630904934 5617262085 62 d-14    /
      data bth1cs( 12) / -.8659006076 8383779940 1033989539 94 d-15    /
      data bth1cs( 13) / +.1485604914 1536749003 4236890606 83 d-15    /
      data bth1cs( 14) / -.2751952981 5904085805 3712121250 09 d-16    /
      data bth1cs( 15) / +.5455079609 0481089625 0362236409 23 d-17    /
      data bth1cs( 16) / -.1148653450 1983642749 5436310271 77 d-17    /
      data bth1cs( 17) / +.2553521337 7973900223 1990525335 22 d-18    /
      data bth1cs( 18) / -.5962149019 7413450395 7682879078 49 d-19    /
      data bth1cs( 19) / +.1455662290 2372718620 2883020058 33 d-19    /
      data bth1cs( 20) / -.3702218542 2450538201 5797760195 93 d-20    /
      data bth1cs( 21) / +.9776307412 5345357664 1684345179 24 d-21    /
      data bth1cs( 22) / -.2672682163 9668488468 7237753930 52 d-21    /
      data bth1cs( 23) / +.7545330038 4983271794 0381906557 64 d-22    /
      data bth1cs( 24) / -.2194789991 9802744897 8923833716 47 d-22    /
      data bth1cs( 25) / +.6564839462 3955262178 9069998174 93 d-23    /
      data bth1cs( 26) / -.2015560429 8370207570 7840768695 19 d-23    /
      data bth1cs( 27) / +.6341776855 6776143492 1446671856 70 d-24    /
      data bth1cs( 28) / -.2041927788 5337895634 8137699555 91 d-24    /
      data bth1cs( 29) / +.6719146422 0720567486 6589800185 51 d-25    /
      data bth1cs( 30) / -.2256907911 0207573595 7090036873 36 d-25    /
      data bth1cs( 31) / +.7729771989 2989706370 9269598719 29 d-26    /
      data bth1cs( 32) / -.2696744451 2294640913 2114240809 20 d-26    /
      data bth1cs( 33) / +.9574934451 8502698072 2955219336 27 d-27    /
      data bth1cs( 34) / -.3456916844 8890113000 1756808276 27 d-27    /
      data bth1cs( 35) / +.1268123481 7398436504 2119862383 74 d-27    /
      data bth1cs( 36) / -.4723253663 0722639860 4649937134 45 d-28    /
      data bth1cs( 37) / +.1785000847 8186376177 8586197964 17 d-28    /
      data bth1cs( 38) / -.6840436100 4510395406 2152235667 46 d-29    /
      data bth1cs( 39) / +.2656602867 1720419358 2934226722 12 d-29    /
      data bth1cs( 40) / -.1045040252 7914452917 7141614846 70 d-29    /
      data bth1cs( 41) / +.4161829082 5377144306 8619171970 64 d-30    /
      data bth1cs( 42) / -.1677163920 3643714856 5013478828 87 d-30    /
      data bth1cs( 43) / +.6836199777 6664389173 5359280285 28 d-31    /
      data bth1cs( 44) / -.2817224786 1233641166 7395746228 10 d-31    /
c
      data pi4 / 0.7853981633 9744830961 5660845819 876 d0 /
      data nbm1, nbt12, nbm12, nbth1, xmax / 4*0, 0.d0 /
c
      if (nbm1.ne.0) go to 10
      eta = 0.1*sngl(d1mach(3))
      nbm1 = initds (bm1cs, 37, eta)
      nbt12 = initds (bt12cs, 39, eta)
      nbm12 = initds (bm12cs, 40, eta)
      nbth1 = initds (bth1cs, 44, eta)
c
      xmax = 1.0d0/d1mach(4)
c
 10   if (x.lt.4.d0) call seteru (22hd9b1mp  x must be ge 4, 22, 1, 2)
c
      if (x.gt.8.d0) go to 20
      z = (128.d0/(x*x) - 5.d0)/3.d0
      ampl = (0.75d0 + dcsevl (z, bm1cs, nbm1))/dsqrt(x)
      theta = x - 3.d0*pi4 + dcsevl (z, bt12cs, nbt12)/x
      return
c
 20   if (x.gt.xmax) call seteru (
     1  37hd9b1mp  no precision because x is big, 37, 2, 2)
c
      z = 128.d0/(x*x) - 1.d0
      ampl = (0.75d0 + dcsevl (z, bm12cs, nbm12))/dsqrt(x)
      theta = x - 3.d0*pi4 + dcsevl (z, bth1cs, nbth1)/x
      return
c
      end
      double precision function dbsk0e (x)
c july 1980 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, bk0cs(16), ak0cs(38), ak02cs(33),
     1  xsml, y, d1mach, dcsevl, dbesi0, dexp, dlog, dsqrt
c     external d1mach, dbesi0, dcsevl, dexp, dlog, dsqrt, initds
c
c series for bk0        on the interval  0.          to  4.00000e+00
c                                        with weighted error   3.08e-33
c                                         log weighted error  32.51
c                               significant figures required  32.05
c                                    decimal places required  33.11
c
      data bk0 cs(  1) / -.3532739323 3902768720 1140060063 153 d-1    /
      data bk0 cs(  2) / +.3442898999 2462848688 6344927529 213 d+0    /
      data bk0 cs(  3) / +.3597993651 5361501626 5721303687 231 d-1    /
      data bk0 cs(  4) / +.1264615411 4469259233 8479508673 447 d-2    /
      data bk0 cs(  5) / +.2286212103 1194517860 8269830297 585 d-4    /
      data bk0 cs(  6) / +.2534791079 0261494573 0790013428 354 d-6    /
      data bk0 cs(  7) / +.1904516377 2202088589 7214059381 366 d-8    /
      data bk0 cs(  8) / +.1034969525 7633624585 1008317853 089 d-10   /
      data bk0 cs(  9) / +.4259816142 7910825765 2445327170 133 d-13   /
      data bk0 cs( 10) / +.1374465435 8807508969 4238325440 000 d-15   /
      data bk0 cs( 11) / +.3570896528 5083735909 9688597333 333 d-18   /
      data bk0 cs( 12) / +.7631643660 1164373766 7498666666 666 d-21   /
      data bk0 cs( 13) / +.1365424988 4407818590 8053333333 333 d-23   /
      data bk0 cs( 14) / +.2075275266 9066680831 9999999999 999 d-26   /
      data bk0 cs( 15) / +.2712814218 0729856000 0000000000 000 d-29   /
      data bk0 cs( 16) / +.3082593887 9146666666 6666666666 666 d-32   /
c
c series for ak0        on the interval  1.25000e-01 to  5.00000e-01
c                                        with weighted error   2.85e-32
c                                         log weighted error  31.54
c                               significant figures required  30.19
c                                    decimal places required  32.33
c
      data ak0 cs(  1) / -.7643947903 3279414240 8297827008 8 d-1      /
      data ak0 cs(  2) / -.2235652605 6998190520 2309555079 1 d-1      /
      data ak0 cs(  3) / +.7734181154 6938582353 0061817404 7 d-3      /
      data ak0 cs(  4) / -.4281006688 8860994644 5214643541 6 d-4      /
      data ak0 cs(  5) / +.3081700173 8629747436 5001482666 0 d-5      /
      data ak0 cs(  6) / -.2639367222 0096649740 6744889272 3 d-6      /
      data ak0 cs(  7) / +.2563713036 4034692062 9408826574 2 d-7      /
      data ak0 cs(  8) / -.2742705549 9002012638 5721191524 4 d-8      /
      data ak0 cs(  9) / +.3169429658 0974995920 8083287340 3 d-9      /
      data ak0 cs( 10) / -.3902353286 9621841416 0106571796 2 d-10     /
      data ak0 cs( 11) / +.5068040698 1885754020 5009212728 6 d-11     /
      data ak0 cs( 12) / -.6889574741 0078706795 4171355798 4 d-12     /
      data ak0 cs( 13) / +.9744978497 8259176913 8820133683 1 d-13     /
      data ak0 cs( 14) / -.1427332841 8845485053 8985534012 2 d-13     /
      data ak0 cs( 15) / +.2156412571 0214630395 5806297652 7 d-14     /
      data ak0 cs( 16) / -.3349654255 1495627721 8878205853 0 d-15     /
      data ak0 cs( 17) / +.5335260216 9529116921 4528039260 1 d-16     /
      data ak0 cs( 18) / -.8693669980 8907538076 3962237883 7 d-17     /
      data ak0 cs( 19) / +.1446404347 8622122278 8776344234 6 d-17     /
      data ak0 cs( 20) / -.2452889825 5001296824 0467875157 3 d-18     /
      data ak0 cs( 21) / +.4233754526 2321715728 2170634240 0 d-19     /
      data ak0 cs( 22) / -.7427946526 4544641956 9534129493 3 d-20     /
      data ak0 cs( 23) / +.1323150529 3926668662 7796746240 0 d-20     /
      data ak0 cs( 24) / -.2390587164 7396494513 3598146559 9 d-21     /
      data ak0 cs( 25) / +.4376827585 9232261401 6571255466 6 d-22     /
      data ak0 cs( 26) / -.8113700607 3451180593 3901141333 3 d-23     /
      data ak0 cs( 27) / +.1521819913 8321729583 1037815466 6 d-23     /
      data ak0 cs( 28) / -.2886041941 4833977702 3595861333 3 d-24     /
      data ak0 cs( 29) / +.5530620667 0547179799 9261013333 3 d-25     /
      data ak0 cs( 30) / -.1070377329 2498987285 9163306666 6 d-25     /
      data ak0 cs( 31) / +.2091086893 1423843002 9632853333 3 d-26     /
      data ak0 cs( 32) / -.4121713723 6462038274 1026133333 3 d-27     /
      data ak0 cs( 33) / +.8193483971 1213076401 3568000000 0 d-28     /
      data ak0 cs( 34) / -.1642000275 4592977267 8075733333 3 d-28     /
      data ak0 cs( 35) / +.3316143281 4802271958 9034666666 6 d-29     /
      data ak0 cs( 36) / -.6746863644 1452959410 8586666666 6 d-30     /
      data ak0 cs( 37) / +.1382429146 3184246776 3541333333 3 d-30     /
      data ak0 cs( 38) / -.2851874167 3598325708 1173333333 3 d-31     /
c
c series for ak02       on the interval  0.          to  1.25000e-01
c                                        with weighted error   2.30e-32
c                                         log weighted error  31.64
c                               significant figures required  29.68
c                                    decimal places required  32.40
c
      data ak02cs(  1) / -.1201869826 3075922398 3934621245 2 d-1      /
      data ak02cs(  2) / -.9174852691 0256953106 5256107571 3 d-2      /
      data ak02cs(  3) / +.1444550931 7750058210 4884387805 7 d-3      /
      data ak02cs(  4) / -.4013614175 4357097286 7102107787 9 d-5      /
      data ak02cs(  5) / +.1567831810 8523106725 9034899033 3 d-6      /
      data ak02cs(  6) / -.7770110438 5217377103 1579975446 0 d-8      /
      data ak02cs(  7) / +.4611182576 1797178825 3313052958 6 d-9      /
      data ak02cs(  8) / -.3158592997 8605657705 2666580330 9 d-10     /
      data ak02cs(  9) / +.2435018039 3650411278 3588781432 9 d-11     /
      data ak02cs( 10) / -.2074331387 3983478977 0985337350 6 d-12     /
      data ak02cs( 11) / +.1925787280 5899170847 4273650469 3 d-13     /
      data ak02cs( 12) / -.1927554805 8389561036 0034718221 8 d-14     /
      data ak02cs( 13) / +.2062198029 1978182782 8523786964 4 d-15     /
      data ak02cs( 14) / -.2341685117 5792424026 0364019507 1 d-16     /
      data ak02cs( 15) / +.2805902810 6430422468 1517882845 8 d-17     /
      data ak02cs( 16) / -.3530507631 1618079458 1548246357 3 d-18     /
      data ak02cs( 17) / +.4645295422 9351082674 2421633706 6 d-19     /
      data ak02cs( 18) / -.6368625941 3442664739 2205346133 3 d-20     /
      data ak02cs( 19) / +.9069521310 9865155676 2234880000 0 d-21     /
      data ak02cs( 20) / -.1337974785 4236907398 4500531199 9 d-21     /
      data ak02cs( 21) / +.2039836021 8599523155 2208896000 0 d-22     /
      data ak02cs( 22) / -.3207027481 3678405000 6086997333 3 d-23     /
      data ak02cs( 23) / +.5189744413 6623099636 2635946666 6 d-24     /
      data ak02cs( 24) / -.8629501497 5405721929 6460799999 9 d-25     /
      data ak02cs( 25) / +.1472161183 1025598552 0803840000 0 d-25     /
      data ak02cs( 26) / -.2573069023 8670112838 1235199999 9 d-26     /
      data ak02cs( 27) / +.4601774086 6435165873 7664000000 0 d-27     /
      data ak02cs( 28) / -.8411555324 2010937371 3066666666 6 d-28     /
      data ak02cs( 29) / +.1569806306 6353689393 0154666666 6 d-28     /
      data ak02cs( 30) / -.2988226453 0057577889 7919999999 9 d-29     /
      data ak02cs( 31) / +.5796831375 2168365206 1866666666 6 d-30     /
      data ak02cs( 32) / -.1145035994 3476813321 5573333333 3 d-30     /
      data ak02cs( 33) / +.2301266594 2496828020 0533333333 3 d-31     /
      data ntk0, ntak0, ntak02, xsml / 3*0, 0.0d0 /
c
      if (ntk0.ne.0) go to 10
      eta = 0.1*sngl(d1mach(3))
      ntk0 = initds (bk0cs, 16, eta)
      ntak0 = initds (ak0cs, 38, eta)
      ntak02 = initds (ak02cs, 33, eta)
      xsml = dsqrt (4.0d0*d1mach(3))
c
 10   if (x.le.0.d0) call seteru (29hdbsk0e  x is zero or negative, 29,
     1  1, 2)
      if (x.gt.2.0d0) go to 20
c
      y = 0.d0
      if (x.gt.xsml) y = x*x
      dbsk0e = dexp(x)*(-dlog(0.5d0*x)*dbesi0(x) - 0.25d0 +
     1  dcsevl (.5d0*y-1.d0, bk0cs, ntk0))
      return
c
 20   if (x.le.8.d0) dbsk0e = (1.25d0 + dcsevl ((16.d0/x-5.d0)/3.d0,
     1  ak0cs, ntak0))/dsqrt(x)
      if (x.gt.8.d0) dbsk0e = (1.25d0 +
     1  dcsevl (16.d0/x-1.d0, ak02cs, ntak02))/dsqrt(x)
c
      return
      end
      double precision function dbsk1e (x)
c july 1980 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, bk1cs(16), ak1cs(38), ak12cs(33), xmin,
     1  xsml, y, d1mach, dcsevl, dbesi1, dexp, dlog, dsqrt
c     external d1mach, dbesi1, dcsevl, dexp, dlog, dsqrt, initds
c
c series for bk1        on the interval  0.          to  4.00000e+00
c                                        with weighted error   9.16e-32
c                                         log weighted error  31.04
c                               significant figures required  30.61
c                                    decimal places required  31.64
c
      data bk1 cs(  1) / +.2530022733 8947770532 5311208685 33 d-1     /
      data bk1 cs(  2) / -.3531559607 7654487566 7238316918 01 d+0     /
      data bk1 cs(  3) / -.1226111808 2265714823 4790679300 42 d+0     /
      data bk1 cs(  4) / -.6975723859 6398643501 8129202960 83 d-2     /
      data bk1 cs(  5) / -.1730288957 5130520630 1765073689 79 d-3     /
      data bk1 cs(  6) / -.2433406141 5659682349 6007350301 64 d-5     /
      data bk1 cs(  7) / -.2213387630 7347258558 3152525451 26 d-7     /
      data bk1 cs(  8) / -.1411488392 6335277610 9583302126 08 d-9     /
      data bk1 cs(  9) / -.6666901694 1993290060 8537512643 73 d-12    /
      data bk1 cs( 10) / -.2427449850 5193659339 2631968648 53 d-14    /
      data bk1 cs( 11) / -.7023863479 3862875971 7837971200 00 d-17    /
      data bk1 cs( 12) / -.1654327515 5100994675 4910293333 33 d-19    /
      data bk1 cs( 13) / -.3233834745 9944491991 8933333333 33 d-22    /
      data bk1 cs( 14) / -.5331275052 9265274999 4666666666 66 d-25    /
      data bk1 cs( 15) / -.7513040716 2157226666 6666666666 66 d-28    /
      data bk1 cs( 16) / -.9155085717 6541866666 6666666666 66 d-31    /
c
c series for ak1        on the interval  1.25000e-01 to  5.00000e-01
c                                        with weighted error   3.07e-32
c                                         log weighted error  31.51
c                               significant figures required  30.71
c                                    decimal places required  32.30
c
      data ak1 cs(  1) / +.2744313406 9738829695 2576662272 66 d+0     /
      data ak1 cs(  2) / +.7571989953 1993678170 8923781492 90 d-1     /
      data ak1 cs(  3) / -.1441051556 4754061229 8531161756 25 d-2     /
      data ak1 cs(  4) / +.6650116955 1257479394 2513854770 36 d-4     /
      data ak1 cs(  5) / -.4369984709 5201407660 5808450891 67 d-5     /
      data ak1 cs(  6) / +.3540277499 7630526799 4171390085 34 d-6     /
      data ak1 cs(  7) / -.3311163779 2932920208 9826882457 04 d-7     /
      data ak1 cs(  8) / +.3445977581 9010534532 3114997709 92 d-8     /
      data ak1 cs(  9) / -.3898932347 4754271048 9819374927 58 d-9     /
      data ak1 cs( 10) / +.4720819750 4658356400 9474493390 05 d-10    /
      data ak1 cs( 11) / -.6047835662 8753562345 3735915628 90 d-11    /
      data ak1 cs( 12) / +.8128494874 8658747888 1938379856 63 d-12    /
      data ak1 cs( 13) / -.1138694574 7147891428 9239159510 42 d-12    /
      data ak1 cs( 14) / +.1654035840 8462282325 9729482050 90 d-13    /
      data ak1 cs( 15) / -.2480902567 7068848221 5160104405 33 d-14    /
      data ak1 cs( 16) / +.3829237890 7024096948 4292272991 57 d-15    /
      data ak1 cs( 17) / -.6064734104 0012418187 7682103773 86 d-16    /
      data ak1 cs( 18) / +.9832425623 2648616038 1940046506 66 d-17    /
      data ak1 cs( 19) / -.1628416873 8284380035 6666201156 26 d-17    /
      data ak1 cs( 20) / +.2750153649 6752623718 2841203370 66 d-18    /
      data ak1 cs( 21) / -.4728966646 3953250924 2810695680 00 d-19    /
      data ak1 cs( 22) / +.8268150002 8109932722 3920503466 66 d-20    /
      data ak1 cs( 23) / -.1468140513 6624956337 1939648853 33 d-20    /
      data ak1 cs( 24) / +.2644763926 9208245978 0858948266 66 d-21    /
      data ak1 cs( 25) / -.4829015756 4856387897 9698688000 00 d-22    /
      data ak1 cs( 26) / +.8929302074 3610130180 6563327999 99 d-23    /
      data ak1 cs( 27) / -.1670839716 8972517176 9977514666 66 d-23    /
      data ak1 cs( 28) / +.3161645603 4040694931 3686186666 66 d-24    /
      data ak1 cs( 29) / -.6046205531 2274989106 5064106666 66 d-25    /
      data ak1 cs( 30) / +.1167879894 2042732700 7184213333 33 d-25    /
      data ak1 cs( 31) / -.2277374158 2653996232 8678400000 00 d-26    /
      data ak1 cs( 32) / +.4481109730 0773675795 3058133333 33 d-27    /
      data ak1 cs( 33) / -.8893288476 9020194062 3360000000 00 d-28    /
      data ak1 cs( 34) / +.1779468001 8850275131 3920000000 00 d-28    /
      data ak1 cs( 35) / -.3588455596 7329095821 9946666666 66 d-29    /
      data ak1 cs( 36) / +.7290629049 2694257991 6799999999 99 d-30    /
      data ak1 cs( 37) / -.1491844984 5546227073 0240000000 00 d-30    /
      data ak1 cs( 38) / +.3073657387 2934276300 7999999999 99 d-31    /
c
c series for ak12       on the interval  0.          to  1.25000e-01
c                                        with weighted error   2.41e-32
c                                         log weighted error  31.62
c                               significant figures required  30.25
c                                    decimal places required  32.38
c
      data ak12cs(  1) / +.6379308343 7390010366 0048853410 2 d-1      /
      data ak12cs(  2) / +.2832887813 0497209358 3503028470 8 d-1      /
      data ak12cs(  3) / -.2475370673 9052503454 1454556673 2 d-3      /
      data ak12cs(  4) / +.5771972451 6072488204 7097662576 3 d-5      /
      data ak12cs(  5) / -.2068939219 5365483027 4553319655 2 d-6      /
      data ak12cs(  6) / +.9739983441 3818041803 0921309788 7 d-8      /
      data ak12cs(  7) / -.5585336140 3806249846 8889551112 9 d-9      /
      data ak12cs(  8) / +.3732996634 0461852402 2121285473 1 d-10     /
      data ak12cs(  9) / -.2825051961 0232254451 3506575492 8 d-11     /
      data ak12cs( 10) / +.2372019002 4841441736 4349695548 6 d-12     /
      data ak12cs( 11) / -.2176677387 9917539792 6830166793 8 d-13     /
      data ak12cs( 12) / +.2157914161 6160324539 3956268970 6 d-14     /
      data ak12cs( 13) / -.2290196930 7182692759 9155133815 4 d-15     /
      data ak12cs( 14) / +.2582885729 8232749619 1993956522 6 d-16     /
      data ak12cs( 15) / -.3076752641 2684631876 2109817344 0 d-17     /
      data ak12cs( 16) / +.3851487721 2804915970 9489684479 9 d-18     /
      data ak12cs( 17) / -.5044794897 6415289771 1728250880 0 d-19     /
      data ak12cs( 18) / +.6888673850 4185442370 1829222399 9 d-20     /
      data ak12cs( 19) / -.9775041541 9501183030 0213248000 0 d-21     /
      data ak12cs( 20) / +.1437416218 5238364610 0165973333 3 d-21     /
      data ak12cs( 21) / -.2185059497 3443473734 9973333333 3 d-22     /
      data ak12cs( 22) / +.3426245621 8092206316 4538880000 0 d-23     /
      data ak12cs( 23) / -.5531064394 2464082325 0124800000 0 d-24     /
      data ak12cs( 24) / +.9176601505 6859954037 8282666666 6 d-25     /
      data ak12cs( 25) / -.1562287203 6180249114 4874666666 6 d-25     /
      data ak12cs( 26) / +.2725419375 4843331323 4943999999 9 d-26     /
      data ak12cs( 27) / -.4865674910 0748279923 7802666666 6 d-27     /
      data ak12cs( 28) / +.8879388552 7235025873 5786666666 6 d-28     /
      data ak12cs( 29) / -.1654585918 0392575489 3653333333 3 d-28     /
      data ak12cs( 30) / +.3145111321 3578486743 0399999999 9 d-29     /
      data ak12cs( 31) / -.6092998312 1931276124 1600000000 0 d-30     /
      data ak12cs( 32) / +.1202021939 3698158346 2399999999 9 d-30     /
      data ak12cs( 33) / -.2412930801 4594088413 8666666666 6 d-31     /
c
      data ntk1, ntak1, ntak12, xmin, xsml / 3*0, 2*0.d0 /
c
      if (ntk1.ne.0) go to 10
      eta = 0.1*sngl(d1mach(3))
      ntk1 = initds (bk1cs, 16, eta)
      ntak1 = initds (ak1cs, 38, eta)
      ntak12 = initds (ak12cs, 33, eta)
c
      xmin = dexp (dmax1(dlog(d1mach(1)), -dlog(d1mach(2))) + 0.01d0)
      xsml = dsqrt (4.0d0*d1mach(3))
c
 10   if (x.le.0.d0) call seteru (29hdbsk1e  x is zero or negative, 29,
     1  1, 2)
      if (x.gt.2.0d0) go to 20
c
      if (x.lt.xmin) call seteru (31hdbsk1e  x so small k1 overflows,
     1  31, 2, 2)
      y = 0.d0
      if (x.gt.xsml) y = x*x
      dbsk1e = dexp(x)*(dlog(0.5d0*x)*dbesi1(x) + (0.75d0 +
     1  dcsevl (0.5d0*y-1.d0, bk1cs, ntk1))/x )
      return
c
 20   if (x.le.8.d0) dbsk1e = (1.25d0 + dcsevl ((16.d0/x-5.d0)/3.d0,
     1  ak1cs, ntak1))/dsqrt(x)
      if (x.gt.8.d0) dbsk1e = (1.25d0 +
     1  dcsevl (16.d0/x-1.d0, ak12cs, ntak12))/dsqrt(x)
c
      return
      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)
c     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
c     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
      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 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