1*DECK ZBESI 2 SUBROUTINE ZBESI (ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR) 3C***BEGIN PROLOGUE ZBESI 4C***PURPOSE Compute a sequence of the Bessel functions I(a,z) for 5C complex argument z and real nonnegative orders a=b,b+1, 6C b+2,... where b>0. A scaling option is available to 7C help avoid overflow. 8C***LIBRARY SLATEC 9C***CATEGORY C10B4 10C***TYPE COMPLEX (CBESI-C, ZBESI-C) 11C***KEYWORDS BESSEL FUNCTIONS OF COMPLEX ARGUMENT, I BESSEL FUNCTIONS, 12C MODIFIED BESSEL FUNCTIONS 13C***AUTHOR Amos, D. E., (SNL) 14C***DESCRIPTION 15C 16C ***A DOUBLE PRECISION ROUTINE*** 17C On KODE=1, ZBESI computes an N-member sequence of complex 18C Bessel functions CY(L)=I(FNU+L-1,Z) for real nonnegative 19C orders FNU+L-1, L=1,...,N and complex Z in the cut plane 20C -pi<arg(Z)<=pi where Z=ZR+i*ZI. On KODE=2, CBESI returns 21C the scaled functions 22C 23C CY(L) = exp(-abs(X))*I(FNU+L-1,Z), L=1,...,N and X=Re(Z) 24C 25C which removes the exponential growth in both the left and 26C right half-planes as Z goes to infinity. 27C 28C Input 29C ZR - DOUBLE PRECISION real part of argument Z 30C ZI - DOUBLE PRECISION imag part of argument Z 31C FNU - DOUBLE PRECISION initial order, FNU>=0 32C KODE - A parameter to indicate the scaling option 33C KODE=1 returns 34C CY(L)=I(FNU+L-1,Z), L=1,...,N 35C =2 returns 36C CY(L)=exp(-abs(X))*I(FNU+L-1,Z), L=1,...,N 37C where X=Re(Z) 38C N - Number of terms in the sequence, N>=1 39C 40C Output 41C CYR - DOUBLE PRECISION real part of result vector 42C CYI - DOUBLE PRECISION imag part of result vector 43C NZ - Number of underflows set to zero 44C NZ=0 Normal return 45C NZ>0 CY(L)=0, L=N-NZ+1,...,N 46C IERR - Error flag 47C IERR=0 Normal return - COMPUTATION COMPLETED 48C IERR=1 Input error - NO COMPUTATION 49C IERR=2 Overflow - NO COMPUTATION 50C (Re(Z) too large on KODE=1) 51C IERR=3 Precision warning - COMPUTATION COMPLETED 52C (Result has half precision or less 53C because abs(Z) or FNU+N-1 is large) 54C IERR=4 Precision error - NO COMPUTATION 55C (Result has no precision because 56C abs(Z) or FNU+N-1 is too large) 57C IERR=5 Algorithmic error - NO COMPUTATION 58C (Termination condition not met) 59C 60C *Long Description: 61C 62C The computation of I(a,z) is carried out by the power series 63C for small abs(z), the asymptotic expansion for large abs(z), 64C the Miller algorithm normalized by the Wronskian and a 65C Neumann series for intermediate magnitudes of z, and the 66C uniform asymptotic expansions for I(a,z) and J(a,z) for 67C large orders a. Backward recurrence is used to generate 68C sequences or reduce orders when necessary. 69C 70C The calculations above are done in the right half plane and 71C continued into the left half plane by the formula 72C 73C I(a,z*exp(t)) = exp(t*a)*I(a,z), Re(z)>0 74C t = i*pi or -i*pi 75C 76C For negative orders, the formula 77C 78C I(-a,z) = I(a,z) + (2/pi)*sin(pi*a)*K(a,z) 79C 80C can be used. However, for large orders close to integers the 81C the function changes radically. When a is a large positive 82C integer, the magnitude of I(-a,z)=I(a,z) is a large 83C negative power of ten. But when a is not an integer, 84C K(a,z) dominates in magnitude with a large positive power of 85C ten and the most that the second term can be reduced is by 86C unit roundoff from the coefficient. Thus, wide changes can 87C occur within unit roundoff of a large integer for a. Here, 88C large means a>abs(z). 89C 90C In most complex variable computation, one must evaluate ele- 91C mentary functions. When the magnitude of Z or FNU+N-1 is 92C large, losses of significance by argument reduction occur. 93C Consequently, if either one exceeds U1=SQRT(0.5/UR), then 94C losses exceeding half precision are likely and an error flag 95C IERR=3 is triggered where UR=MAX(D1MACH(4),1.0D-18) is double 96C precision unit roundoff limited to 18 digits precision. Also, 97C if either is larger than U2=0.5/UR, then all significance is 98C lost and IERR=4. In order to use the INT function, arguments 99C must be further restricted not to exceed the largest machine 100C integer, U3=I1MACH(9). Thus, the magnitude of Z and FNU+N-1 101C is restricted by MIN(U2,U3). In IEEE arithmetic, U1,U2, and 102C U3 approximate 2.0E+3, 4.2E+6, 2.1E+9 in single precision 103C and 4.7E+7, 2.3E+15 and 2.1E+9 in double precision. This 104C makes U2 limiting in single precision and U3 limiting in 105C double precision. This means that one can expect to retain, 106C in the worst cases on IEEE machines, no digits in single pre- 107C cision and only 6 digits in double precision. Similar con- 108C siderations hold for other machines. 109C 110C The approximate relative error in the magnitude of a complex 111C Bessel function can be expressed as P*10**S where P=MAX(UNIT 112C ROUNDOFF,1.0E-18) is the nominal precision and 10**S repre- 113C sents the increase in error due to argument reduction in the 114C elementary functions. Here, S=MAX(1,ABS(LOG10(ABS(Z))), 115C ABS(LOG10(FNU))) approximately (i.e., S=MAX(1,ABS(EXPONENT OF 116C ABS(Z),ABS(EXPONENT OF FNU)) ). However, the phase angle may 117C have only absolute accuracy. This is most likely to occur 118C when one component (in magnitude) is larger than the other by 119C several orders of magnitude. If one component is 10**K larger 120C than the other, then one can expect only MAX(ABS(LOG10(P))-K, 121C 0) significant digits; or, stated another way, when K exceeds 122C the exponent of P, no significant digits remain in the smaller 123C component. However, the phase angle retains absolute accuracy 124C because, in complex arithmetic with precision P, the smaller 125C component will not (as a rule) decrease below P times the 126C magnitude of the larger component. In these extreme cases, 127C the principal phase angle is on the order of +P, -P, PI/2-P, 128C or -PI/2+P. 129C 130C***REFERENCES 1. M. Abramowitz and I. A. Stegun, Handbook of Mathe- 131C matical Functions, National Bureau of Standards 132C Applied Mathematics Series 55, U. S. Department 133C of Commerce, Tenth Printing (1972) or later. 134C 2. D. E. Amos, Computation of Bessel Functions of 135C Complex Argument, Report SAND83-0086, Sandia National 136C Laboratories, Albuquerque, NM, May 1983. 137C 3. D. E. Amos, Computation of Bessel Functions of 138C Complex Argument and Large Order, Report SAND83-0643, 139C Sandia National Laboratories, Albuquerque, NM, May 140C 1983. 141C 4. D. E. Amos, A Subroutine Package for Bessel Functions 142C of a Complex Argument and Nonnegative Order, Report 143C SAND85-1018, Sandia National Laboratory, Albuquerque, 144C NM, May 1985. 145C 5. D. E. Amos, A portable package for Bessel functions 146C of a complex argument and nonnegative order, ACM 147C Transactions on Mathematical Software, 12 (September 148C 1986), pp. 265-273. 149C 150C***ROUTINES CALLED D1MACH, I1MACH, ZABS, ZBINU 151C***REVISION HISTORY (YYMMDD) 152C 830501 DATE WRITTEN 153C 890801 REVISION DATE from Version 3.2 154C 910415 Prologue converted to Version 4.0 format. (BAB) 155C 920128 Category corrected. (WRB) 156C 920811 Prologue revised. (DWL) 157C***END PROLOGUE ZBESI 158C COMPLEX CONE,CSGN,CW,CY,CZERO,Z,ZN 159 DOUBLE PRECISION AA, ALIM, ARG, CONEI, CONER, CSGNI, CSGNR, CYI, 160 * CYR, DIG, ELIM, FNU, FNUL, PI, RL, R1M5, STR, TOL, ZI, ZNI, ZNR, 161 * ZR, D1MACH, AZ, BB, FN, ZABS, ASCLE, RTOL, ATOL, STI 162 INTEGER I, IERR, INU, K, KODE, K1,K2,N,NZ,NN, I1MACH 163 DIMENSION CYR(N), CYI(N) 164 EXTERNAL ZABS 165 DATA PI /3.14159265358979324D0/ 166 DATA CONER, CONEI /1.0D0,0.0D0/ 167C 168C***FIRST EXECUTABLE STATEMENT ZBESI 169 IERR = 0 170 NZ=0 171 IF (FNU.LT.0.0D0) IERR=1 172 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 173 IF (N.LT.1) IERR=1 174 IF (IERR.NE.0) RETURN 175C----------------------------------------------------------------------- 176C SET PARAMETERS RELATED TO MACHINE CONSTANTS. 177C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. 178C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. 179C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND 180C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR 181C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. 182C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. 183C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). 184C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU. 185C----------------------------------------------------------------------- 186 TOL = MAX(D1MACH(4),1.0D-18) 187 K1 = I1MACH(15) 188 K2 = I1MACH(16) 189 R1M5 = D1MACH(5) 190 K = MIN(ABS(K1),ABS(K2)) 191 ELIM = 2.303D0*(K*R1M5-3.0D0) 192 K1 = I1MACH(14) - 1 193 AA = R1M5*K1 194 DIG = MIN(AA,18.0D0) 195 AA = AA*2.303D0 196 ALIM = ELIM + MAX(-AA,-41.45D0) 197 RL = 1.2D0*DIG + 3.0D0 198 FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) 199C----------------------------------------------------------------------- 200C TEST FOR PROPER RANGE 201C----------------------------------------------------------------------- 202 AZ = ZABS(ZR,ZI) 203 FN = FNU+(N-1) 204 AA = 0.5D0/TOL 205 BB=I1MACH(9)*0.5D0 206 AA = MIN(AA,BB) 207 IF (AZ.GT.AA) GO TO 260 208 IF (FN.GT.AA) GO TO 260 209 AA = SQRT(AA) 210 IF (AZ.GT.AA) IERR=3 211 IF (FN.GT.AA) IERR=3 212 ZNR = ZR 213 ZNI = ZI 214 CSGNR = CONER 215 CSGNI = CONEI 216 IF (ZR.GE.0.0D0) GO TO 40 217 ZNR = -ZR 218 ZNI = -ZI 219C----------------------------------------------------------------------- 220C CALCULATE CSGN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE 221C WHEN FNU IS LARGE 222C----------------------------------------------------------------------- 223 INU = FNU 224 ARG = (FNU-INU)*PI 225 IF (ZI.LT.0.0D0) ARG = -ARG 226 CSGNR = COS(ARG) 227 CSGNI = SIN(ARG) 228 IF (MOD(INU,2).EQ.0) GO TO 40 229 CSGNR = -CSGNR 230 CSGNI = -CSGNI 231 40 CONTINUE 232 CALL ZBINU(ZNR, ZNI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL, TOL, 233 * ELIM, ALIM) 234 IF (NZ.LT.0) GO TO 120 235 IF (ZR.GE.0.0D0) RETURN 236C----------------------------------------------------------------------- 237C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE 238C----------------------------------------------------------------------- 239 NN = N - NZ 240 IF (NN.EQ.0) RETURN 241 RTOL = 1.0D0/TOL 242 ASCLE = D1MACH(1)*RTOL*1.0D+3 243 DO 50 I=1,NN 244C STR = CYR(I)*CSGNR - CYI(I)*CSGNI 245C CYI(I) = CYR(I)*CSGNI + CYI(I)*CSGNR 246C CYR(I) = STR 247 AA = CYR(I) 248 BB = CYI(I) 249 ATOL = 1.0D0 250 IF (MAX(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 55 251 AA = AA*RTOL 252 BB = BB*RTOL 253 ATOL = TOL 254 55 CONTINUE 255 STR = AA*CSGNR - BB*CSGNI 256 STI = AA*CSGNI + BB*CSGNR 257 CYR(I) = STR*ATOL 258 CYI(I) = STI*ATOL 259 CSGNR = -CSGNR 260 CSGNI = -CSGNI 261 50 CONTINUE 262 RETURN 263 120 CONTINUE 264 IF(NZ.EQ.(-2)) GO TO 130 265 NZ = 0 266 IERR=2 267 RETURN 268 130 CONTINUE 269 NZ=0 270 IERR=5 271 RETURN 272 260 CONTINUE 273 NZ=0 274 IERR=4 275 RETURN 276 END 277