00001 SUBROUTINE ZBESK(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR) 00002 C***BEGIN PROLOGUE ZBESK 00003 C***DATE WRITTEN 830501 (YYMMDD) 00004 C***REVISION DATE 890801 (YYMMDD) 00005 C***CATEGORY NO. B5K 00006 C***KEYWORDS K-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION, 00007 C MODIFIED BESSEL FUNCTION OF THE SECOND KIND, 00008 C BESSEL FUNCTION OF THE THIRD KIND 00009 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES 00010 C***PURPOSE TO COMPUTE K-BESSEL FUNCTIONS OF COMPLEX ARGUMENT 00011 C***DESCRIPTION 00012 C 00013 C ***A DOUBLE PRECISION ROUTINE*** 00014 C 00015 C ON KODE=1, CBESK COMPUTES AN N MEMBER SEQUENCE OF COMPLEX 00016 C BESSEL FUNCTIONS CY(J)=K(FNU+J-1,Z) FOR REAL, NONNEGATIVE 00017 C ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z.NE.CMPLX(0.0,0.0) 00018 C IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESK 00019 C RETURNS THE SCALED K FUNCTIONS, 00020 C 00021 C CY(J)=EXP(Z)*K(FNU+J-1,Z) , J=1,...,N, 00022 C 00023 C WHICH REMOVE THE EXPONENTIAL BEHAVIOR IN BOTH THE LEFT AND 00024 C RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND 00025 C NOTATION ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL 00026 C FUNCTIONS (REF. 1). 00027 C 00028 C INPUT ZR,ZI,FNU ARE DOUBLE PRECISION 00029 C ZR,ZI - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0), 00030 C -PI.LT.ARG(Z).LE.PI 00031 C FNU - ORDER OF INITIAL K FUNCTION, FNU.GE.0.0D0 00032 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 00033 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION 00034 C KODE= 1 RETURNS 00035 C CY(I)=K(FNU+I-1,Z), I=1,...,N 00036 C = 2 RETURNS 00037 C CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N 00038 C 00039 C OUTPUT CYR,CYI ARE DOUBLE PRECISION 00040 C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS 00041 C CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE 00042 C CY(I)=K(FNU+I-1,Z), I=1,...,N OR 00043 C CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N 00044 C DEPENDING ON KODE 00045 C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW. 00046 C NZ= 0 , NORMAL RETURN 00047 C NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO DUE 00048 C TO UNDERFLOW, CY(I)=CMPLX(0.0D0,0.0D0), 00049 C I=1,...,N WHEN X.GE.0.0. WHEN X.LT.0.0 00050 C NZ STATES ONLY THE NUMBER OF UNDERFLOWS 00051 C IN THE SEQUENCE. 00052 C 00053 C IERR - ERROR FLAG 00054 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED 00055 C IERR=1, INPUT ERROR - NO COMPUTATION 00056 C IERR=2, OVERFLOW - NO COMPUTATION, FNU IS 00057 C TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH 00058 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE 00059 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT 00060 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE 00061 C ACCURACY 00062 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- 00063 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- 00064 C CANCE BY ARGUMENT REDUCTION 00065 C IERR=5, ERROR - NO COMPUTATION, 00066 C ALGORITHM TERMINATION CONDITION NOT MET 00067 C 00068 C***LONG DESCRIPTION 00069 C 00070 C EQUATIONS OF THE REFERENCE ARE IMPLEMENTED FOR SMALL ORDERS 00071 C DNU AND DNU+1.0 IN THE RIGHT HALF PLANE X.GE.0.0. FORWARD 00072 C RECURRENCE GENERATES HIGHER ORDERS. K IS CONTINUED TO THE LEFT 00073 C HALF PLANE BY THE RELATION 00074 C 00075 C K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z) 00076 C MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1 00077 C 00078 C WHERE I(FNU,Z) IS THE I BESSEL FUNCTION. 00079 C 00080 C FOR LARGE ORDERS, FNU.GT.FNUL, THE K FUNCTION IS COMPUTED 00081 C BY MEANS OF ITS UNIFORM ASYMPTOTIC EXPANSIONS. 00082 C 00083 C FOR NEGATIVE ORDERS, THE FORMULA 00084 C 00085 C K(-FNU,Z) = K(FNU,Z) 00086 C 00087 C CAN BE USED. 00088 C 00089 C CBESK ASSUMES THAT A SIGNIFICANT DIGIT SINH(X) FUNCTION IS 00090 C AVAILABLE. 00091 C 00092 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- 00093 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS 00094 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. 00095 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN 00096 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG 00097 C IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS 00098 C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. 00099 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS 00100 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS 00101 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE 00102 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS 00103 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 00104 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION 00105 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION 00106 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN 00107 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT 00108 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS 00109 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. 00110 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. 00111 C 00112 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX 00113 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT 00114 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- 00115 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE 00116 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), 00117 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF 00118 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY 00119 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN 00120 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY 00121 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER 00122 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, 00123 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS 00124 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER 00125 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY 00126 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER 00127 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE 00128 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, 00129 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, 00130 C OR -PI/2+P. 00131 C 00132 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ 00133 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF 00134 C COMMERCE, 1955. 00135 C 00136 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 00137 C BY D. E. AMOS, SAND83-0083, MAY, 1983. 00138 C 00139 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 00140 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983. 00141 C 00142 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 00143 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- 00144 C 1018, MAY, 1985 00145 C 00146 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 00147 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. 00148 C MATH. SOFTWARE, 1986 00149 C 00150 C***ROUTINES CALLED ZACON,ZBKNU,ZBUNK,ZUOIK,XZABS,I1MACH,D1MACH 00151 C***END PROLOGUE ZBESK 00152 C 00153 C COMPLEX CY,Z 00154 DOUBLE PRECISION AA, ALIM, ALN, ARG, AZ, CYI, CYR, DIG, ELIM, FN, 00155 * FNU, FNUL, RL, R1M5, TOL, UFL, ZI, ZR, D1MACH, XZABS, BB 00156 INTEGER IERR, K, KODE, K1, K2, MR, N, NN, NUF, NW, NZ, I1MACH 00157 DIMENSION CYR(N), CYI(N) 00158 C***FIRST EXECUTABLE STATEMENT ZBESK 00159 IERR = 0 00160 NZ=0 00161 IF (ZI.EQ.0.0E0 .AND. ZR.EQ.0.0E0) IERR=1 00162 IF (FNU.LT.0.0D0) IERR=1 00163 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 00164 IF (N.LT.1) IERR=1 00165 IF (IERR.NE.0) RETURN 00166 NN = N 00167 C----------------------------------------------------------------------- 00168 C SET PARAMETERS RELATED TO MACHINE CONSTANTS. 00169 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. 00170 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. 00171 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND 00172 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR 00173 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. 00174 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. 00175 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). 00176 C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU 00177 C----------------------------------------------------------------------- 00178 TOL = DMAX1(D1MACH(4),1.0D-18) 00179 K1 = I1MACH(15) 00180 K2 = I1MACH(16) 00181 R1M5 = D1MACH(5) 00182 K = MIN0(IABS(K1),IABS(K2)) 00183 ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) 00184 K1 = I1MACH(14) - 1 00185 AA = R1M5*DBLE(FLOAT(K1)) 00186 DIG = DMIN1(AA,18.0D0) 00187 AA = AA*2.303D0 00188 ALIM = ELIM + DMAX1(-AA,-41.45D0) 00189 FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) 00190 RL = 1.2D0*DIG + 3.0D0 00191 C----------------------------------------------------------------------------- 00192 C TEST FOR PROPER RANGE 00193 C----------------------------------------------------------------------- 00194 AZ = XZABS(ZR,ZI) 00195 FN = FNU + DBLE(FLOAT(NN-1)) 00196 AA = 0.5D0/TOL 00197 BB=DBLE(FLOAT(I1MACH(9)))*0.5D0 00198 AA = DMIN1(AA,BB) 00199 IF (AZ.GT.AA) GO TO 260 00200 IF (FN.GT.AA) GO TO 260 00201 AA = DSQRT(AA) 00202 IF (AZ.GT.AA) IERR=3 00203 IF (FN.GT.AA) IERR=3 00204 C----------------------------------------------------------------------- 00205 C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE 00206 C----------------------------------------------------------------------- 00207 C UFL = DEXP(-ELIM) 00208 UFL = D1MACH(1)*1.0D+3 00209 IF (AZ.LT.UFL) GO TO 180 00210 IF (FNU.GT.FNUL) GO TO 80 00211 IF (FN.LE.1.0D0) GO TO 60 00212 IF (FN.GT.2.0D0) GO TO 50 00213 IF (AZ.GT.TOL) GO TO 60 00214 ARG = 0.5D0*AZ 00215 ALN = -FN*DLOG(ARG) 00216 IF (ALN.GT.ELIM) GO TO 180 00217 GO TO 60 00218 50 CONTINUE 00219 CALL ZUOIK(ZR, ZI, FNU, KODE, 2, NN, CYR, CYI, NUF, TOL, ELIM, 00220 * ALIM) 00221 IF (NUF.LT.0) GO TO 180 00222 NZ = NZ + NUF 00223 NN = NN - NUF 00224 C----------------------------------------------------------------------- 00225 C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK 00226 C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I 00227 C----------------------------------------------------------------------- 00228 IF (NN.EQ.0) GO TO 100 00229 60 CONTINUE 00230 IF (ZR.LT.0.0D0) GO TO 70 00231 C----------------------------------------------------------------------- 00232 C RIGHT HALF PLANE COMPUTATION, REAL(Z).GE.0. 00233 C----------------------------------------------------------------------- 00234 CALL ZBKNU(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, TOL, ELIM, ALIM) 00235 IF (NW.LT.0) GO TO 200 00236 NZ=NW 00237 RETURN 00238 C----------------------------------------------------------------------- 00239 C LEFT HALF PLANE COMPUTATION 00240 C PI/2.LT.ARG(Z).LE.PI AND -PI.LT.ARG(Z).LT.-PI/2. 00241 C----------------------------------------------------------------------- 00242 70 CONTINUE 00243 IF (NZ.NE.0) GO TO 180 00244 MR = 1 00245 IF (ZI.LT.0.0D0) MR = -1 00246 CALL ZACON(ZR, ZI, FNU, KODE, MR, NN, CYR, CYI, NW, RL, FNUL, 00247 * TOL, ELIM, ALIM) 00248 IF (NW.LT.0) GO TO 200 00249 NZ=NW 00250 RETURN 00251 C----------------------------------------------------------------------- 00252 C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL 00253 C----------------------------------------------------------------------- 00254 80 CONTINUE 00255 MR = 0 00256 IF (ZR.GE.0.0D0) GO TO 90 00257 MR = 1 00258 IF (ZI.LT.0.0D0) MR = -1 00259 90 CONTINUE 00260 CALL ZBUNK(ZR, ZI, FNU, KODE, MR, NN, CYR, CYI, NW, TOL, ELIM, 00261 * ALIM) 00262 IF (NW.LT.0) GO TO 200 00263 NZ = NZ + NW 00264 RETURN 00265 100 CONTINUE 00266 IF (ZR.LT.0.0D0) GO TO 180 00267 RETURN 00268 180 CONTINUE 00269 NZ = 0 00270 IERR=2 00271 RETURN 00272 200 CONTINUE 00273 IF(NW.EQ.(-1)) GO TO 180 00274 NZ=0 00275 IERR=5 00276 RETURN 00277 260 CONTINUE 00278 NZ=0 00279 IERR=4 00280 RETURN 00281 END
1.7.1
