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, ZABS, 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 = ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0))) 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 CII = 1.0D0 204 INU = INT(SNGL(FNU)) 205 INUH = INU/2 206 IR = INU - 2*INUH 207 ARG = (FNU-DBLE(FLOAT(INU-IR)))*HPI 208 CSGNR = DCOS(ARG) 209 CSGNI = DSIN(ARG) 210 IF (MOD(INUH,2).EQ.0) GO TO 40 211 CSGNR = -CSGNR 212 CSGNI = -CSGNI 213 40 CONTINUE 214C----------------------------------------------------------------------- 215C ZN IS IN THE RIGHT HALF PLANE 216C----------------------------------------------------------------------- 217 ZNR = ZI 218 ZNI = -ZR 219 IF (ZI.GE.0.0D0) GO TO 50 220 ZNR = -ZNR 221 ZNI = -ZNI 222 CSGNI = -CSGNI 223 CII = -CII 224 50 CONTINUE 225 CALL ZBINU(ZNR, ZNI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL, TOL, 226 * ELIM, ALIM) 227 IF (NZ.LT.0) GO TO 130 228 NL = N - NZ 229 IF (NL.EQ.0) RETURN 230 RTOL = 1.0D0/TOL 231 ASCLE = D1MACH(1)*RTOL*1.0D+3 232 DO 60 I=1,NL 233C STR = CYR(I)*CSGNR - CYI(I)*CSGNI 234C CYI(I) = CYR(I)*CSGNI + CYI(I)*CSGNR 235C CYR(I) = STR 236 AA = CYR(I) 237 BB = CYI(I) 238 ATOL = 1.0D0 239 IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 55 240 AA = AA*RTOL 241 BB = BB*RTOL 242 ATOL = TOL 243 55 CONTINUE 244 STR = AA*CSGNR - BB*CSGNI 245 STI = AA*CSGNI + BB*CSGNR 246 CYR(I) = STR*ATOL 247 CYI(I) = STI*ATOL 248 STR = -CSGNI*CII 249 CSGNI = CSGNR*CII 250 CSGNR = STR 251 60 CONTINUE 252 RETURN 253 130 CONTINUE 254 IF(NZ.EQ.(-2)) GO TO 140 255 NZ = 0 256 IERR = 2 257 RETURN 258 140 CONTINUE 259 NZ=0 260 IERR=5 261 RETURN 262 260 CONTINUE 263 NZ=0 264 IERR=4 265 RETURN 266 END 267