relibc_openlibm/slatec/dbkisr.f

*DECK DBKISR
      SUBROUTINE DBKISR (X, N, SUM, IERR)
C***BEGIN PROLOGUE  DBKISR
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DBSKIN
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (BKISR-S, DBKISR-D)
C***AUTHOR  Amos, D. E., (SNLA)
C***DESCRIPTION
C
C     DBKISR computes repeated integrals of the K0 Bessel function
C     by the series for N=0,1, and 2.
C
C***SEE ALSO  DBSKIN
C***ROUTINES CALLED  D1MACH, DPSIXN
C***REVISION HISTORY  (YYMMDD)
C   820601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C   910722  Updated AUTHOR section.  (ALS)
C***END PROLOGUE  DBKISR
      INTEGER I, IERR, K, KK, KKN, K1, N, NP
      DOUBLE PRECISION AK, ATOL, BK, C, FK, FN, HX, HXS, POL, PR, SUM,
     * TKP, TOL, TRM, X, XLN
      DOUBLE PRECISION DPSIXN, D1MACH
      DIMENSION C(2)
      SAVE C
C
      DATA C(1), C(2) /1.57079632679489662D+00,1.0D0/
C***FIRST EXECUTABLE STATEMENT  DBKISR
      IERR=0
      TOL = MAX(D1MACH(4),1.0D-18)
      IF (X.LT.TOL) GO TO 50
      PR = 1.0D0
      POL = 0.0D0
      IF (N.EQ.0) GO TO 20
      DO 10 I=1,N
        POL = -POL*X + C(I)
        PR = PR*X/I
   10 CONTINUE
   20 CONTINUE
      HX = X*0.5D0
      HXS = HX*HX
      XLN = LOG(HX)
      NP = N + 1
      TKP = 3.0D0
      FK = 2.0D0
      FN = N
      BK = 4.0D0
      AK = 2.0D0/((FN+1.0D0)*(FN+2.0D0))
      SUM = AK*(DPSIXN(N+3)-DPSIXN(3)+DPSIXN(2)-XLN)
      ATOL = SUM*TOL*0.75D0
      DO 30 K=2,20
        AK = AK*(HXS/BK)*((TKP+1.0D0)/(TKP+FN+1.0D0))*(TKP/(TKP+FN))
        K1 = K + 1
        KK = K1 + K
        KKN = KK + N
        TRM = (DPSIXN(K1)+DPSIXN(KKN)-DPSIXN(KK)-XLN)*AK
        SUM = SUM + TRM
        IF (ABS(TRM).LE.ATOL) GO TO 40
        TKP = TKP + 2.0D0
        BK = BK + TKP
        FK = FK + 1.0D0
   30 CONTINUE
      GO TO 80
   40 CONTINUE
      SUM = (SUM*HXS+DPSIXN(NP)-XLN)*PR
      IF (N.EQ.1) SUM = -SUM
      SUM = POL + SUM
      RETURN
C-----------------------------------------------------------------------
C     SMALL X CASE, X.LT.WORD TOLERANCE
C-----------------------------------------------------------------------
   50 CONTINUE
      IF (N.GT.0) GO TO 60
      HX = X*0.5D0
      SUM = DPSIXN(1) - LOG(HX)
      RETURN
   60 CONTINUE
      SUM = C(N)
      RETURN
   80 CONTINUE
      IERR=2
      RETURN
      END