amos/zbesy.f

   1       SUBROUTINE ZBESY(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, CWRKR, CWRKI,
   2      *                 IERR)
   3 C***BEGIN PROLOGUE  ZBESY
   4 C***DATE WRITTEN   830501   (YYMMDD)
   5 C***REVISION DATE  890801   (YYMMDD)
   6 C***CATEGORY NO.  B5K
   7 C***KEYWORDS  Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT,
   8 C             BESSEL FUNCTION OF SECOND KIND
   9 C***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
  10 C***PURPOSE  TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT
  11 C***DESCRIPTION
  12 C
  13 C                      ***A DOUBLE PRECISION ROUTINE***
  14 C
  15 C         ON KODE=1, CBESY COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
  16 C         BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE
  17 C         ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE
  18 C         -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESY RETURNS THE SCALED
  19 C         FUNCTIONS
  20 C
  21 C         CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z)   I = 1,...,N , Y=AIMAG(Z)
  22 C
  23 C         WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND
  24 C         LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION
  25 C         ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS
  26 C         (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 Y FUNCTION, FNU.GE.0.0D0
  32 C           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
  33 C                    KODE= 1  RETURNS
  34 C                             CY(I)=Y(FNU+I-1,Z), I=1,...,N
  35 C                        = 2  RETURNS
  36 C                             CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N
  37 C                             WHERE Y=AIMAG(Z)
  38 C           N      - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
  39 C           CWRKR, - DOUBLE PRECISION WORK VECTORS OF DIMENSION AT
  40 C           CWRKI    AT LEAST N
  41 C
  42 C         OUTPUT     CYR,CYI ARE DOUBLE PRECISION
  43 C           CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
  44 C                    CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
  45 C                    CY(I)=Y(FNU+I-1,Z)  OR
  46 C                    CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y))  I=1,...,N
  47 C                    DEPENDING ON KODE.
  48 C           NZ     - NZ=0 , A NORMAL RETURN
  49 C                    NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO
  50 C                    UNDERFLOW (GENERALLY ON KODE=2)
  51 C           IERR   - ERROR FLAG
  52 C                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
  53 C                    IERR=1, INPUT ERROR   - NO COMPUTATION
  54 C                    IERR=2, OVERFLOW      - NO COMPUTATION, FNU IS
  55 C                            TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH
  56 C                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
  57 C                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
  58 C                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
  59 C                            ACCURACY
  60 C                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
  61 C                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
  62 C                            CANCE BY ARGUMENT REDUCTION
  63 C                    IERR=5, ERROR              - NO COMPUTATION,
  64 C                            ALGORITHM TERMINATION CONDITION NOT MET
  65 C
  66 C***LONG DESCRIPTION
  67 C
  68 C         THE COMPUTATION IS CARRIED OUT BY THE FORMULA
  69 C
  70 C         Y(FNU,Z)=0.5*(H(1,FNU,Z)-H(2,FNU,Z))/I
  71 C
  72 C         WHERE I**2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z)
  73 C         AND H(2,FNU,Z) ARE CALCULATED IN CBESH.
  74 C
  75 C         FOR NEGATIVE ORDERS,THE FORMULA
  76 C
  77 C              Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU)
  78 C
  79 C         CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD
  80 C         INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE
  81 C         POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)*
  82 C         SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS
  83 C         NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A
  84 C         LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM
  85 C         CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS,
  86 C         WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF
  87 C         ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z).
  88 C
  89 C         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
  90 C         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
  91 C         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
  92 C         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
  93 C         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
  94 C         IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
  95 C         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
  96 C         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
  97 C         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
  98 C         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
  99 C         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
 100 C         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
 101 C         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
 102 C         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
 103 C         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
 104 C         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
 105 C         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
 106 C         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
 107 C         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
 108 C
 109 C         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
 110 C         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
 111 C         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
 112 C         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
 113 C         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
 114 C         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
 115 C         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
 116 C         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
 117 C         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
 118 C         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
 119 C         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
 120 C         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
 121 C         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
 122 C         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
 123 C         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
 124 C         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
 125 C         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
 126 C         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
 127 C         OR -PI/2+P.
 128 C
 129 C***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
 130 C                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
 131 C                 COMMERCE, 1955.
 132 C
 133 C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
 134 C                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
 135 C
 136 C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
 137 C                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
 138 C
 139 C               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
 140 C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
 141 C                 1018, MAY, 1985
 142 C
 143 C               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
 144 C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
 145 C                 MATH. SOFTWARE, 1986
 146 C
 147 C***ROUTINES CALLED  ZBESH,I1MACH,D1MACH
 148 C***END PROLOGUE  ZBESY
 149 C
 150 C     COMPLEX CWRK,CY,C1,C2,EX,HCI,Z,ZU,ZV
 151       DOUBLE PRECISION CWRKI, CWRKR, CYI, CYR, C1I, C1R, C2I, C2R,
 152      * ELIM, EXI, EXR, EY, FNU, HCII, STI, STR, TAY, ZI, ZR, DEXP,
 153      * D1MACH, ASCLE, RTOL, ATOL, AA, BB, TOL
 154       INTEGER I, IERR, K, KODE, K1, K2, N, NZ, NZ1, NZ2, I1MACH
 155       DIMENSION CYR(N), CYI(N), CWRKR(N), CWRKI(N)
 156 C***FIRST EXECUTABLE STATEMENT  ZBESY
 157       IERR = 0
 158       NZ=0
 159       IF (ZR.EQ.0.0D0 .AND. ZI.EQ.0.0D0) IERR=1
 160       IF (FNU.LT.0.0D0) IERR=1
 161       IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
 162       IF (N.LT.1) IERR=1
 163       IF (IERR.NE.0) RETURN
 164       HCII = 0.5D0
 165       CALL ZBESH(ZR, ZI, FNU, KODE, 1, N, CYR, CYI, NZ1, IERR)
 166       IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170
 167       CALL ZBESH(ZR, ZI, FNU, KODE, 2, N, CWRKR, CWRKI, NZ2, IERR)
 168       IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170
 169       NZ = MIN0(NZ1,NZ2)
 170       IF (KODE.EQ.2) GO TO 60
 171       DO 50 I=1,N
 172         STR = CWRKR(I) - CYR(I)
 173         STI = CWRKI(I) - CYI(I)
 174         CYR(I) = -STI*HCII
 175         CYI(I) = STR*HCII
 176    50 CONTINUE
 177       RETURN
 178    60 CONTINUE
 179       TOL = DMAX1(D1MACH(4),1.0D-18)
 180       K1 = I1MACH(15)
 181       K2 = I1MACH(16)
 182       K = MIN0(IABS(K1),IABS(K2))
 183       R1M5 = D1MACH(5)
 184 C-----------------------------------------------------------------------
 185 C     ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT
 186 C-----------------------------------------------------------------------
 187       ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0)
 188       EXR = DCOS(ZR)
 189       EXI = DSIN(ZR)
 190       EY = 0.0D0
 191       TAY = DABS(ZI+ZI)
 192       IF (TAY.LT.ELIM) EY = DEXP(-TAY)
 193       IF (ZI.LT.0.0D0) GO TO 90
 194       C1R = EXR*EY
 195       C1I = EXI*EY
 196       C2R = EXR
 197       C2I = -EXI
 198    70 CONTINUE
 199       NZ = 0
 200       RTOL = 1.0D0/TOL
 201       ASCLE = D1MACH(1)*RTOL*1.0D+3
 202       DO 80 I=1,N
 203 C       STR = C1R*CYR(I) - C1I*CYI(I)
 204 C       STI = C1R*CYI(I) + C1I*CYR(I)
 205 C       STR = -STR + C2R*CWRKR(I) - C2I*CWRKI(I)
 206 C       STI = -STI + C2R*CWRKI(I) + C2I*CWRKR(I)
 207 C       CYR(I) = -STI*HCII
 208 C       CYI(I) = STR*HCII
 209         AA = CWRKR(I)
 210         BB = CWRKI(I)
 211         ATOL = 1.0D0
 212         IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 75
 213           AA = AA*RTOL
 214           BB = BB*RTOL
 215           ATOL = TOL
 216    75   CONTINUE
 217         STR = (AA*C2R - BB*C2I)*ATOL
 218         STI = (AA*C2I + BB*C2R)*ATOL
 219         AA = CYR(I)
 220         BB = CYI(I)
 221         ATOL = 1.0D0
 222         IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 85
 223           AA = AA*RTOL
 224           BB = BB*RTOL
 225           ATOL = TOL
 226    85   CONTINUE
 227         STR = STR - (AA*C1R - BB*C1I)*ATOL
 228         STI = STI - (AA*C1I + BB*C1R)*ATOL
 229         CYR(I) = -STI*HCII
 230         CYI(I) =  STR*HCII
 231         IF (STR.EQ.0.0D0 .AND. STI.EQ.0.0D0 .AND. EY.EQ.0.0D0) NZ = NZ
 232      *   + 1
 233    80 CONTINUE
 234       RETURN
 235    90 CONTINUE
 236       C1R = EXR
 237       C1I = EXI
 238       C2R = EXR*EY
 239       C2I = -EXI*EY
 240       GO TO 70
 241   170 CONTINUE
 242       NZ = 0
 243       RETURN
 244       END