00001 SUBROUTINE CBESY(Z, FNU, KODE, N, CY, NZ, CWRK, IERR) 00002 C***BEGIN PROLOGUE CBESY 00003 C***DATE WRITTEN 830501 (YYMMDD) 00004 C***REVISION DATE 890801 (YYMMDD) 00005 C***CATEGORY NO. B5K 00006 C***KEYWORDS Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT, 00007 C BESSEL FUNCTION OF SECOND KIND 00008 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES 00009 C***PURPOSE TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT 00010 C***DESCRIPTION 00011 C 00012 C ON KODE=1, CBESY COMPUTES AN N MEMBER SEQUENCE OF COMPLEX 00013 C BESSEL FUNCTIONS CY(I)=Y(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, CBESY RETURNS THE SCALED 00016 C FUNCTIONS 00017 C 00018 C CY(I)=EXP(-ABS(Y))*Y(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), Z.NE.CMPLX(0.,0.),-PI.LT.ARG(Z).LE.PI 00027 C FNU - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0E0 00028 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION 00029 C KODE= 1 RETURNS 00030 C CY(I)=Y(FNU+I-1,Z), I=1,...,N 00031 C = 2 RETURNS 00032 C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N 00033 C WHERE Y=AIMAG(Z) 00034 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 00035 C CWRK - A COMPLEX WORK VECTOR OF DIMENSION AT LEAST N 00036 C 00037 C OUTPUT 00038 C CY - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN 00039 C VALUES FOR THE SEQUENCE 00040 C CY(I)=Y(FNU+I-1,Z) OR 00041 C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)) I=1,...,N 00042 C DEPENDING ON KODE. 00043 C NZ - NZ=0 , A NORMAL RETURN 00044 C NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO 00045 C UNDERFLOW (GENERALLY ON KODE=2) 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, FNU+N-1 IS 00050 C TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH 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 Y(FNU,Z)=0.5*(H(1,FNU,Z)-H(2,FNU,Z))/I 00066 C 00067 C WHERE I**2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z) 00068 C AND H(2,FNU,Z) ARE CALCULATED IN CBESH. 00069 C 00070 C FOR NEGATIVE ORDERS,THE FORMULA 00071 C 00072 C Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU) 00073 C 00074 C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD 00075 C INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE 00076 C POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)* 00077 C SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS 00078 C NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A 00079 C LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM 00080 C CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, 00081 C WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF 00082 C ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z). 00083 C 00084 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- 00085 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS 00086 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. 00087 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN 00088 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG 00089 C IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO 00090 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS 00091 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS 00092 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE 00093 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS 00094 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 00095 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION 00096 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION 00097 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN 00098 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT 00099 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS 00100 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. 00101 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. 00102 C 00103 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX 00104 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT 00105 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- 00106 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE 00107 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), 00108 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF 00109 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY 00110 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN 00111 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY 00112 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER 00113 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, 00114 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS 00115 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER 00116 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY 00117 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER 00118 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE 00119 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, 00120 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, 00121 C OR -PI/2+P. 00122 C 00123 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ 00124 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF 00125 C COMMERCE, 1955. 00126 C 00127 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 00128 C BY D. E. AMOS, SAND83-0083, MAY, 1983. 00129 C 00130 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 00131 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 00132 C 00133 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 00134 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- 00135 C 1018, MAY, 1985 00136 C 00137 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 00138 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. 00139 C MATH. SOFTWARE, 1986 00140 C 00141 C***ROUTINES CALLED CBESH,I1MACH,R1MACH 00142 C***END PROLOGUE CBESY 00143 C 00144 COMPLEX CWRK, CY, C1, C2, EX, HCI, Z, ZU, ZV 00145 REAL ELIM, EY, FNU, R1, R2, TAY, XX, YY, R1MACH, ASCLE, RTOL, 00146 * ATOL, AA, BB 00147 INTEGER I, IERR, K, KODE, K1, K2, N, NZ, NZ1, NZ2, I1MACH 00148 DIMENSION CY(N), CWRK(N) 00149 C***FIRST EXECUTABLE STATEMENT CBESY 00150 XX = REAL(Z) 00151 YY = AIMAG(Z) 00152 IERR = 0 00153 NZ=0 00154 IF (XX.EQ.0.0E0 .AND. YY.EQ.0.0E0) IERR=1 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 HCI = CMPLX(0.0E0,0.5E0) 00160 CALL CBESH(Z, FNU, KODE, 1, N, CY, NZ1, IERR) 00161 IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170 00162 CALL CBESH(Z, FNU, KODE, 2, N, CWRK, NZ2, IERR) 00163 IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170 00164 NZ = MIN0(NZ1,NZ2) 00165 IF (KODE.EQ.2) GO TO 60 00166 DO 50 I=1,N 00167 CY(I) = HCI*(CWRK(I)-CY(I)) 00168 50 CONTINUE 00169 RETURN 00170 60 CONTINUE 00171 TOL = AMAX1(R1MACH(4),1.0E-18) 00172 K1 = I1MACH(12) 00173 K2 = I1MACH(13) 00174 K = MIN0(IABS(K1),IABS(K2)) 00175 R1M5 = R1MACH(5) 00176 C----------------------------------------------------------------------- 00177 C ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT 00178 C----------------------------------------------------------------------- 00179 ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) 00180 R1 = COS(XX) 00181 R2 = SIN(XX) 00182 EX = CMPLX(R1,R2) 00183 EY = 0.0E0 00184 TAY = ABS(YY+YY) 00185 IF (TAY.LT.ELIM) EY = EXP(-TAY) 00186 IF (YY.LT.0.0E0) GO TO 90 00187 C1 = EX*CMPLX(EY,0.0E0) 00188 C2 = CONJG(EX) 00189 70 CONTINUE 00190 NZ = 0 00191 RTOL = 1.0E0/TOL 00192 ASCLE = R1MACH(1)*RTOL*1.0E+3 00193 DO 80 I=1,N 00194 C CY(I) = HCI*(C2*CWRK(I)-C1*CY(I)) 00195 ZV = CWRK(I) 00196 AA=REAL(ZV) 00197 BB=AIMAG(ZV) 00198 ATOL=1.0E0 00199 IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 75 00200 ZV = ZV*CMPLX(RTOL,0.0E0) 00201 ATOL = TOL 00202 75 CONTINUE 00203 ZV = ZV*C2*HCI 00204 ZV = ZV*CMPLX(ATOL,0.0E0) 00205 ZU=CY(I) 00206 AA=REAL(ZU) 00207 BB=AIMAG(ZU) 00208 ATOL=1.0E0 00209 IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 85 00210 ZU = ZU*CMPLX(RTOL,0.0E0) 00211 ATOL = TOL 00212 85 CONTINUE 00213 ZU = ZU*C1*HCI 00214 ZU = ZU*CMPLX(ATOL,0.0E0) 00215 CY(I) = ZV - ZU 00216 IF (CY(I).EQ.CMPLX(0.0E0,0.0E0) .AND. EY.EQ.0.0E0) NZ = NZ + 1 00217 80 CONTINUE 00218 RETURN 00219 90 CONTINUE 00220 C1 = EX 00221 C2 = CONJG(EX)*CMPLX(EY,0.0E0) 00222 GO TO 70 00223 170 CONTINUE 00224 NZ = 0 00225 RETURN 00226 END
1.7.1
