amos/zbesk.f

   1       SUBROUTINE ZBESK(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR)
   2 C***BEGIN PROLOGUE  ZBESK
   3 C***DATE WRITTEN   830501   (YYMMDD)
   4 C***REVISION DATE  890801   (YYMMDD)
   5 C***CATEGORY NO.  B5K
   6 C***KEYWORDS  K-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION,
   7 C             MODIFIED BESSEL FUNCTION OF THE SECOND KIND,
   8 C             BESSEL FUNCTION OF THE THIRD KIND
   9 C***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
  10 C***PURPOSE  TO COMPUTE K-BESSEL FUNCTIONS OF COMPLEX ARGUMENT
  11 C***DESCRIPTION
  12 C
  13 C                      ***A DOUBLE PRECISION ROUTINE***
  14 C
  15 C         ON KODE=1, CBESK COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
  16 C         BESSEL FUNCTIONS CY(J)=K(FNU+J-1,Z) FOR REAL, NONNEGATIVE
  17 C         ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z.NE.CMPLX(0.0,0.0)
  18 C         IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESK
  19 C         RETURNS THE SCALED K FUNCTIONS,
  20 C
  21 C         CY(J)=EXP(Z)*K(FNU+J-1,Z) , J=1,...,N,
  22 C
  23 C         WHICH REMOVE THE EXPONENTIAL BEHAVIOR IN BOTH THE LEFT AND
  24 C         RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND
  25 C         NOTATION ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL
  26 C         FUNCTIONS (REF. 1).
  27 C
  28 C         INPUT      ZR,ZI,FNU ARE DOUBLE PRECISION
  29 C           ZR,ZI  - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0),
  30 C                    -PI.LT.ARG(Z).LE.PI
  31 C           FNU    - ORDER OF INITIAL K FUNCTION, FNU.GE.0.0D0
  32 C           N      - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
  33 C           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
  34 C                    KODE= 1  RETURNS
  35 C                             CY(I)=K(FNU+I-1,Z), I=1,...,N
  36 C                        = 2  RETURNS
  37 C                             CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N
  38 C
  39 C         OUTPUT     CYR,CYI ARE DOUBLE PRECISION
  40 C           CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
  41 C                    CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
  42 C                    CY(I)=K(FNU+I-1,Z), I=1,...,N OR
  43 C                    CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N
  44 C                    DEPENDING ON KODE
  45 C           NZ     - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW.
  46 C                    NZ= 0   , NORMAL RETURN
  47 C                    NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO DUE
  48 C                              TO UNDERFLOW, CY(I)=CMPLX(0.0D0,0.0D0),
  49 C                              I=1,...,N WHEN X.GE.0.0. WHEN X.LT.0.0
  50 C                              NZ STATES ONLY THE NUMBER OF UNDERFLOWS
  51 C                              IN THE SEQUENCE.
  52 C
  53 C           IERR   - ERROR FLAG
  54 C                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
  55 C                    IERR=1, INPUT ERROR   - NO COMPUTATION
  56 C                    IERR=2, OVERFLOW      - NO COMPUTATION, FNU IS
  57 C                            TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH
  58 C                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
  59 C                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
  60 C                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
  61 C                            ACCURACY
  62 C                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
  63 C                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
  64 C                            CANCE BY ARGUMENT REDUCTION
  65 C                    IERR=5, ERROR              - NO COMPUTATION,
  66 C                            ALGORITHM TERMINATION CONDITION NOT MET
  67 C
  68 C***LONG DESCRIPTION
  69 C
  70 C         EQUATIONS OF THE REFERENCE ARE IMPLEMENTED FOR SMALL ORDERS
  71 C         DNU AND DNU+1.0 IN THE RIGHT HALF PLANE X.GE.0.0. FORWARD
  72 C         RECURRENCE GENERATES HIGHER ORDERS. K IS CONTINUED TO THE LEFT
  73 C         HALF PLANE BY THE RELATION
  74 C
  75 C         K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z)
  76 C         MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1
  77 C
  78 C         WHERE I(FNU,Z) IS THE I BESSEL FUNCTION.
  79 C
  80 C         FOR LARGE ORDERS, FNU.GT.FNUL, THE K FUNCTION IS COMPUTED
  81 C         BY MEANS OF ITS UNIFORM ASYMPTOTIC EXPANSIONS.
  82 C
  83 C         FOR NEGATIVE ORDERS, THE FORMULA
  84 C
  85 C                       K(-FNU,Z) = K(FNU,Z)
  86 C
  87 C         CAN BE USED.
  88 C
  89 C         CBESK ASSUMES THAT A SIGNIFICANT DIGIT SINH(X) FUNCTION IS
  90 C         AVAILABLE.
  91 C
  92 C         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
  93 C         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
  94 C         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
  95 C         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
  96 C         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
  97 C         IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
  98 C         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
  99 C         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
 100 C         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
 101 C         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
 102 C         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
 103 C         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
 104 C         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
 105 C         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
 106 C         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
 107 C         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
 108 C         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
 109 C         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
 110 C         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
 111 C
 112 C         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
 113 C         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
 114 C         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
 115 C         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
 116 C         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
 117 C         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
 118 C         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
 119 C         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
 120 C         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
 121 C         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
 122 C         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
 123 C         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
 124 C         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
 125 C         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
 126 C         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
 127 C         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
 128 C         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
 129 C         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
 130 C         OR -PI/2+P.
 131 C
 132 C***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
 133 C                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
 134 C                 COMMERCE, 1955.
 135 C
 136 C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
 137 C                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
 138 C
 139 C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
 140 C                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983.
 141 C
 142 C               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
 143 C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
 144 C                 1018, MAY, 1985
 145 C
 146 C               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
 147 C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
 148 C                 MATH. SOFTWARE, 1986
 149 C
 150 C***ROUTINES CALLED  ZACON,ZBKNU,ZBUNK,ZUOIK,AZABS,I1MACH,D1MACH
 151 C***END PROLOGUE  ZBESK
 152 C
 153 C     COMPLEX CY,Z
 154       DOUBLE PRECISION AA, ALIM, ALN, ARG, AZ, CYI, CYR, DIG, ELIM, FN,
 155      * FNU, FNUL, RL, R1M5, TOL, UFL, ZI, ZR, D1MACH, AZABS, BB
 156       INTEGER IERR, K, KODE, K1, K2, MR, N, NN, NUF, NW, NZ, I1MACH
 157       DIMENSION CYR(N), CYI(N)
 158 C***FIRST EXECUTABLE STATEMENT  ZBESK
 159       IERR = 0
 160       NZ=0
 161       IF (ZI.EQ.0.0E0 .AND. ZR.EQ.0.0E0) IERR=1
 162       IF (FNU.LT.0.0D0) IERR=1
 163       IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
 164       IF (N.LT.1) IERR=1
 165       IF (IERR.NE.0) RETURN
 166       NN = N
 167 C-----------------------------------------------------------------------
 168 C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
 169 C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
 170 C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
 171 C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
 172 C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
 173 C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
 174 C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
 175 C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
 176 C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
 177 C-----------------------------------------------------------------------
 178       TOL = DMAX1(D1MACH(4),1.0D-18)
 179       K1 = I1MACH(15)
 180       K2 = I1MACH(16)
 181       R1M5 = D1MACH(5)
 182       K = MIN0(IABS(K1),IABS(K2))
 183       ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0)
 184       K1 = I1MACH(14) - 1
 185       AA = R1M5*DBLE(FLOAT(K1))
 186       DIG = DMIN1(AA,18.0D0)
 187       AA = AA*2.303D0
 188       ALIM = ELIM + DMAX1(-AA,-41.45D0)
 189       FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0)
 190       RL = 1.2D0*DIG + 3.0D0
 191 C-----------------------------------------------------------------------------
 192 C     TEST FOR PROPER RANGE
 193 C-----------------------------------------------------------------------
 194       AZ = AZABS(ZR,ZI)
 195       FN = FNU + DBLE(FLOAT(NN-1))
 196       AA = 0.5D0/TOL
 197       BB=DBLE(FLOAT(I1MACH(9)))*0.5D0
 198       AA = DMIN1(AA,BB)
 199       IF (AZ.GT.AA) GO TO 260
 200       IF (FN.GT.AA) GO TO 260
 201       AA = DSQRT(AA)
 202       IF (AZ.GT.AA) IERR=3
 203       IF (FN.GT.AA) IERR=3
 204 C-----------------------------------------------------------------------
 205 C     OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE
 206 C-----------------------------------------------------------------------
 207 C     UFL = DEXP(-ELIM)
 208       UFL = D1MACH(1)*1.0D+3
 209       IF (AZ.LT.UFL) GO TO 180
 210       IF (FNU.GT.FNUL) GO TO 80
 211       IF (FN.LE.1.0D0) GO TO 60
 212       IF (FN.GT.2.0D0) GO TO 50
 213       IF (AZ.GT.TOL) GO TO 60
 214       ARG = 0.5D0*AZ
 215       ALN = -FN*DLOG(ARG)
 216       IF (ALN.GT.ELIM) GO TO 180
 217       GO TO 60
 218    50 CONTINUE
 219       CALL ZUOIK(ZR, ZI, FNU, KODE, 2, NN, CYR, CYI, NUF, TOL, ELIM,
 220      * ALIM)
 221       IF (NUF.LT.0) GO TO 180
 222       NZ = NZ + NUF
 223       NN = NN - NUF
 224 C-----------------------------------------------------------------------
 225 C     HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK
 226 C     IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
 227 C-----------------------------------------------------------------------
 228       IF (NN.EQ.0) GO TO 100
 229    60 CONTINUE
 230       IF (ZR.LT.0.0D0) GO TO 70
 231 C-----------------------------------------------------------------------
 232 C     RIGHT HALF PLANE COMPUTATION, REAL(Z).GE.0.
 233 C-----------------------------------------------------------------------
 234       CALL ZBKNU(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, TOL, ELIM, ALIM)
 235       IF (NW.LT.0) GO TO 200
 236       NZ=NW
 237       RETURN
 238 C-----------------------------------------------------------------------
 239 C     LEFT HALF PLANE COMPUTATION
 240 C     PI/2.LT.ARG(Z).LE.PI AND -PI.LT.ARG(Z).LT.-PI/2.
 241 C-----------------------------------------------------------------------
 242    70 CONTINUE
 243       IF (NZ.NE.0) GO TO 180
 244       MR = 1
 245       IF (ZI.LT.0.0D0) MR = -1
 246       CALL ZACON(ZR, ZI, FNU, KODE, MR, NN, CYR, CYI, NW, RL, FNUL,
 247      * TOL, ELIM, ALIM)
 248       IF (NW.LT.0) GO TO 200
 249       NZ=NW
 250       RETURN
 251 C-----------------------------------------------------------------------
 252 C     UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL
 253 C-----------------------------------------------------------------------
 254    80 CONTINUE
 255       MR = 0
 256       IF (ZR.GE.0.0D0) GO TO 90
 257       MR = 1
 258       IF (ZI.LT.0.0D0) MR = -1
 259    90 CONTINUE
 260       CALL ZBUNK(ZR, ZI, FNU, KODE, MR, NN, CYR, CYI, NW, TOL, ELIM,
 261      * ALIM)
 262       IF (NW.LT.0) GO TO 200
 263       NZ = NZ + NW
 264       RETURN
 265   100 CONTINUE
 266       IF (ZR.LT.0.0D0) GO TO 180
 267       RETURN
 268   180 CONTINUE
 269       NZ = 0
 270       IERR=2
 271       RETURN
 272   200 CONTINUE
 273       IF(NW.EQ.(-1)) GO TO 180
 274       NZ=0
 275       IERR=5
 276       RETURN
 277   260 CONTINUE
 278       NZ=0
 279       IERR=4
 280       RETURN
 281       END