Logo Search packages:      
Sourcecode: maxima version File versions  Download package

dbesk.f

*DECK DBESK
      SUBROUTINE DBESK (X, FNU, KODE, N, Y, NZ)
C***BEGIN PROLOGUE  DBESK
C***PURPOSE  Implement forward recursion on the three term recursion
C            relation for a sequence of non-negative order Bessel
C            functions K/SUB(FNU+I-1)/(X), or scaled Bessel functions
C            EXP(X)*K/SUB(FNU+I-1)/(X), I=1,...,N for real, positive
C            X and non-negative orders FNU.
C***LIBRARY   SLATEC
C***CATEGORY  C10B3
C***TYPE      DOUBLE PRECISION (BESK-S, DBESK-D)
C***KEYWORDS  K BESSEL FUNCTION, SPECIAL FUNCTIONS
C***AUTHOR  Amos, D. E., (SNLA)
C***DESCRIPTION
C
C     Abstract  **** a double precision routine ****
C         DBESK implements forward recursion on the three term
C         recursion relation for a sequence of non-negative order Bessel
C         functions K/sub(FNU+I-1)/(X), or scaled Bessel functions
C         EXP(X)*K/sub(FNU+I-1)/(X), I=1,..,N for real X .GT. 0.0D0 and
C         non-negative orders FNU.  If FNU .LT. NULIM, orders FNU and
C         FNU+1 are obtained from DBSKNU to start the recursion.  If
C         FNU .GE. NULIM, the uniform asymptotic expansion is used for
C         orders FNU and FNU+1 to start the recursion.  NULIM is 35 or
C         70 depending on whether N=1 or N .GE. 2.  Under and overflow
C         tests are made on the leading term of the asymptotic expansion
C         before any extensive computation is done.
C
C         The maximum number of significant digits obtainable
C         is the smaller of 14 and the number of digits carried in
C         double precision arithmetic.
C
C     Description of Arguments
C
C         Input      X,FNU are double precision
C           X      - X .GT. 0.0D0
C           FNU    - order of the initial K function, FNU .GE. 0.0D0
C           KODE   - a parameter to indicate the scaling option
C                    KODE=1 returns Y(I)=       K/sub(FNU+I-1)/(X),
C                                        I=1,...,N
C                    KODE=2 returns Y(I)=EXP(X)*K/sub(FNU+I-1)/(X),
C                                        I=1,...,N
C           N      - number of members in the sequence, N .GE. 1
C
C         Output     Y is double precision
C           Y      - a vector whose first N components contain values
C                    for the sequence
C                    Y(I)=       k/sub(FNU+I-1)/(X), I=1,...,N  or
C                    Y(I)=EXP(X)*K/sub(FNU+I-1)/(X), I=1,...,N
C                    depending on KODE
C           NZ     - number of components of Y set to zero due to
C                    underflow with KODE=1,
C                    NZ=0   , normal return, computation completed
C                    NZ .NE. 0, first NZ components of Y set to zero
C                             due to underflow, Y(I)=0.0D0, I=1,...,NZ
C
C     Error Conditions
C         Improper input arguments - a fatal error
C         Overflow - a fatal error
C         Underflow with KODE=1 -  a non-fatal error (NZ .NE. 0)
C
C***REFERENCES  F. W. J. Olver, Tables of Bessel Functions of Moderate
C                 or Large Orders, NPL Mathematical Tables 6, Her
C                 Majesty





































































































































































































