1 SUBROUTINE CBESY(Z, FNU, KODE, N, CY, NZ, CWRK, IERR) 2C***BEGIN PROLOGUE CBESY 3C***DATE WRITTEN 830501 (YYMMDD) 4C***REVISION DATE 890801 (YYMMDD) 5C***CATEGORY NO. B5K 6C***KEYWORDS Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT, 7C BESSEL FUNCTION OF SECOND KIND 8C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES 9C***PURPOSE TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT 10C***DESCRIPTION 11C 12C ON KODE=1, CBESY COMPUTES AN N MEMBER SEQUENCE OF COMPLEX 13C BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE 14C ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE 15C -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESY RETURNS THE SCALED 16C FUNCTIONS 17C 18C CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z) I = 1,...,N , Y=AIMAG(Z) 19C 20C WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND 21C LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION 22C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS 23C (REF. 1). 24C 25C INPUT 26C Z - Z=CMPLX(X,Y), Z.NE.CMPLX(0.,0.),-PI.LT.ARG(Z).LE.PI 27C FNU - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0E0 28C KODE - A PARAMETER TO INDICATE THE SCALING OPTION 29C KODE= 1 RETURNS 30C CY(I)=Y(FNU+I-1,Z), I=1,...,N 31C = 2 RETURNS 32C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N 33C WHERE Y=AIMAG(Z) 34C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 35C CWRK - A COMPLEX WORK VECTOR OF DIMENSION AT LEAST N 36C 37C OUTPUT 38C CY - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN 39C VALUES FOR THE SEQUENCE 40C CY(I)=Y(FNU+I-1,Z) OR 41C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)) I=1,...,N 42C DEPENDING ON KODE. 43C NZ - NZ=0 , A NORMAL RETURN 44C NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO 45C UNDERFLOW (GENERALLY ON KODE=2) 46C IERR - ERROR FLAG 47C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED 48C IERR=1, INPUT ERROR - NO COMPUTATION 49C IERR=2, OVERFLOW - NO COMPUTATION, FNU+N-1 IS 50C TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH 51C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE 52C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT 53C REDUCTION PRODUCE LESS THAN HALF OF MACHINE 54C ACCURACY 55C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- 56C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- 57C CANCE BY ARGUMENT REDUCTION 58C IERR=5, ERROR - NO COMPUTATION, 59C ALGORITHM TERMINATION CONDITION NOT MET 60C 61C***LONG DESCRIPTION 62C 63C THE COMPUTATION IS CARRIED OUT BY THE FORMULA 64C 65C Y(FNU,Z)=0.5*(H(1,FNU,Z)-H(2,FNU,Z))/I 66C 67C WHERE I**2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z) 68C AND H(2,FNU,Z) ARE CALCULATED IN CBESH. 69C 70C FOR NEGATIVE ORDERS,THE FORMULA 71C 72C Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU) 73C 74C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD 75C INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE 76C POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)* 77C SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS 78C NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A 79C LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM 80C CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, 81C WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF 82C ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z). 83C 84C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- 85C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS 86C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. 87C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN 88C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG 89C IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO 90C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS 91C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS 92C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE 93C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS 94C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 95C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION 96C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION 97C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN 98C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT 99C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS 100C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. 101C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. 102C 103C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX 104C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT 105C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- 106C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE 107C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), 108C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF 109C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY 110C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN 111C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY 112C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER 113C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, 114C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS 115C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER 116C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY 117C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER 118C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE 119C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, 120C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, 121C OR -PI/2+P. 122C 123C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ 124C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF 125C COMMERCE, 1955. 126C 127C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 128C BY D. E. AMOS, SAND83-0083, MAY, 1983. 129C 130C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 131C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 132C 133C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 134C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- 135C 1018, MAY, 1985 136C 137C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 138C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. 139C MATH. SOFTWARE, 1986 140C 141C***ROUTINES CALLED CBESH,I1MACH,R1MACH 142C***END PROLOGUE CBESY 143C 144 COMPLEX CWRK, CY, C1, C2, EX, HCI, Z, ZU, ZV 145 REAL ELIM, EY, FNU, R1, R2, TAY, XX, YY, R1MACH, ASCLE, RTOL, 146 * ATOL, AA, BB 147 INTEGER I, IERR, K, KODE, K1, K2, N, NZ, NZ1, NZ2, I1MACH 148 DIMENSION CY(N), CWRK(N) 149C***FIRST EXECUTABLE STATEMENT CBESY 150 XX = REAL(Z) 151 YY = AIMAG(Z) 152 IERR = 0 153 NZ=0 154 IF (XX.EQ.0.0E0 .AND. YY.EQ.0.0E0) IERR=1 155 IF (FNU.LT.0.0E0) IERR=1 156 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 157 IF (N.LT.1) IERR=1 158 IF (IERR.NE.0) RETURN 159 HCI = CMPLX(0.0E0,0.5E0) 160 CALL CBESH(Z, FNU, KODE, 1, N, CY, NZ1, IERR) 161 IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170 162 CALL CBESH(Z, FNU, KODE, 2, N, CWRK, NZ2, IERR) 163 IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170 164 NZ = MIN0(NZ1,NZ2) 165 IF (KODE.EQ.2) GO TO 60 166 DO 50 I=1,N 167 CY(I) = HCI*(CWRK(I)-CY(I)) 168 50 CONTINUE 169 RETURN 170 60 CONTINUE 171 TOL = AMAX1(R1MACH(4),1.0E-18) 172 K1 = I1MACH(12) 173 K2 = I1MACH(13) 174 K = MIN0(IABS(K1),IABS(K2)) 175 R1M5 = R1MACH(5) 176C----------------------------------------------------------------------- 177C ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT 178C----------------------------------------------------------------------- 179 ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) 180 R1 = COS(XX) 181 R2 = SIN(XX) 182 EX = CMPLX(R1,R2) 183 EY = 0.0E0 184 TAY = ABS(YY+YY) 185 IF (TAY.LT.ELIM) EY = EXP(-TAY) 186 IF (YY.LT.0.0E0) GO TO 90 187 C1 = EX*CMPLX(EY,0.0E0) 188 C2 = CONJG(EX) 189 70 CONTINUE 190 NZ = 0 191 RTOL = 1.0E0/TOL 192 ASCLE = R1MACH(1)*RTOL*1.0E+3 193 DO 80 I=1,N 194C CY(I) = HCI*(C2*CWRK(I)-C1*CY(I)) 195 ZV = CWRK(I) 196 AA=REAL(ZV) 197 BB=AIMAG(ZV) 198 ATOL=1.0E0 199 IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 75 200 ZV = ZV*CMPLX(RTOL,0.0E0) 201 ATOL = TOL 202 75 CONTINUE 203 ZV = ZV*C2*HCI 204 ZV = ZV*CMPLX(ATOL,0.0E0) 205 ZU=CY(I) 206 AA=REAL(ZU) 207 BB=AIMAG(ZU) 208 ATOL=1.0E0 209 IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 85 210 ZU = ZU*CMPLX(RTOL,0.0E0) 211 ATOL = TOL 212 85 CONTINUE 213 ZU = ZU*C1*HCI 214 ZU = ZU*CMPLX(ATOL,0.0E0) 215 CY(I) = ZV - ZU 216 IF (CY(I).EQ.CMPLX(0.0E0,0.0E0) .AND. EY.EQ.0.0E0) NZ = NZ + 1 217 80 CONTINUE 218 RETURN 219 90 CONTINUE 220 C1 = EX 221 C2 = CONJG(EX)*CMPLX(EY,0.0E0) 222 GO TO 70 223 170 CONTINUE 224 NZ = 0 225 RETURN 226 END 227