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