's Stationery Office, London, 1962.C               N. M. Temme, On the numerical evaluation of the modifiedC                 Bessel function of the third kind, Journal ofC                 Computational Physics 19, (1975), pp. 324-337.C***ROUTINES CALLED  D1MACH, DASYIK, DBESK0, DBESK1, DBSK0E, DBSK1E,C                    DBSKNU, I1MACH, XERMSGC***REVISION HISTORY  (YYMMDD)C   790201  DATE WRITTENC   890531  Changed all specific intrinsics to generic.  (WRB)C   890911  Removed unnecessary intrinsics.  (WRB)C   890911  REVISION DATE from Version 3.2C   891214  Prologue converted to Version 4.0 format.  (BAB)C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)C   920501  Reformatted the REFERENCES section.  (WRB)C***END PROLOGUE  DBESKC      INTEGER I, J, K, KODE, MZ, N, NB, ND, NN, NUD, NULIM, NZ      INTEGER I1MACH      DOUBLE PRECISION CN,DNU,ELIM,ETX,FLGIK,FN,FNN,FNU,GLN,GNU,RTZ,     1 S, S1, S2, T, TM, TRX, W, X, XLIM, Y, ZN      DOUBLE PRECISION DBESK0, DBESK1, DBSK1E, DBSK0E, D1MACH      DIMENSION W(2), NULIM(2), Y(*)      SAVE NULIM      DATA NULIM(1),NULIM(2) / 35 , 70 /C***FIRST EXECUTABLE STATEMENT  DBESK      NN = -I1MACH(15)      ELIM = 2.303D0*(NN*D1MACH(5)-3.0D0)      XLIM = D1MACH(1)*1.0D+3      IF (KODE.LT.1 .OR. KODE.GT.2) GO TO 280      IF (FNU.LT.0.0D0) GO TO 290      IF (X.LE.0.0D0) GO TO 300      IF (X.LT.XLIM) GO TO 320      IF (N.LT.1) GO TO 310      ETX = KODE - 1CC     ND IS A DUMMY VARIABLE FOR NC     GNU IS A DUMMY VARIABLE FOR FNUC     NZ = NUMBER OF UNDERFLOWS ON KODE=1C      ND = N      NZ = 0      NUD = INT(FNU)      DNU = FNU - NUD      GNU = FNU      NN = MIN(2,ND)      FN = FNU + N - 1      FNN = FN      IF (FN.LT.2.0D0) GO TO 150CC     OVERFLOW TEST  (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION)C     FOR THE LAST ORDER, FNU+N-1.GE.NULIMC      ZN = X/FN      IF (ZN.EQ.0.0D0) GO TO 320      RTZ = SQRT(1.0D0+ZN*ZN)      GLN = LOG((1.0D0+RTZ)/ZN)      T = RTZ*(1.0D0-ETX) + ETX/(ZN+RTZ)      CN = -FN*(T-GLN)      IF (CN.GT.ELIM) GO TO 320      IF (NUD.LT.NULIM(NN)) GO TO 30      IF (NN.EQ.1) GO TO 20   10 CONTINUECC     UNDERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION)C     FOR THE FIRST ORDER, FNU.GE.NULIMC      FN = GNU      ZN = X/FN      RTZ = SQRT(1.0D0+ZN*ZN)      GLN = LOG((1.0D0+RTZ)/ZN)      T = RTZ*(1.0D0-ETX) + ETX/(ZN+RTZ)      CN = -FN*(T-GLN)   20 CONTINUE      IF (CN.LT.-ELIM) GO TO 230CC     ASYMPTOTIC EXPANSION FOR ORDERS FNU AND FNU+1.GE.NULIMC      FLGIK = -1.0D0      CALL DASYIK(X,GNU,KODE,FLGIK,RTZ,CN,NN,Y)      IF (NN.EQ.1) GO TO 240      TRX = 2.0D0/X      TM = (GNU+GNU+2.0D0)/X      GO TO 130C   30 CONTINUE      IF (KODE.EQ.2) GO TO 40CC     UNDERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION IN X)C     FOR ORDER DNUC      IF (X.GT.ELIM) GO TO 230   40 CONTINUE      IF (DNU.NE.0.0D0) GO TO 80      IF (KODE.EQ.2) GO TO 50      S1 = DBESK0(X)      GO TO 60   50 S1 = DBSK0E(X)   60 CONTINUE      IF (NUD.EQ.0 .AND. ND.EQ.1) GO TO 120      IF (KODE.EQ.2) GO TO 70      S2 = DBESK1(X)      GO TO 90   70 S2 = DBSK1E(X)      GO TO 90   80 CONTINUE      NB = 2      IF (NUD.EQ.0 .AND. ND.EQ.1) NB = 1      CALL DBSKNU(X, DNU, KODE, NB, W, NZ)      S1 = W(1)      IF (NB.EQ.1) GO TO 120      S2 = W(2)   90 CONTINUE      TRX = 2.0D0/X      TM = (DNU+DNU+2.0D0)/XC     FORWARD RECUR FROM DNU TO FNU+1 TO GET Y(1) AND Y(2)      IF (ND.EQ.1) NUD = NUD - 1      IF (NUD.GT.0) GO TO 100      IF (ND.GT.1) GO TO 120      S1 = S2      GO TO 120  100 CONTINUE      DO 110 I=1,NUD        S = S2        S2 = TM*S2 + S1        S1 = S        TM = TM + TRX  110 CONTINUE      IF (ND.EQ.1) S1 = S2  120 CONTINUE      Y(1) = S1      IF (ND.EQ.1) GO TO 240      Y(2) = S2  130 CONTINUE      IF (ND.EQ.2) GO TO 240C     FORWARD RECUR FROM FNU+2 TO FNU+N-1      DO 140 I=3,ND        Y(I) = TM*Y(I-1) + Y(I-2)        TM = TM + TRX  140 CONTINUE      GO TO 240C  150 CONTINUEC     UNDERFLOW TEST FOR KODE=1      IF (KODE.EQ.2) GO TO 160      IF (X.GT.ELIM) GO TO 230  160 CONTINUEC     OVERFLOW TEST      IF (FN.LE.1.0D0) GO TO 170      IF (-FN*(LOG(X)-0.693D0).GT.ELIM) GO TO 320  170 CONTINUE      IF (DNU.EQ.0.0D0) GO TO 180      CALL DBSKNU(X, FNU, KODE, ND, Y, MZ)      GO TO 240  180 CONTINUE      J = NUD      IF (J.EQ.1) GO TO 210      J = J + 1      IF (KODE.EQ.2) GO TO 190      Y(J) = DBESK0(X)      GO TO 200  190 Y(J) = DBSK0E(X)  200 IF (ND.EQ.1) GO TO 240      J = J + 1  210 IF (KODE.EQ.2) GO TO 220      Y(J) = DBESK1(X)      GO TO 240  220 Y(J) = DBSK1E(X)      GO TO 240CC     UPDATE PARAMETERS ON UNDERFLOWC  230 CONTINUE      NUD = NUD + 1      ND = ND - 1      IF (ND.EQ.0) GO TO 240      NN = MIN(2,ND)      GNU = GNU + 1.0D0      IF (FNN.LT.2.0D0) GO TO 230      IF (NUD.LT.NULIM(NN)) GO TO 230      GO TO 10  240 CONTINUE      NZ = N - ND      IF (NZ.EQ.0) RETURN      IF (ND.EQ.0) GO TO 260      DO 250 I=1,ND        J = N - I + 1        K = ND - I + 1        Y(J) = Y(K)  250 CONTINUE  260 CONTINUE      DO 270 I=1,NZ        Y(I) = 0.0D0  270 CONTINUE      RETURNCCC  280 CONTINUE      CALL XERMSG ('SLATEC', 'DBESK
',     +   'SCALING OPTION, KODE, NOT 1 OR 2


', 2, 1)      RETURN  290 CONTINUE      CALL XERMSG ('SLATEC', 'DBESK', 'ORDER, FNU, LESS THAN ZERO



', 2,     +   1)      RETURN  300 CONTINUE      CALL XERMSG ('SLATEC', 'DBESK', 'X LESS THAN OR EQUAL TO ZERO



',     +   2, 1)      RETURN  310 CONTINUE      CALL XERMSG ('SLATEC', 'DBESK', 'N LESS THAN ONE


', 2, 1)      RETURN  320 CONTINUE      CALL XERMSG ('SLATEC', 'DBESK
',     +   'OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL



Generated by  Doxygen 1.6.0   Back to index