00001 SUBROUTINE CBESI(Z, FNU, KODE, N, CY, NZ, IERR) 00002 C***BEGIN PROLOGUE CBESI 00003 C***DATE WRITTEN 830501 (YYMMDD) 00004 C***REVISION DATE 890801 (YYMMDD) 00005 C***CATEGORY NO. B5K 00006 C***KEYWORDS I-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION, 00007 C MODIFIED BESSEL FUNCTION OF THE FIRST KIND 00008 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES 00009 C***PURPOSE TO COMPUTE I-BESSEL FUNCTIONS OF COMPLEX ARGUMENT 00010 C***DESCRIPTION 00011 C 00012 C ON KODE=1, CBESI COMPUTES AN N MEMBER SEQUENCE OF COMPLEX 00013 C BESSEL FUNCTIONS CY(J)=I(FNU+J-1,Z) FOR REAL, NONNEGATIVE 00014 C ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z IN THE CUT PLANE 00015 C -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESI RETURNS THE SCALED 00016 C FUNCTIONS 00017 C 00018 C CY(J)=EXP(-ABS(X))*I(FNU+J-1,Z) J = 1,...,N , X=REAL(Z) 00019 C 00020 C WITH THE EXPONENTIAL GROWTH REMOVED IN BOTH THE LEFT AND 00021 C RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND 00022 C NOTATION ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL 00023 C FUNCTIONS (REF.1) 00024 C 00025 C INPUT 00026 C Z - Z=CMPLX(X,Y), -PI.LT.ARG(Z).LE.PI 00027 C FNU - ORDER OF INITIAL I FUNCTION, FNU.GE.0.0E0 00028 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION 00029 C KODE= 1 RETURNS 00030 C CY(J)=I(FNU+J-1,Z), J=1,...,N 00031 C = 2 RETURNS 00032 C CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X)), J=1,...,N 00033 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 00034 C 00035 C OUTPUT 00036 C CY - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN 00037 C VALUES FOR THE SEQUENCE 00038 C CY(J)=I(FNU+J-1,Z) OR 00039 C CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X)) J=1,...,N 00040 C DEPENDING ON KODE, X=REAL(Z) 00041 C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW, 00042 C NZ= 0 , NORMAL RETURN 00043 C NZ.GT.0 , LAST NZ COMPONENTS OF CY SET TO ZERO 00044 C DUE TO UNDERFLOW, CY(J)=CMPLX(0.0,0.0), 00045 C J = N-NZ+1,...,N 00046 C IERR - ERROR FLAG 00047 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED 00048 C IERR=1, INPUT ERROR - NO COMPUTATION 00049 C IERR=2, OVERFLOW - NO COMPUTATION, REAL(Z) TOO 00050 C LARGE ON KODE=1 00051 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE 00052 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT 00053 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE 00054 C ACCURACY 00055 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- 00056 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- 00057 C CANCE BY ARGUMENT REDUCTION 00058 C IERR=5, ERROR - NO COMPUTATION, 00059 C ALGORITHM TERMINATION CONDITION NOT MET 00060 C 00061 C***LONG DESCRIPTION 00062 C 00063 C THE COMPUTATION IS CARRIED OUT BY THE POWER SERIES FOR 00064 C SMALL CABS(Z), THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z), 00065 C THE MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN AND A 00066 C NEUMANN SERIES FOR IMTERMEDIATE MAGNITUDES, AND THE 00067 C UNIFORM ASYMPTOTIC EXPANSIONS FOR I(FNU,Z) AND J(FNU,Z) 00068 C FOR LARGE ORDERS. BACKWARD RECURRENCE IS USED TO GENERATE 00069 C SEQUENCES OR REDUCE ORDERS WHEN NECESSARY. 00070 C 00071 C THE CALCULATIONS ABOVE ARE DONE IN THE RIGHT HALF PLANE AND 00072 C CONTINUED INTO THE LEFT HALF PLANE BY THE FORMULA 00073 C 00074 C I(FNU,Z*EXP(M*PI)) = EXP(M*PI*FNU)*I(FNU,Z) REAL(Z).GT.0.0 00075 C M = +I OR -I, I**2=-1 00076 C 00077 C FOR NEGATIVE ORDERS,THE FORMULA 00078 C 00079 C I(-FNU,Z) = I(FNU,Z) + (2/PI)*SIN(PI*FNU)*K(FNU,Z) 00080 C 00081 C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO INTEGERS, THE 00082 C THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE POSITIVE 00083 C INTEGER,THE MAGNITUDE OF I(-FNU,Z)=I(FNU,Z) IS A LARGE 00084 C NEGATIVE POWER OF TEN. BUT WHEN FNU IS NOT AN INTEGER, 00085 C K(FNU,Z) DOMINATES IN MAGNITUDE WITH A LARGE POSITIVE POWER OF 00086 C TEN AND THE MOST THAT THE SECOND TERM CAN BE REDUCED IS BY 00087 C UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, WIDE CHANGES CAN 00088 C OCCUR WITHIN UNIT ROUNDOFF OF A LARGE INTEGER FOR FNU. HERE, 00089 C LARGE MEANS FNU.GT.CABS(Z). 00090 C 00091 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- 00092 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS 00093 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. 00094 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN 00095 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG 00096 C IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO 00097 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS 00098 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS 00099 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE 00100 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS 00101 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 00102 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION 00103 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION 00104 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN 00105 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT 00106 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS 00107 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. 00108 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. 00109 C 00110 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX 00111 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT 00112 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- 00113 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE 00114 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), 00115 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF 00116 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY 00117 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN 00118 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY 00119 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER 00120 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, 00121 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS 00122 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER 00123 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY 00124 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER 00125 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE 00126 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, 00127 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, 00128 C OR -PI/2+P. 00129 C 00130 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ 00131 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF 00132 C COMMERCE, 1955. 00133 C 00134 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 00135 C BY D. E. AMOS, SAND83-0083, MAY, 1983. 00136 C 00137 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 00138 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 00139 C 00140 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 00141 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- 00142 C 1018, MAY, 1985 00143 C 00144 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 00145 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. 00146 C MATH. SOFTWARE, 1986 00147 C 00148 C***ROUTINES CALLED CBINU,I1MACH,R1MACH 00149 C***END PROLOGUE CBESI 00150 COMPLEX CONE, CSGN, CY, Z, ZN 00151 REAL AA, ALIM, ARG, DIG, ELIM, FNU, FNUL, PI, RL, R1M5, S1, S2, 00152 * TOL, XX, YY, R1MACH, AZ, FN, BB, ASCLE, RTOL, ATOL 00153 INTEGER I, IERR, INU, K, KODE, K1, K2, N, NN, NZ, I1MACH 00154 DIMENSION CY(N) 00155 DATA PI /3.14159265358979324E0/ 00156 DATA CONE / (1.0E0,0.0E0) / 00157 C 00158 C***FIRST EXECUTABLE STATEMENT CBESI 00159 IERR = 0 00160 NZ=0 00161 IF (FNU.LT.0.0E0) IERR=1 00162 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 00163 IF (N.LT.1) IERR=1 00164 IF (IERR.NE.0) RETURN 00165 XX = REAL(Z) 00166 YY = AIMAG(Z) 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 = AMAX1(R1MACH(4),1.0E-18) 00179 K1 = I1MACH(12) 00180 K2 = I1MACH(13) 00181 R1M5 = R1MACH(5) 00182 K = MIN0(IABS(K1),IABS(K2)) 00183 ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) 00184 K1 = I1MACH(11) - 1 00185 AA = R1M5*FLOAT(K1) 00186 DIG = AMIN1(AA,18.0E0) 00187 AA = AA*2.303E0 00188 ALIM = ELIM + AMAX1(-AA,-41.45E0) 00189 RL = 1.2E0*DIG + 3.0E0 00190 FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) 00191 AZ = CABS(Z) 00192 C----------------------------------------------------------------------- 00193 C TEST FOR RANGE 00194 C----------------------------------------------------------------------- 00195 AA = 0.5E0/TOL 00196 BB=FLOAT(I1MACH(9))*0.5E0 00197 AA=AMIN1(AA,BB) 00198 IF(AZ.GT.AA) GO TO 140 00199 FN=FNU+FLOAT(N-1) 00200 IF(FN.GT.AA) GO TO 140 00201 AA=SQRT(AA) 00202 IF(AZ.GT.AA) IERR=3 00203 IF(FN.GT.AA) IERR=3 00204 ZN = Z 00205 CSGN = CONE 00206 IF (XX.GE.0.0E0) GO TO 40 00207 ZN = -Z 00208 C----------------------------------------------------------------------- 00209 C CALCULATE CSGN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE 00210 C WHEN FNU IS LARGE 00211 C----------------------------------------------------------------------- 00212 INU = INT(FNU) 00213 ARG = (FNU-FLOAT(INU))*PI 00214 IF (YY.LT.0.0E0) ARG = -ARG 00215 S1 = COS(ARG) 00216 S2 = SIN(ARG) 00217 CSGN = CMPLX(S1,S2) 00218 IF (MOD(INU,2).EQ.1) CSGN = -CSGN 00219 40 CONTINUE 00220 CALL CBINU(ZN, FNU, KODE, N, CY, NZ, RL, FNUL, TOL, ELIM, ALIM) 00221 IF (NZ.LT.0) GO TO 120 00222 IF (XX.GE.0.0E0) RETURN 00223 C----------------------------------------------------------------------- 00224 C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE 00225 C----------------------------------------------------------------------- 00226 NN = N - NZ 00227 IF (NN.EQ.0) RETURN 00228 RTOL = 1.0E0/TOL 00229 ASCLE = R1MACH(1)*RTOL*1.0E+3 00230 DO 50 I=1,NN 00231 C CY(I) = CY(I)*CSGN 00232 ZN=CY(I) 00233 AA=REAL(ZN) 00234 BB=AIMAG(ZN) 00235 ATOL=1.0E0 00236 IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 55 00237 ZN = ZN*CMPLX(RTOL,0.0E0) 00238 ATOL = TOL 00239 55 CONTINUE 00240 ZN = ZN*CSGN 00241 CY(I) = ZN*CMPLX(ATOL,0.0E0) 00242 CSGN = -CSGN 00243 50 CONTINUE 00244 RETURN 00245 120 CONTINUE 00246 IF(NZ.EQ.(-2)) GO TO 130 00247 NZ = 0 00248 IERR=2 00249 RETURN 00250 130 CONTINUE 00251 NZ=0 00252 IERR=5 00253 RETURN 00254 140 CONTINUE 00255 NZ=0 00256 IERR=4 00257 RETURN 00258 END
1.7.1
