1 SUBROUTINE ZBESJ(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR) 2C***BEGIN PROLOGUE ZBESJ 3C***DATE WRITTEN 830501 (YYMMDD) 4C***REVISION DATE 890801 (YYMMDD) 5C***CATEGORY NO. B5K 6C***KEYWORDS J-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT, 7C BESSEL FUNCTION OF FIRST KIND 8C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES 9C***PURPOSE TO COMPUTE THE J-BESSEL FUNCTION OF A COMPLEX ARGUMENT 10C***DESCRIPTION 11C 12C ***A DOUBLE PRECISION ROUTINE*** 13C ON KODE=1, CBESJ COMPUTES AN N MEMBER SEQUENCE OF COMPLEX 14C BESSEL FUNCTIONS CY(I)=J(FNU+I-1,Z) FOR REAL, NONNEGATIVE 15C ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE 16C -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESJ RETURNS THE SCALED 17C FUNCTIONS 18C 19C CY(I)=EXP(-ABS(Y))*J(FNU+I-1,Z) I = 1,...,N , Y=AIMAG(Z) 20C 21C WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND 22C LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION 23C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS 24C (REF. 1). 25C 26C INPUT ZR,ZI,FNU ARE DOUBLE PRECISION 27C ZR,ZI - Z=CMPLX(ZR,ZI), -PI.LT.ARG(Z).LE.PI 28C FNU - ORDER OF INITIAL J FUNCTION, FNU.GE.0.0D0 29C KODE - A PARAMETER TO INDICATE THE SCALING OPTION 30C KODE= 1 RETURNS 31C CY(I)=J(FNU+I-1,Z), I=1,...,N 32C = 2 RETURNS 33C CY(I)=J(FNU+I-1,Z)EXP(-ABS(Y)), I=1,...,N 34C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 35C 36C OUTPUT CYR,CYI ARE DOUBLE PRECISION 37C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS 38C CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE 39C CY(I)=J(FNU+I-1,Z) OR 40C CY(I)=J(FNU+I-1,Z)EXP(-ABS(Y)) I=1,...,N 41C DEPENDING ON KODE, Y=AIMAG(Z). 42C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW, 43C NZ= 0 , NORMAL RETURN 44C NZ.GT.0 , LAST NZ COMPONENTS OF CY SET ZERO DUE 45C TO UNDERFLOW, CY(I)=CMPLX(0.0D0,0.0D0), 46C I = N-NZ+1,...,N 47C IERR - ERROR FLAG 48C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED 49C IERR=1, INPUT ERROR - NO COMPUTATION 50C IERR=2, OVERFLOW - NO COMPUTATION, AIMAG(Z) 51C TOO LARGE ON KODE=1 52C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE 53C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT 54C REDUCTION PRODUCE LESS THAN HALF OF MACHINE 55C ACCURACY 56C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- 57C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- 58C CANCE BY ARGUMENT REDUCTION 59C IERR=5, ERROR - NO COMPUTATION, 60C ALGORITHM TERMINATION CONDITION NOT MET 61C 62C***LONG DESCRIPTION 63C 64C THE COMPUTATION IS CARRIED OUT BY THE FORMULA 65C 66C J(FNU,Z)=EXP( FNU*PI*I/2)*I(FNU,-I*Z) AIMAG(Z).GE.0.0 67C 68C J(FNU,Z)=EXP(-FNU*PI*I/2)*I(FNU, I*Z) AIMAG(Z).LT.0.0 69C 70C WHERE I**2 = -1 AND I(FNU,Z) IS THE I BESSEL FUNCTION. 71C 72C FOR NEGATIVE ORDERS,THE FORMULA 73C 74C J(-FNU,Z) = J(FNU,Z)*COS(PI*FNU) - Y(FNU,Z)*SIN(PI*FNU) 75C 76C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO INTEGERS, THE 77C THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE POSITIVE 78C INTEGER,THE MAGNITUDE OF J(-FNU,Z)=J(FNU,Z)*COS(PI*FNU) IS A 79C LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS NOT AN INTEGER, 80C Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A LARGE POSITIVE POWER OF 81C TEN AND THE MOST THAT THE SECOND TERM CAN BE REDUCED IS BY 82C UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, WIDE CHANGES CAN 83C OCCUR WITHIN UNIT ROUNDOFF OF A LARGE INTEGER FOR FNU. HERE, 84C LARGE MEANS FNU.GT.CABS(Z). 85C 86C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- 87C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS 88C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. 89C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN 90C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG 91C IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS 92C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. 93C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS 94C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS 95C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE 96C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS 97C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 98C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION 99C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION 100C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN 101C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT 102C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS 103C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. 104C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. 105C 106C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX 107C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT 108C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- 109C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE 110C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), 111C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF 112C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY 113C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN 114C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY 115C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER 116C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, 117C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS 118C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER 119C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY 120C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER 121C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE 122C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, 123C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, 124C OR -PI/2+P. 125C 126C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ 127C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF 128C COMMERCE, 1955. 129C 130C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 131C BY D. E. AMOS, SAND83-0083, MAY, 1983. 132C 133C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 134C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 135C 136C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 137C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- 138C 1018, MAY, 1985 139C 140C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 141C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. 142C MATH. SOFTWARE, 1986 143C 144C***ROUTINES CALLED ZBINU,I1MACH,D1MACH 145C***END PROLOGUE ZBESJ 146C 147C COMPLEX CI,CSGN,CY,Z,ZN 148 DOUBLE PRECISION AA, ALIM, ARG, CII, CSGNI, CSGNR, CYI, CYR, DIG, 149 * ELIM, FNU, FNUL, HPI, RL, R1M5, STR, TOL, ZI, ZNI, ZNR, ZR, 150 * D1MACH, BB, FN, AZ, XZABS, ASCLE, RTOL, ATOL, STI 151 INTEGER I, IERR, INU, INUH, IR, K, KODE, K1, K2, N, NL, NZ, I1MACH 152 DIMENSION CYR(N), CYI(N) 153 DATA HPI /1.57079632679489662D0/ 154C 155C***FIRST EXECUTABLE STATEMENT ZBESJ 156 IERR = 0 157 NZ=0 158 IF (FNU.LT.0.0D0) IERR=1 159 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 160 IF (N.LT.1) IERR=1 161 IF (IERR.NE.0) RETURN 162C----------------------------------------------------------------------- 163C SET PARAMETERS RELATED TO MACHINE CONSTANTS. 164C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. 165C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. 166C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND 167C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR 168C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. 169C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. 170C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). 171C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU. 172C----------------------------------------------------------------------- 173 TOL = DMAX1(D1MACH(4),1.0D-18) 174 K1 = I1MACH(15) 175 K2 = I1MACH(16) 176 R1M5 = D1MACH(5) 177 K = MIN0(IABS(K1),IABS(K2)) 178 ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) 179 K1 = I1MACH(14) - 1 180 AA = R1M5*DBLE(FLOAT(K1)) 181 DIG = DMIN1(AA,18.0D0) 182 AA = AA*2.303D0 183 ALIM = ELIM + DMAX1(-AA,-41.45D0) 184 RL = 1.2D0*DIG + 3.0D0 185 FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) 186C----------------------------------------------------------------------- 187C TEST FOR PROPER RANGE 188C----------------------------------------------------------------------- 189 AZ = XZABS(ZR,ZI) 190 FN = FNU+DBLE(FLOAT(N-1)) 191 AA = 0.5D0/TOL 192 BB=DBLE(FLOAT(I1MACH(9)))*0.5D0 193 AA = DMIN1(AA,BB) 194 IF (AZ.GT.AA) GO TO 260 195 IF (FN.GT.AA) GO TO 260 196 AA = DSQRT(AA) 197 IF (AZ.GT.AA) IERR=3 198 IF (FN.GT.AA) IERR=3 199C----------------------------------------------------------------------- 200C CALCULATE CSGN=EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE 201C WHEN FNU IS LARGE 202C----------------------------------------------------------------------- 203 35 CONTINUE 204 CII = 1.0D0 205 INU = INT(SNGL(FNU)) 206 INUH = INU/2 207 IR = INU - 2*INUH 208 ARG = (FNU-DBLE(FLOAT(INU-IR)))*HPI 209 CSGNR = DCOS(ARG) 210 CSGNI = DSIN(ARG) 211 IF (MOD(INUH,2).EQ.0) GO TO 40 212 CSGNR = -CSGNR 213 CSGNI = -CSGNI 214 40 CONTINUE 215C----------------------------------------------------------------------- 216C ZN IS IN THE RIGHT HALF PLANE 217C----------------------------------------------------------------------- 218 ZNR = ZI 219 ZNI = -ZR 220 IF (ZI.GE.0.0D0) GO TO 50 221 ZNR = -ZNR 222 ZNI = -ZNI 223 CSGNI = -CSGNI 224 CII = -CII 225 50 CONTINUE 226 CALL ZBINU(ZNR, ZNI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL, TOL, 227 * ELIM, ALIM) 228 IF (NZ.LT.0) GO TO 130 229 NL = N - NZ 230 IF (NL.EQ.0) RETURN 231 RTOL = 1.0D0/TOL 232 ASCLE = D1MACH(1)*RTOL*1.0D+3 233 DO 60 I=1,NL 234C STR = CYR(I)*CSGNR - CYI(I)*CSGNI 235C CYI(I) = CYR(I)*CSGNI + CYI(I)*CSGNR 236C CYR(I) = STR 237 AA = CYR(I) 238 BB = CYI(I) 239 ATOL = 1.0D0 240 IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 55 241 AA = AA*RTOL 242 BB = BB*RTOL 243 ATOL = TOL 244 55 CONTINUE 245 STR = AA*CSGNR - BB*CSGNI 246 STI = AA*CSGNI + BB*CSGNR 247 CYR(I) = STR*ATOL 248 CYI(I) = STI*ATOL 249 STR = -CSGNI*CII 250 CSGNI = CSGNR*CII 251 CSGNR = STR 252 60 CONTINUE 253 RETURN 254 130 CONTINUE 255 IF(NZ.EQ.(-2)) GO TO 140 256 NZ = 0 257 IERR = 2 258 RETURN 259 140 CONTINUE 260 NZ=0 261 IERR=5 262 RETURN 263 260 CONTINUE 264 IERR=4 265 GO TO 35 266 RETURN 267 END 268