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