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