C REMARK ON ALGORITHM 644, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 21, NO. 4, December, 1995, P. 388--393. C C This file contains 6 files separated by lines of the form C C*** filename C C The filenames in this file are: C C Readme cbsubs.f cqcbes.f C machcon.f zbsubs.f zqcbes.f C C C*** Readme INSTALLATION HINTS FOR ALGORITHM 644 D. E. AMOS Sandia National Laboratories ----------------------------------------------------------------- Algorithm 644 [2],[3] computes all major Bessel functions of a complex argument and nonnegative order. The single precision callable routines are denoted by CBESH, CBESI, CBESJ, CBESK, CBESY, CAIRY AND CBIRY for the H, I, J, K, Y, and Airy functions Ai and Bi. The corresponding double precision routines start with Z in place of C. CBESH and ZBESH compute H functions of kinds 1 and 2 and an option exists for the derivatives of the Airy functions in CAIRY, ZAIRY, CBIRY, and ZBIRY. Scaling to remove exponential underflow or overflow and generation of sequences in increasing orders are auxilliary options where appropriate. It is convenient to consider the package in four parts: (1) Quick check PROGRAMS (CQC?? or ZQC??) which exercise the main callable routines in single and double precision arithmetic: CQCBH and ZQCBH to test CBESH and ZBESH CQCBI and ZQCBI to test CBESI and ZBESI CQCBJ and ZQCBJ to test CBESJ and ZBESJ CQCBK and ZQCBK to test CBESK and ZBESK CQCBY and ZQCBY to test CBESY and ZBESY CQCAI and ZQCAI to test CAIRY, ZAIRY, CBIRY and ZBIRY and auxilliary subroutines CBESYH and ZBESYH The documentation for the latest versions of these quick checks is found in reference [4]. These routines need not be a permanent part of the installation once they have been compiled and successfully executed. (2) The main callable subroutines which comprise the algorithm in single and double precision arithmetic: CBESH and ZBESH for H functions of kinds 1 and 2 CBESI and ZBESI for I functions CBESJ and ZBESJ for J functions CBESK and ZBESK for K functions CBESY and ZBESY for Y functions CAIRY and ZAIRY for Airy function Ai and Ai' CBIRY and ZBIRY for Airy function Bi and Bi' GAMLN and DGAMLN for lngamma(x), x>0 The Z routines manipulate double precision ordered pairs as double precision complex numbers and return double precison ordered pairs as answers. The lower level routines ZMLT,ZDIV,ZSQRT,ZEXP,ZLOG,ZSHCH,ZABS provide a means of doing double precision complex arithmetic on double precision ordered pairs. These can also be called by user-coded routines to accomplish double precision complex arithmetic. There are no ZADD or ZSUB routines -- double precision complex additions and subtractions are simple enough to do in-line. Each subroutine or function has a prologue which defines the calling sequence. These routines, along with those of (3) and (4), make up the permanent installation. (3) Included in the algorithm are three (callable) function subroutines which define the machine to the Bessel function package. These subroutines consist mainly of comment lines which contain (inactive) FORTRAN statements for machine constants such as real and double precision unit roundoff, real and double precision overflow and underflow limits, largest integer constant, base of integer, real and double precision arithmetic, etc. for a variety of machines. The prologue of each of these functions INTEGER I1MACH(I) I=1,16 for integer constants REAL R1MACH(I) I=1,5 for real constants DOUBLE PRECISION D1MACH(I) I=1,5 for double precision constants defines the return value for each of the calling indices. TO USE THE PACKAGE, YOU MUST DEFINE THE MACHINE ON WHICH THE PACKAGE IS TO BE USED BY EDITING EACH OF THESE FUNCTIONS. With an editor, search out your machine from among those listed in each ?1MACH function. A keyword like IBM, SUN, VAX, CDC or PC is often helpful in locating the right machine. When your machine has been located, replace the C in column 1 with a space to make each of the FORTRAN statements active for the selected machine. If your machine is not found, the correct constants can be computed from the definitions in the prologue. It often happens that many PC constants are close enough to the IBM or SUN constants to work with the package. Reference [5] should help in generating a set of machine constants. One does not need constants accurate to the last bit and leaving a little slack on arithmetic limits can sometimes avoid hardware or software problems. (4) Lower level non-callable routines The latest version containing all improvements is dated 930101. REFERENCES 1. Abramowitz, M. and Stegun, I. A., Handbook of Mathematical Functions, NBS Applied Math Series 55, U.S. Dept. of Commerce, Washington, D.C., 1955 2. Amos, D. E., Algorithm 644, A Portable Package For Bessel Functions of a Complex Argument and Nonnegative Order, ACM Transactions on Mathematical Software, Vol. 12, No. 3, September 1986, Pages 265-273 3. Amos, D. E., Remark on Algorithm 644, ACM Transactions on Mathematical Software, Vol. 16, No. 4, December 1990, Page 404 4. Amos, D. E., Remark on Algorithm 644 (Improvements in Algorithm 644), ACM Transactions on Mathematical Software, (Estimated date 1995) 5. Cody, W. J., Algorithm 665, MACHAR: A Subroutine to Dynamically Determine Machine Parameters, ACM Transactions on Mathematical Software, Vol. 14, No. 4, December 1988, Pages 303-311 C*** cbsubs.f ----------------------------------------------------------------- C>>> CBSUBS.FOR: Single precision subroutines ----------------------------------------------------------------- SUBROUTINE CBESH(Z, FNU, KODE, M, N, CY, NZ, IERR) C***BEGIN PROLOGUE CBESH C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C***CATEGORY NO. B5K C***KEYWORDS H-BESSEL FUNCTIONS,BESSEL FUNCTIONS OF COMPLEX ARGUMENT, C BESSEL FUNCTIONS OF THIRD KIND,HANKEL FUNCTIONS C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES C***PURPOSE TO COMPUTE THE H-BESSEL FUNCTIONS OF A COMPLEX ARGUMENT C***DESCRIPTION C C ON KODE=1, CBESH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX C HANKEL (BESSEL) FUNCTIONS CY(J)=H(M,FNU+J-1,Z) FOR KINDS M=1 C OR 2, REAL, NONNEGATIVE ORDERS FNU+J-1, J=1,...,N, AND COMPLEX C Z.NE.CMPLX(0.0E0,0.0E0) IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI. C ON KODE=2, CBESH COMPUTES THE SCALED HANKEL FUNCTIONS C C CY(I)=H(M,FNU+J-1,Z)*EXP(-MM*Z*I) MM=3-2M, I**2=-1. C C WHICH REMOVES THE EXPONENTIAL BEHAVIOR IN BOTH THE UPPER C AND LOWER HALF PLANES. DEFINITIONS AND NOTATION ARE FOUND IN C THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS (REF. 1). C C INPUT C Z - Z=CMPLX(X,Y), Z.NE.CMPLX(0.,0.),-PI.LT.ARG(Z).LE.PI C FNU - ORDER OF INITIAL H FUNCTION, FNU.GE.0.0E0 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION C KODE= 1 RETURNS C CY(J)=H(M,FNU+J-1,Z), J=1,...,N C = 2 RETURNS C CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M)) C J=1,...,N , I**2=-1 C M - KIND OF HANKEL FUNCTION, M=1 OR 2 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 C C OUTPUT C CY - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN C VALUES FOR THE SEQUENCE C CY(J)=H(M,FNU+J-1,Z) OR C CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M)) J=1,...,N C DEPENDING ON KODE, I**2=-1. C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW, C NZ= 0 , NORMAL RETURN C NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO C DUE TO UNDERFLOW, CY(J)=CMPLX(0.0,0.0) C J=1,...,NZ WHEN Y.GT.0.0 AND M=1 OR C Y.LT.0.0 AND M=2. FOR THE COMPLMENTARY C HALF PLANES, NZ STATES ONLY THE NUMBER C OF UNDERFLOWS. C IERR -ERROR FLAG C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED C IERR=1, INPUT ERROR - NO COMPUTATION C IERR=2, OVERFLOW - NO COMPUTATION, FNU+N-1 TOO C LARGE OR CABS(Z) TOO SMALL OR BOTH C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT C REDUCTION PRODUCE LESS THAN HALF OF MACHINE C ACCURACY C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- C CANCE BY ARGUMENT REDUCTION C IERR=5, ERROR - NO COMPUTATION, C ALGORITHM TERMINATION CONDITION NOT MET C C***LONG DESCRIPTION C C THE COMPUTATION IS CARRIED OUT BY THE RELATION C C H(M,FNU,Z)=(1/MP)*EXP(-MP*FNU)*K(FNU,Z*EXP(-MP)) C MP=MM*HPI*I, MM=3-2*M, HPI=PI/2, I**2=-1 C C FOR M=1 OR 2 WHERE THE K BESSEL FUNCTION IS COMPUTED FOR THE C RIGHT HALF PLANE RE(Z).GE.0.0. THE K FUNCTION IS CONTINUED C TO THE LEFT HALF PLANE BY THE RELATION C C K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z) C MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1 C C WHERE I(FNU,Z) IS THE I BESSEL FUNCTION. C C EXPONENTIAL DECAY OF H(M,FNU,Z) OCCURS IN THE UPPER HALF Z C PLANE FOR M=1 AND THE LOWER HALF Z PLANE FOR M=2. EXPONENTIAL C GROWTH OCCURS IN THE COMPLEMENTARY HALF PLANES. SCALING C BY EXP(-MM*Z*I) REMOVES THE EXPONENTIAL BEHAVIOR IN THE C WHOLE Z PLANE FOR Z TO INFINITY. C C FOR NEGATIVE ORDERS,THE FORMULAE C C H(1,-FNU,Z) = H(1,FNU,Z)*CEXP( PI*FNU*I) C H(2,-FNU,Z) = H(2,FNU,Z)*CEXP(-PI*FNU*I) C I**2=-1 C C CAN BE USED. C C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG C IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. C C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, C OR -PI/2+P. C C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF C COMMERCE, 1955. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C BY D. E. AMOS, SAND83-0083, MAY, 1983. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 C C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- C 1018, MAY, 1985 C C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, ACM C TRANS. MATH. SOFTWARE, VOL. 12, NO. 3, SEPTEMBER 1986, C PP 265-273. C C***ROUTINES CALLED CACON,CBKNU,CBUNK,CUOIK,I1MACH,R1MACH C***END PROLOGUE CBESH C COMPLEX CY, Z, ZN, ZT, CSGN REAL AA, ALIM, ALN, ARG, AZ, CPN, DIG, ELIM, FMM, FN, FNU, FNUL, * HPI, RHPI, RL, R1M5, SGN, SPN, TOL, UFL, XN, XX, YN, YY, R1MACH, * BB, ASCLE, RTOL, ATOL INTEGER I, IERR, INU, INUH, IR, K, KODE, K1, K2, M, * MM, MR, N, NN, NUF, NW, NZ, I1MACH DIMENSION CY(N) C DATA HPI /1.57079632679489662E0/ C C***FIRST EXECUTABLE STATEMENT CBESH NZ=0 XX = REAL(Z) YY = AIMAG(Z) IERR = 0 IF (XX.EQ.0.0E0 .AND. YY.EQ.0.0E0) IERR=1 IF (FNU.LT.0.0E0) IERR=1 IF (M.LT.1 .OR. M.GT.2) IERR=1 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 IF (N.LT.1) IERR=1 IF (IERR.NE.0) RETURN NN = N C----------------------------------------------------------------------- C SET PARAMETERS RELATED TO MACHINE CONSTANTS. C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU C----------------------------------------------------------------------- TOL = AMAX1(R1MACH(4),1.0E-18) K1 = I1MACH(12) K2 = I1MACH(13) R1M5 = R1MACH(5) K = MIN0(IABS(K1),IABS(K2)) ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) K1 = I1MACH(11) - 1 AA = R1M5*FLOAT(K1) DIG = AMIN1(AA,18.0E0) AA = AA*2.303E0 ALIM = ELIM + AMAX1(-AA,-41.45E0) FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) RL = 1.2E0*DIG + 3.0E0 FN = FNU + FLOAT(NN-1) MM = 3 - M - M FMM = FLOAT(MM) ZN = Z*CMPLX(0.0E0,-FMM) XN = REAL(ZN) YN = AIMAG(ZN) AZ = CABS(Z) C----------------------------------------------------------------------- C TEST FOR RANGE C----------------------------------------------------------------------- AA = 0.5E0/TOL BB=FLOAT(I1MACH(9))*0.5E0 AA=AMIN1(AA,BB) IF(AZ.GT.AA) GO TO 240 IF(FN.GT.AA) GO TO 240 AA=SQRT(AA) IF(AZ.GT.AA) IERR=3 IF(FN.GT.AA) IERR=3 C----------------------------------------------------------------------- C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE C----------------------------------------------------------------------- UFL = R1MACH(1)*1.0E+3 IF (AZ.LT.UFL) GO TO 220 IF (FNU.GT.FNUL) GO TO 90 IF (FN.LE.1.0E0) GO TO 70 IF (FN.GT.2.0E0) GO TO 60 IF (AZ.GT.TOL) GO TO 70 ARG = 0.5E0*AZ ALN = -FN*ALOG(ARG) IF (ALN.GT.ELIM) GO TO 220 GO TO 70 60 CONTINUE CALL CUOIK(ZN, FNU, KODE, 2, NN, CY, NUF, TOL, ELIM, ALIM) IF (NUF.LT.0) GO TO 220 NZ = NZ + NUF NN = NN - NUF C----------------------------------------------------------------------- C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I C----------------------------------------------------------------------- IF (NN.EQ.0) GO TO 130 70 CONTINUE IF ((XN.LT.0.0E0) .OR. (XN.EQ.0.0E0 .AND. YN.LT.0.0E0 .AND. * M.EQ.2)) GO TO 80 C----------------------------------------------------------------------- C RIGHT HALF PLANE COMPUTATION, XN.GE.0. .AND. (XN.NE.0. .OR. C YN.GE.0. .OR. M=1) C----------------------------------------------------------------------- CALL CBKNU(ZN, FNU, KODE, NN, CY, NZ, TOL, ELIM, ALIM) GO TO 110 C----------------------------------------------------------------------- C LEFT HALF PLANE COMPUTATION C----------------------------------------------------------------------- 80 CONTINUE MR = -MM CALL CACON(ZN, FNU, KODE, MR, NN, CY, NW, RL, FNUL, TOL, ELIM, * ALIM) IF (NW.LT.0) GO TO 230 NZ=NW GO TO 110 90 CONTINUE C----------------------------------------------------------------------- C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL C----------------------------------------------------------------------- MR = 0 IF ((XN.GE.0.0E0) .AND. (XN.NE.0.0E0 .OR. YN.GE.0.0E0 .OR. * M.NE.2)) GO TO 100 MR = -MM IF (XN.EQ.0.0E0 .AND. YN.LT.0.0E0) ZN = -ZN 100 CONTINUE CALL CBUNK(ZN, FNU, KODE, MR, NN, CY, NW, TOL, ELIM, ALIM) IF (NW.LT.0) GO TO 230 NZ = NZ + NW 110 CONTINUE C----------------------------------------------------------------------- C H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT) C C ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2 C----------------------------------------------------------------------- SGN = SIGN(HPI,-FMM) C----------------------------------------------------------------------- C CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE C WHEN FNU IS LARGE C----------------------------------------------------------------------- INU = INT(FNU) INUH = INU/2 IR = INU - 2*INUH ARG = (FNU-FLOAT(INU-IR))*SGN RHPI = 1.0E0/SGN CPN = RHPI*COS(ARG) SPN = RHPI*SIN(ARG) C ZN = CMPLX(-SPN,CPN) CSGN = CMPLX(-SPN,CPN) C IF (MOD(INUH,2).EQ.1) ZN = -ZN IF (MOD(INUH,2).EQ.1) CSGN = -CSGN ZT = CMPLX(0.0E0,-FMM) RTOL = 1.0E0/TOL ASCLE = UFL*RTOL DO 120 I=1,NN C CY(I) = CY(I)*ZN C ZN = ZN*ZT ZN=CY(I) AA=REAL(ZN) BB=AIMAG(ZN) ATOL=1.0E0 IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 125 ZN = ZN*CMPLX(RTOL,0.0E0) ATOL = TOL 125 CONTINUE ZN = ZN*CSGN CY(I) = ZN*CMPLX(ATOL,0.0E0) CSGN = CSGN*ZT 120 CONTINUE RETURN 130 CONTINUE IF (XN.LT.0.0E0) GO TO 220 RETURN 220 CONTINUE IERR=2 NZ=0 RETURN 230 CONTINUE IF(NW.EQ.(-1)) GO TO 220 NZ=0 IERR=5 RETURN 240 CONTINUE NZ=0 IERR=4 RETURN END SUBROUTINE CBESI(Z, FNU, KODE, N, CY, NZ, IERR) C***BEGIN PROLOGUE CBESI C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C***CATEGORY NO. B5K C***KEYWORDS I-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION, C MODIFIED BESSEL FUNCTION OF THE FIRST KIND C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES C***PURPOSE TO COMPUTE I-BESSEL FUNCTIONS OF COMPLEX ARGUMENT C***DESCRIPTION C C ON KODE=1, CBESI COMPUTES AN N MEMBER SEQUENCE OF COMPLEX C BESSEL FUNCTIONS CY(J)=I(FNU+J-1,Z) FOR REAL, NONNEGATIVE C ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z IN THE CUT PLANE C -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESI RETURNS THE SCALED C FUNCTIONS C C CY(J)=EXP(-ABS(X))*I(FNU+J-1,Z) J = 1,...,N , X=REAL(Z) C C WITH THE EXPONENTIAL GROWTH REMOVED IN BOTH THE LEFT AND C RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND C NOTATION ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL C FUNCTIONS (REF.1) C C INPUT C Z - Z=CMPLX(X,Y), -PI.LT.ARG(Z).LE.PI C FNU - ORDER OF INITIAL I FUNCTION, FNU.GE.0.0E0 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION C KODE= 1 RETURNS C CY(J)=I(FNU+J-1,Z), J=1,...,N C = 2 RETURNS C CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X)), J=1,...,N C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 C C OUTPUT C CY - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN C VALUES FOR THE SEQUENCE C CY(J)=I(FNU+J-1,Z) OR C CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X)) J=1,...,N C DEPENDING ON KODE, X=REAL(Z) C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW, C NZ= 0 , NORMAL RETURN C NZ.GT.0 , LAST NZ COMPONENTS OF CY SET TO ZERO C DUE TO UNDERFLOW, CY(J)=CMPLX(0.0,0.0), C J = N-NZ+1,...,N C IERR - ERROR FLAG C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED C IERR=1, INPUT ERROR - NO COMPUTATION C IERR=2, OVERFLOW - NO COMPUTATION, REAL(Z) TOO C LARGE ON KODE=1 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT C REDUCTION PRODUCE LESS THAN HALF OF MACHINE C ACCURACY C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- C CANCE BY ARGUMENT REDUCTION C IERR=5, ERROR - NO COMPUTATION, C ALGORITHM TERMINATION CONDITION NOT MET C C***LONG DESCRIPTION C C THE COMPUTATION IS CARRIED OUT BY THE POWER SERIES FOR C SMALL CABS(Z), THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z), C THE MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN AND A C NEUMANN SERIES FOR IMTERMEDIATE MAGNITUDES, AND THE C UNIFORM ASYMPTOTIC EXPANSIONS FOR I(FNU,Z) AND J(FNU,Z) C FOR LARGE ORDERS. BACKWARD RECURRENCE IS USED TO GENERATE C SEQUENCES OR REDUCE ORDERS WHEN NECESSARY. C C THE CALCULATIONS ABOVE ARE DONE IN THE RIGHT HALF PLANE AND C CONTINUED INTO THE LEFT HALF PLANE BY THE FORMULA C C I(FNU,Z*EXP(M*PI)) = EXP(M*PI*FNU)*I(FNU,Z) REAL(Z).GT.0.0 C M = +I OR -I, I**2=-1 C C FOR NEGATIVE ORDERS,THE FORMULA C C I(-FNU,Z) = I(FNU,Z) + (2/PI)*SIN(PI*FNU)*K(FNU,Z) C C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO INTEGERS, THE C THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE POSITIVE C INTEGER,THE MAGNITUDE OF I(-FNU,Z)=I(FNU,Z) IS A LARGE C NEGATIVE POWER OF TEN. BUT WHEN FNU IS NOT AN INTEGER, C K(FNU,Z) DOMINATES IN MAGNITUDE WITH A LARGE POSITIVE POWER OF C TEN AND THE MOST THAT THE SECOND TERM CAN BE REDUCED IS BY C UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, WIDE CHANGES CAN C OCCUR WITHIN UNIT ROUNDOFF OF A LARGE INTEGER FOR FNU. HERE, C LARGE MEANS FNU.GT.CABS(Z). C C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG C IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. C C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, C OR -PI/2+P. C C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF C COMMERCE, 1955. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C BY D. E. AMOS, SAND83-0083, MAY, 1983. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 C C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- C 1018, MAY, 1985 C C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, ACM C TRANS. MATH. SOFTWARE, VOL. 12, NO. 3, SEPTEMBER 1986, C PP 265-273. C C***ROUTINES CALLED CBINU,I1MACH,R1MACH C***END PROLOGUE CBESI COMPLEX CONE, CSGN, CY, Z, ZN REAL AA, ALIM, ARG, DIG, ELIM, FNU, FNUL, PI, RL, R1M5, S1, S2, * TOL, XX, YY, R1MACH, AZ, FN, BB, ASCLE, RTOL, ATOL INTEGER I, IERR, INU, K, KODE, K1, K2, N, NN, NZ, I1MACH DIMENSION CY(N) DATA PI /3.14159265358979324E0/ DATA CONE / (1.0E0,0.0E0) / C C***FIRST EXECUTABLE STATEMENT CBESI IERR = 0 NZ=0 IF (FNU.LT.0.0E0) IERR=1 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 IF (N.LT.1) IERR=1 IF (IERR.NE.0) RETURN XX = REAL(Z) YY = AIMAG(Z) C----------------------------------------------------------------------- C SET PARAMETERS RELATED TO MACHINE CONSTANTS. C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU. C----------------------------------------------------------------------- TOL = AMAX1(R1MACH(4),1.0E-18) K1 = I1MACH(12) K2 = I1MACH(13) R1M5 = R1MACH(5) K = MIN0(IABS(K1),IABS(K2)) ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) K1 = I1MACH(11) - 1 AA = R1M5*FLOAT(K1) DIG = AMIN1(AA,18.0E0) AA = AA*2.303E0 ALIM = ELIM + AMAX1(-AA,-41.45E0) RL = 1.2E0*DIG + 3.0E0 FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) AZ = CABS(Z) C----------------------------------------------------------------------- C TEST FOR RANGE C----------------------------------------------------------------------- AA = 0.5E0/TOL BB=FLOAT(I1MACH(9))*0.5E0 AA=AMIN1(AA,BB) IF(AZ.GT.AA) GO TO 140 FN=FNU+FLOAT(N-1) IF(FN.GT.AA) GO TO 140 AA=SQRT(AA) IF(AZ.GT.AA) IERR=3 IF(FN.GT.AA) IERR=3 ZN = Z CSGN = CONE IF (XX.GE.0.0E0) GO TO 40 ZN = -Z C----------------------------------------------------------------------- C CALCULATE CSGN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE C WHEN FNU IS LARGE C----------------------------------------------------------------------- INU = INT(FNU) ARG = (FNU-FLOAT(INU))*PI IF (YY.LT.0.0E0) ARG = -ARG S1 = COS(ARG) S2 = SIN(ARG) CSGN = CMPLX(S1,S2) IF (MOD(INU,2).EQ.1) CSGN = -CSGN 40 CONTINUE CALL CBINU(ZN, FNU, KODE, N, CY, NZ, RL, FNUL, TOL, ELIM, ALIM) IF (NZ.LT.0) GO TO 120 IF (XX.GE.0.0E0) RETURN C----------------------------------------------------------------------- C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE C----------------------------------------------------------------------- NN = N - NZ IF (NN.EQ.0) RETURN RTOL = 1.0E0/TOL ASCLE = R1MACH(1)*RTOL*1.0E+3 DO 50 I=1,NN C CY(I) = CY(I)*CSGN ZN=CY(I) AA=REAL(ZN) BB=AIMAG(ZN) ATOL=1.0E0 IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 55 ZN = ZN*CMPLX(RTOL,0.0E0) ATOL = TOL 55 CONTINUE ZN = ZN*CSGN CY(I) = ZN*CMPLX(ATOL,0.0E0) CSGN = -CSGN 50 CONTINUE RETURN 120 CONTINUE IF(NZ.EQ.(-2)) GO TO 130 NZ = 0 IERR=2 RETURN 130 CONTINUE NZ=0 IERR=5 RETURN 140 CONTINUE NZ=0 IERR=4 RETURN END SUBROUTINE CBESJ(Z, FNU, KODE, N, CY, NZ, IERR) C***BEGIN PROLOGUE CBESJ C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C***CATEGORY NO. B5K C***KEYWORDS J-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT, C BESSEL FUNCTION OF FIRST KIND C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES C***PURPOSE TO COMPUTE THE J-BESSEL FUNCTION OF A COMPLEX ARGUMENT C***DESCRIPTION C C ON KODE=1, CBESJ COMPUTES AN N MEMBER SEQUENCE OF COMPLEX C BESSEL FUNCTIONS CY(I)=J(FNU+I-1,Z) FOR REAL, NONNEGATIVE C ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE C -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESJ RETURNS THE SCALED C FUNCTIONS C C CY(I)=EXP(-ABS(Y))*J(FNU+I-1,Z) I = 1,...,N , Y=AIMAG(Z) C C WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND C LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS C (REF. 1). C C INPUT C Z - Z=CMPLX(X,Y), -PI.LT.ARG(Z).LE.PI C FNU - ORDER OF INITIAL J FUNCTION, FNU.GE.0.0E0 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION C KODE= 1 RETURNS C CY(I)=J(FNU+I-1,Z), I=1,...,N C = 2 RETURNS C CY(I)=J(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,... C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 C C OUTPUT C CY - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN C VALUES FOR THE SEQUENCE C CY(I)=J(FNU+I-1,Z) OR C CY(I)=J(FNU+I-1,Z)*EXP(-ABS(Y)) I=1,...,N C DEPENDING ON KODE, Y=AIMAG(Z). C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW, C NZ= 0 , NORMAL RETURN C NZ.GT.0 , LAST NZ COMPONENTS OF CY SET TO ZERO C DUE TO UNDERFLOW, CY(I)=CMPLX(0.0,0.0), C I = N-NZ+1,...,N C IERR - ERROR FLAG C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED C IERR=1, INPUT ERROR - NO COMPUTATION C IERR=2, OVERFLOW - NO COMPUTATION, AIMAG(Z) C TOO LARGE ON KODE=1 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT C REDUCTION PRODUCE LESS THAN HALF OF MACHINE C ACCURACY C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- C CANCE BY ARGUMENT REDUCTION C IERR=5, ERROR - NO COMPUTATION, C ALGORITHM TERMINATION CONDITION NOT MET C C***LONG DESCRIPTION C C THE COMPUTATION IS CARRIED OUT BY THE FORMULA C C J(FNU,Z)=EXP( FNU*PI*I/2)*I(FNU,-I*Z) AIMAG(Z).GE.0.0 C C J(FNU,Z)=EXP(-FNU*PI*I/2)*I(FNU, I*Z) AIMAG(Z).LT.0.0 C C WHERE I**2 = -1 AND I(FNU,Z) IS THE I BESSEL FUNCTION. C C FOR NEGATIVE ORDERS,THE FORMULA C C J(-FNU,Z) = J(FNU,Z)*COS(PI*FNU) - Y(FNU,Z)*SIN(PI*FNU) C C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO INTEGERS, THE C THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE POSITIVE C INTEGER,THE MAGNITUDE OF J(-FNU,Z)=J(FNU,Z)*COS(PI*FNU) IS A C LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS NOT AN INTEGER, C Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A LARGE POSITIVE POWER OF C TEN AND THE MOST THAT THE SECOND TERM CAN BE REDUCED IS BY C UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, WIDE CHANGES CAN C OCCUR WITHIN UNIT ROUNDOFF OF A LARGE INTEGER FOR FNU. HERE, C LARGE MEANS FNU.GT.CABS(Z). C C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG C IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. C C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, C OR -PI/2+P. C C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF C COMMERCE, 1955. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C BY D. E. AMOS, SAND83-0083, MAY, 1983. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 C C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- C 1018, MAY, 1985 C C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, ACM C TRANS. MATH. SOFTWARE, VOL. 12, NO. 3, SEPTEMBER 1986, C PP 265-273. C C***ROUTINES CALLED CBINU,I1MACH,R1MACH C***END PROLOGUE CBESJ C COMPLEX CI, CSGN, CY, Z, ZN REAL AA, ALIM, ARG, DIG, ELIM, FNU, FNUL, HPI, RL, R1, R1M5, R2, * TOL, YY, R1MACH, AZ, FN, BB, ASCLE, RTOL, ATOL INTEGER I, IERR, INU, INUH, IR, KODE, K1, K2, N, NL, NZ, I1MACH, K DIMENSION CY(N) DATA HPI /1.57079632679489662E0/ C C***FIRST EXECUTABLE STATEMENT CBESJ IERR = 0 NZ=0 IF (FNU.LT.0.0E0) IERR=1 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 IF (N.LT.1) IERR=1 IF (IERR.NE.0) RETURN C----------------------------------------------------------------------- C SET PARAMETERS RELATED TO MACHINE CONSTANTS. C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU. C----------------------------------------------------------------------- TOL = AMAX1(R1MACH(4),1.0E-18) K1 = I1MACH(12) K2 = I1MACH(13) R1M5 = R1MACH(5) K = MIN0(IABS(K1),IABS(K2)) ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) K1 = I1MACH(11) - 1 AA = R1M5*FLOAT(K1) DIG = AMIN1(AA,18.0E0) AA = AA*2.303E0 ALIM = ELIM + AMAX1(-AA,-41.45E0) RL = 1.2E0*DIG + 3.0E0 FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) CI = CMPLX(0.0E0,1.0E0) YY = AIMAG(Z) AZ = CABS(Z) C----------------------------------------------------------------------- C TEST FOR RANGE C----------------------------------------------------------------------- AA = 0.5E0/TOL BB=FLOAT(I1MACH(9))*0.5E0 AA=AMIN1(AA,BB) FN=FNU+FLOAT(N-1) IF(AZ.GT.AA) GO TO 140 IF(FN.GT.AA) GO TO 140 AA=SQRT(AA) IF(AZ.GT.AA) IERR=3 IF(FN.GT.AA) IERR=3 C----------------------------------------------------------------------- C CALCULATE CSGN=EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE C WHEN FNU IS LARGE C----------------------------------------------------------------------- INU = INT(FNU) INUH = INU/2 IR = INU - 2*INUH ARG = (FNU-FLOAT(INU-IR))*HPI R1 = COS(ARG) R2 = SIN(ARG) CSGN = CMPLX(R1,R2) IF (MOD(INUH,2).EQ.1) CSGN = -CSGN C----------------------------------------------------------------------- C ZN IS IN THE RIGHT HALF PLANE C----------------------------------------------------------------------- ZN = -Z*CI IF (YY.GE.0.0E0) GO TO 40 ZN = -ZN CSGN = CONJG(CSGN) CI = CONJG(CI) 40 CONTINUE CALL CBINU(ZN, FNU, KODE, N, CY, NZ, RL, FNUL, TOL, ELIM, ALIM) IF (NZ.LT.0) GO TO 120 NL = N - NZ IF (NL.EQ.0) RETURN RTOL = 1.0E0/TOL ASCLE = R1MACH(1)*RTOL*1.0E+3 DO 50 I=1,NL C CY(I)=CY(I)*CSGN ZN=CY(I) AA=REAL(ZN) BB=AIMAG(ZN) ATOL=1.0E0 IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 55 ZN = ZN*CMPLX(RTOL,0.0E0) ATOL = TOL 55 CONTINUE ZN = ZN*CSGN CY(I) = ZN*CMPLX(ATOL,0.0E0) CSGN = CSGN*CI 50 CONTINUE RETURN 120 CONTINUE IF(NZ.EQ.(-2)) GO TO 130 NZ = 0 IERR = 2 RETURN 130 CONTINUE NZ=0 IERR=5 RETURN 140 CONTINUE NZ=0 IERR=4 RETURN END SUBROUTINE CBESK(Z, FNU, KODE, N, CY, NZ, IERR) C***BEGIN PROLOGUE CBESK C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C***CATEGORY NO. B5K C***KEYWORDS K-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION, C MODIFIED BESSEL FUNCTION OF THE SECOND KIND, C BESSEL FUNCTION OF THE THIRD KIND C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES C***PURPOSE TO COMPUTE K-BESSEL FUNCTIONS OF COMPLEX ARGUMENT C***DESCRIPTION C C ON KODE=1, CBESK COMPUTES AN N MEMBER SEQUENCE OF COMPLEX C BESSEL FUNCTIONS CY(J)=K(FNU+J-1,Z) FOR REAL, NONNEGATIVE C ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z.NE.CMPLX(0.0,0.0) C IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESK C RETURNS THE SCALED K FUNCTIONS, C C CY(J)=EXP(Z)*K(FNU+J-1,Z) , J=1,...,N, C C WHICH REMOVE THE EXPONENTIAL BEHAVIOR IN BOTH THE LEFT AND C RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND C NOTATION ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL C FUNCTIONS (REF. 1). C C INPUT C Z - Z=CMPLX(X,Y),Z.NE.CMPLX(0.,0.),-PI.LT.ARG(Z).LE.PI C FNU - ORDER OF INITIAL K FUNCTION, FNU.GE.0.0E0 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION C KODE= 1 RETURNS C CY(I)=K(FNU+I-1,Z), I=1,...,N C = 2 RETURNS C CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N C C OUTPUT C CY - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN C VALUES FOR THE SEQUENCE C CY(I)=K(FNU+I-1,Z), I=1,...,N OR C CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N C DEPENDING ON KODE C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW. C NZ= 0 , NORMAL RETURN C NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO C DUE TO UNDERFLOW, CY(I)=CMPLX(0.0,0.0), C I=1,...,N WHEN X.GE.0.0. WHEN X.LT.0.0 C NZ STATES ONLY THE NUMBER OF UNDERFLOWS C IN THE SEQUENCE. C IERR - ERROR FLAG C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED C IERR=1, INPUT ERROR - NO COMPUTATION C IERR=2, OVERFLOW - NO COMPUTATION, FNU+N-1 IS C TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT C REDUCTION PRODUCE LESS THAN HALF OF MACHINE C ACCURACY C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- C CANCE BY ARGUMENT REDUCTION C IERR=5, ERROR - NO COMPUTATION, C ALGORITHM TERMINATION CONDITION NOT MET C C***LONG DESCRIPTION C C EQUATIONS OF THE REFERENCE ARE IMPLEMENTED FOR SMALL ORDERS C DNU AND DNU+1.0 IN THE RIGHT HALF PLANE X.GE.0.0. FORWARD C RECURRENCE GENERATES HIGHER ORDERS. K IS CONTINUED TO THE LEFT C HALF PLANE BY THE RELATION C C K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z) C MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1 C C WHERE I(FNU,Z) IS THE I BESSEL FUNCTION. C C FOR LARGE ORDERS, FNU.GT.FNUL, THE K FUNCTION IS COMPUTED C BY MEANS OF ITS UNIFORM ASYMPTOTIC EXPANSIONS. C C FOR NEGATIVE ORDERS, THE FORMULA C C K(-FNU,Z) = K(FNU,Z) C C CAN BE USED. C C CBESK ASSUMES THAT A SIGNIFICANT DIGIT SINH(X) FUNCTION IS C AVAILABLE. C C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG C IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. C C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, C OR -PI/2+P. C C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF C COMMERCE, 1955. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C BY D. E. AMOS, SAND83-0083, MAY, 1983. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983. C C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- C 1018, MAY, 1985 C C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, ACM C TRANS. MATH. SOFTWARE, VOL. 12, NO. 3, SEPTEMBER 1986, C PP 265-273. C C***ROUTINES CALLED CACON,CBKNU,CBUNK,CUOIK,I1MACH,R1MACH C***END PROLOGUE CBESK C COMPLEX CY, Z REAL AA, ALIM, ALN, ARG, AZ, DIG, ELIM, FN, FNU, FNUL, RL, R1M5, * TOL, UFL, XX, YY, R1MACH, BB INTEGER IERR, K, KODE, K1, K2, MR, N, NN, NUF, NW, NZ, I1MACH DIMENSION CY(N) C***FIRST EXECUTABLE STATEMENT CBESK IERR = 0 NZ=0 XX = REAL(Z) YY = AIMAG(Z) IF (YY.EQ.0.0E0 .AND. XX.EQ.0.0E0) IERR=1 IF (FNU.LT.0.0E0) IERR=1 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 IF (N.LT.1) IERR=1 IF (IERR.NE.0) RETURN NN = N C----------------------------------------------------------------------- C SET PARAMETERS RELATED TO MACHINE CONSTANTS. C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU C----------------------------------------------------------------------- TOL = AMAX1(R1MACH(4),1.0E-18) K1 = I1MACH(12) K2 = I1MACH(13) R1M5 = R1MACH(5) K = MIN0(IABS(K1),IABS(K2)) ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) K1 = I1MACH(11) - 1 AA = R1M5*FLOAT(K1) DIG = AMIN1(AA,18.0E0) AA = AA*2.303E0 ALIM = ELIM + AMAX1(-AA,-41.45E0) FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) RL = 1.2E0*DIG + 3.0E0 AZ = CABS(Z) FN = FNU + FLOAT(NN-1) C----------------------------------------------------------------------- C TEST FOR RANGE C----------------------------------------------------------------------- AA = 0.5E0/TOL BB=FLOAT(I1MACH(9))*0.5E0 AA=AMIN1(AA,BB) IF(AZ.GT.AA) GO TO 210 IF(FN.GT.AA) GO TO 210 AA=SQRT(AA) IF(AZ.GT.AA) IERR=3 IF(FN.GT.AA) IERR=3 C----------------------------------------------------------------------- C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE C----------------------------------------------------------------------- C UFL = EXP(-ELIM) UFL = R1MACH(1)*1.0E+3 IF (AZ.LT.UFL) GO TO 180 IF (FNU.GT.FNUL) GO TO 80 IF (FN.LE.1.0E0) GO TO 60 IF (FN.GT.2.0E0) GO TO 50 IF (AZ.GT.TOL) GO TO 60 ARG = 0.5E0*AZ ALN = -FN*ALOG(ARG) IF (ALN.GT.ELIM) GO TO 180 GO TO 60 50 CONTINUE CALL CUOIK(Z, FNU, KODE, 2, NN, CY, NUF, TOL, ELIM, ALIM) IF (NUF.LT.0) GO TO 180 NZ = NZ + NUF NN = NN - NUF C----------------------------------------------------------------------- C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I C----------------------------------------------------------------------- IF (NN.EQ.0) GO TO 100 60 CONTINUE IF (XX.LT.0.0E0) GO TO 70 C----------------------------------------------------------------------- C RIGHT HALF PLANE COMPUTATION, REAL(Z).GE.0. C----------------------------------------------------------------------- CALL CBKNU(Z, FNU, KODE, NN, CY, NW, TOL, ELIM, ALIM) IF (NW.LT.0) GO TO 200 NZ=NW RETURN C----------------------------------------------------------------------- C LEFT HALF PLANE COMPUTATION C PI/2.LT.ARG(Z).LE.PI AND -PI.LT.ARG(Z).LT.-PI/2. C----------------------------------------------------------------------- 70 CONTINUE IF (NZ.NE.0) GO TO 180 MR = 1 IF (YY.LT.0.0E0) MR = -1 CALL CACON(Z, FNU, KODE, MR, NN, CY, NW, RL, FNUL, TOL, ELIM, * ALIM) IF (NW.LT.0) GO TO 200 NZ=NW RETURN C----------------------------------------------------------------------- C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL C----------------------------------------------------------------------- 80 CONTINUE MR = 0 IF (XX.GE.0.0E0) GO TO 90 MR = 1 IF (YY.LT.0.0E0) MR = -1 90 CONTINUE CALL CBUNK(Z, FNU, KODE, MR, NN, CY, NW, TOL, ELIM, ALIM) IF (NW.LT.0) GO TO 200 NZ = NZ + NW RETURN 100 CONTINUE IF (XX.LT.0.0E0) GO TO 180 RETURN 180 CONTINUE NZ = 0 IERR=2 RETURN 200 CONTINUE IF(NW.EQ.(-1)) GO TO 180 NZ=0 IERR=5 RETURN 210 CONTINUE NZ=0 IERR=4 RETURN END SUBROUTINE CBESY(Z, FNU, KODE, N, CY, NZ, CWRK, IERR) C***BEGIN PROLOGUE CBESY C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C***CATEGORY NO. B5K C***KEYWORDS Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT, C BESSEL FUNCTION OF SECOND KIND C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES C***PURPOSE TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT C***DESCRIPTION C C ON KODE=1, CBESY COMPUTES AN N MEMBER SEQUENCE OF COMPLEX C BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE C ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE C -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESY RETURNS THE SCALED C FUNCTIONS C C CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z) I = 1,...,N , Y=AIMAG(Z) C C WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND C LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS C (REF. 1). C C INPUT C Z - Z=CMPLX(X,Y), Z.NE.CMPLX(0.,0.),-PI.LT.ARG(Z).LE.PI C FNU - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0E0 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION C KODE= 1 RETURNS C CY(I)=Y(FNU+I-1,Z), I=1,...,N C = 2 RETURNS C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N C WHERE Y=AIMAG(Z) C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 C CWRK - A COMPLEX WORK VECTOR OF DIMENSION AT LEAST N C C OUTPUT C CY - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN C VALUES FOR THE SEQUENCE C CY(I)=Y(FNU+I-1,Z) OR C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)) I=1,...,N C DEPENDING ON KODE. C NZ - NZ=0 , A NORMAL RETURN C NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO C UNDERFLOW (GENERALLY ON KODE=2) C IERR - ERROR FLAG C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED C IERR=1, INPUT ERROR - NO COMPUTATION C IERR=2, OVERFLOW - NO COMPUTATION, FNU+N-1 IS C TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT C REDUCTION PRODUCE LESS THAN HALF OF MACHINE C ACCURACY C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- C CANCE BY ARGUMENT REDUCTION C IERR=5, ERROR - NO COMPUTATION, C ALGORITHM TERMINATION CONDITION NOT MET C C***LONG DESCRIPTION C C THE COMPUTATION IS CARRIED OUT IN TERMS OF THE I(FNU,Z) AND C K(FNU,Z) BESSEL FUNCTIONS IN THE RIGHT HALF PLANE BY C C Y(FNU,Z) = I*CC*I(FNU,ARG) - (2/PI)*CONJG(CC)*K(FNU,ARG) C C Y(FNU,Z) = CONJG(Y(FNU,CONJG(Z))) C C FOR AIMAG(Z).GE.0 AND AIMAG(Z).LT.0 RESPECTIVELY, WHERE C CC=EXP(I*PI*FNU/2), ARG=Z*EXP(-I*PI/2) AND I**2=-1. C C FOR NEGATIVE ORDERS,THE FORMULA C C Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU) C C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD C INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE C POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)* C SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS C NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A C LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM C CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, C WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF C ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z). C C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG C IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. C C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, C OR -PI/2+P. C C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF C COMMERCE, 1955. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C BY D. E. AMOS, SAND83-0083, MAY, 1983. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 C C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- C 1018, MAY, 1985 C C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, ACM C TRANS. MATH. SOFTWARE, VOL. 12, NO. 3, SEPTEMBER 1986, C PP 265-273. C C***ROUTINES CALLED CBESI,CBESK,I1MACH,R1MACH C***END PROLOGUE CBESY C COMPLEX CWRK, CY, CI, CIP, CSGN, CSPN, EX, Z, ZU, ZV, ZZ, ZN REAL ARG, ELIM, EY, FNU, R1, R2, TAY, XX, YY, R1MACH, ASCLE, RTOL, * ATOL, TOL, AA, BB, FFNU, HPI, RHPI, R1M5 INTEGER I, IERR, IFNU, K, KODE, K1, K2, N, NZ, NZ1, NZ2, I1MACH, *I4 DIMENSION CY(N), CWRK(N), CIP(4) DATA CIP(1),CIP(2),CIP(3),CIP(4)/ * (1.0E0,0.0E0) , (0.0E0,1.0E0) , (-1.0E0,0.0E0) , (0.0E0,-1.0E0) / DATA HPI / 1.57079632679489662E0 / C***FIRST EXECUTABLE STATEMENT CBESY XX = REAL(Z) YY = AIMAG(Z) IERR = 0 NZ=0 IF (XX.EQ.0.0E0 .AND. YY.EQ.0.0E0) IERR=1 IF (FNU.LT.0.0E0) IERR=1 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 IF (N.LT.1) IERR=1 IF (IERR.NE.0) RETURN CI = CMPLX(0.0E0,1.0E0) ZZ=Z IF (YY.LT.0.0E0) ZZ=CONJG(Z) ZN = -CI*ZZ CALL CBESI(ZN, FNU, KODE, N, CY, NZ1, IERR) IF (IERR.NE.0.AND.IERR.NE.3) GO TO 90 CALL CBESK(ZN, FNU, KODE, N, CWRK, NZ2, IERR) IF (IERR.NE.0.AND.IERR.NE.3) GO TO 90 NZ = MIN0(NZ1,NZ2) IFNU = INT(FNU) FFNU = FNU - FLOAT(IFNU) ARG = HPI*FFNU CSGN = CMPLX(COS(ARG),SIN(ARG)) I4 = MOD(IFNU,4) + 1 CSGN = CSGN*CIP(I4) RHPI = 1.0E0/HPI CSPN = CONJG(CSGN)*CMPLX(RHPI,0.0E0) CSGN = CSGN*CI IF (KODE.EQ.2) GO TO 60 DO 50 I=1,N CY(I) = CSGN*CY(I)-CSPN*CWRK(I) CSGN = CI*CSGN CSPN = -CI*CSPN 50 CONTINUE IF (YY.LT.0.0E0) THEN DO 55 I=1,N CY(I)=CONJG(CY(I)) 55 CONTINUE ENDIF RETURN 60 CONTINUE R1 = COS(XX) R2 = SIN(XX) EX = CMPLX(R1,R2) TOL = AMAX1(R1MACH(4),1.0E-18) K1 = I1MACH(12) K2 = I1MACH(13) K = MIN0(IABS(K1),IABS(K2)) R1M5 = R1MACH(5) C----------------------------------------------------------------------- C ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT C----------------------------------------------------------------------- ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) EY = 0.0E0 TAY = ABS(YY+YY) IF (TAY.LT.ELIM) EY = EXP(-TAY) CSPN = EX*CMPLX(EY,0.0E0)*CSPN NZ = 0 RTOL = 1.0E0/TOL ASCLE = R1MACH(1)*RTOL*1.0E+3 DO 80 I=1,N C---------------------------------------------------------------------- C CY(I) = CSGN*CY(I)-CSPN*CWRK(I): PRODUCTS ARE COMPUTED IN C SCALED MODE IF CY(I) OR CWRK(I) ARE CLOSE TO UNDERFLOW TO C PREVENT UNDERFLOW IN AN INTERMEDIATE COMPUTATION. C---------------------------------------------------------------------- ZV = CWRK(I) AA=REAL(ZV) BB=AIMAG(ZV) ATOL=1.0E0 IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 75 ZV = ZV*CMPLX(RTOL,0.0E0) ATOL = TOL 75 CONTINUE ZV = ZV*CSPN ZV = ZV*CMPLX(ATOL,0.0E0) ZU = CY(I) AA=REAL(ZU) BB=AIMAG(ZU) ATOL=1.0E0 IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 85 ZU = ZU*CMPLX(RTOL,0.0E0) ATOL = TOL 85 CONTINUE ZU = ZU*CSGN ZU = ZU*CMPLX(ATOL,0.0E0) CY(I) = ZU - ZV IF (YY.LT.0.0E0) CY(I)=CONJG(CY(I)) IF (CY(I).EQ.CMPLX(0.0E0,0.0E0) .AND. EY.EQ.0.0E0) NZ = NZ + 1 CSGN = CI*CSGN CSPN = -CI*CSPN 80 CONTINUE RETURN 90 CONTINUE NZ = 0 RETURN END SUBROUTINE CAIRY(Z, ID, KODE, AI, NZ, IERR) C***BEGIN PROLOGUE CAIRY C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C***CATEGORY NO. B5K C***KEYWORDS AIRY FUNCTION,BESSEL FUNCTIONS OF ORDER ONE THIRD C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES C***PURPOSE TO COMPUTE AIRY FUNCTIONS AI(Z) AND DAI(Z) FOR COMPLEX Z C***DESCRIPTION C C ON KODE=1, CAIRY COMPUTES THE COMPLEX AIRY FUNCTION AI(Z) OR C ITS DERIVATIVE DAI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON C KODE=2, A SCALING OPTION CEXP(ZTA)*AI(Z) OR CEXP(ZTA)* C DAI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL DECAY IN C -PI/3.LT.ARG(Z).LT.PI/3 AND THE EXPONENTIAL GROWTH IN C PI/3.LT.ABS(ARG(Z)).LT.PI WHERE ZTA=(2/3)*Z*CSQRT(Z) C C WHILE THE AIRY FUNCTIONS AI(Z) AND DAI(Z)/DZ ARE ANALYTIC IN C THE WHOLE Z PLANE, THE CORRESPONDING SCALED FUNCTIONS DEFINED C FOR KODE=2 HAVE A CUT ALONG THE NEGATIVE REAL AXIS. C DEFINITIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF C MATHEMATICAL FUNCTIONS (REF. 1). C C INPUT C Z - Z=CMPLX(X,Y) C ID - ORDER OF DERIVATIVE, ID=0 OR ID=1 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION C KODE= 1 RETURNS C AI=AI(Z) ON ID=0 OR C AI=DAI(Z)/DZ ON ID=1 C = 2 RETURNS C AI=CEXP(ZTA)*AI(Z) ON ID=0 OR C AI=CEXP(ZTA)*DAI(Z)/DZ ON ID=1 WHERE C ZTA=(2/3)*Z*CSQRT(Z) C C OUTPUT C AI - COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND C KODE C NZ - UNDERFLOW INDICATOR C NZ= 0 , NORMAL RETURN C NZ= 1 , AI=CMPLX(0.0,0.0) DUE TO UNDERFLOW IN C -PI/3.LT.ARG(Z).LT.PI/3 ON KODE=1 C IERR - ERROR FLAG C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED C IERR=1, INPUT ERROR - NO COMPUTATION C IERR=2, OVERFLOW - NO COMPUTATION, REAL(ZTA) C TOO LARGE WITH KODE=1. C IERR=3, CABS(Z) LARGE - COMPUTATION COMPLETED C LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION C PRODUCE LESS THAN HALF OF MACHINE ACCURACY C IERR=4, CABS(Z) TOO LARGE - NO COMPUTATION C COMPLETE LOSS OF ACCURACY BY ARGUMENT C REDUCTION C IERR=5, ERROR - NO COMPUTATION, C ALGORITHM TERMINATION CONDITION NOT MET C C C***LONG DESCRIPTION C C AI AND DAI ARE COMPUTED FOR CABS(Z).GT.1.0 FROM THE K BESSEL C FUNCTIONS BY C C AI(Z)=C*SQRT(Z)*K(1/3,ZTA) , DAI(Z)=-C*Z*K(2/3,ZTA) C C=1.0/(PI*SQRT(3.0)) C ZTA=(2/3)*Z**(3/2) C C WITH THE POWER SERIES FOR CABS(Z).LE.1.0. C C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z IS LARGE, LOSSES C OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. CONSEQUENTLY, IF C THE MAGNITUDE OF ZETA=(2/3)*Z**1.5 EXCEEDS U1=SQRT(0.5/UR), C THEN LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR C FLAG IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. C ALSO, IF THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN C ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT C FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE C LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA C MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, C AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE C PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE C PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT- C ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG- C NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN C DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN C EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, C NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE C PRECISION ARITHMETIC. SIMILAR CONSIDERATIONS HOLD FOR OTHER C MACHINES. C C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, C OR -PI/2+P. C C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF C COMMERCE, 1955. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 C C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- C 1018, MAY, 1985 C C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, ACM C TRANS. MATH. SOFTWARE, VOL. 12, NO. 3, SEPTEMBER 1986, C PP 265-273. C C***ROUTINES CALLED CACAI,CBKNU,I1MACH,R1MACH C***END PROLOGUE CAIRY COMPLEX AI, CONE, CSQ, CY, S1, S2, TRM1, TRM2, Z, ZTA, Z3 REAL AA, AD, AK, ALIM, ATRM, AZ, AZ3, BK, CK, COEF, C1, C2, DIG, * DK, D1, D2, ELIM, FID, FNU, RL, R1M5, SFAC, TOL, TTH, ZI, ZR, * Z3I, Z3R, R1MACH, BB, ALAZ INTEGER ID, IERR, IFLAG, K, KODE, K1, K2, MR, NN, NZ, I1MACH DIMENSION CY(1) DATA TTH, C1, C2, COEF /6.66666666666666667E-01, * 3.55028053887817240E-01,2.58819403792806799E-01, * 1.83776298473930683E-01/ DATA CONE / (1.0E0,0.0E0) / C***FIRST EXECUTABLE STATEMENT CAIRY IERR = 0 NZ=0 IF (ID.LT.0 .OR. ID.GT.1) IERR=1 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 IF (IERR.NE.0) RETURN AZ = CABS(Z) TOL = AMAX1(R1MACH(4),1.0E-18) FID = FLOAT(ID) IF (AZ.GT.1.0E0) GO TO 60 C----------------------------------------------------------------------- C POWER SERIES FOR CABS(Z).LE.1. C----------------------------------------------------------------------- S1 = CONE S2 = CONE IF (AZ.LT.TOL) GO TO 160 AA = AZ*AZ IF (AA.LT.TOL/AZ) GO TO 40 TRM1 = CONE TRM2 = CONE ATRM = 1.0E0 Z3 = Z*Z*Z AZ3 = AZ*AA AK = 2.0E0 + FID BK = 3.0E0 - FID - FID CK = 4.0E0 - FID DK = 3.0E0 + FID + FID D1 = AK*DK D2 = BK*CK AD = AMIN1(D1,D2) AK = 24.0E0 + 9.0E0*FID BK = 30.0E0 - 9.0E0*FID Z3R = REAL(Z3) Z3I = AIMAG(Z3) DO 30 K=1,25 TRM1 = TRM1*CMPLX(Z3R/D1,Z3I/D1) S1 = S1 + TRM1 TRM2 = TRM2*CMPLX(Z3R/D2,Z3I/D2) S2 = S2 + TRM2 ATRM = ATRM*AZ3/AD D1 = D1 + AK D2 = D2 + BK AD = AMIN1(D1,D2) IF (ATRM.LT.TOL*AD) GO TO 40 AK = AK + 18.0E0 BK = BK + 18.0E0 30 CONTINUE 40 CONTINUE IF (ID.EQ.1) GO TO 50 AI = S1*CMPLX(C1,0.0E0) - Z*S2*CMPLX(C2,0.0E0) IF (KODE.EQ.1) RETURN ZTA = Z*CSQRT(Z)*CMPLX(TTH,0.0E0) AI = AI*CEXP(ZTA) RETURN 50 CONTINUE AI = -S2*CMPLX(C2,0.0E0) IF (AZ.GT.TOL) AI = AI + Z*Z*S1*CMPLX(C1/(1.0E0+FID),0.0E0) IF (KODE.EQ.1) RETURN ZTA = Z*CSQRT(Z)*CMPLX(TTH,0.0E0) AI = AI*CEXP(ZTA) RETURN C----------------------------------------------------------------------- C CASE FOR CABS(Z).GT.1.0 C----------------------------------------------------------------------- 60 CONTINUE FNU = (1.0E0+FID)/3.0E0 C----------------------------------------------------------------------- C SET PARAMETERS RELATED TO MACHINE CONSTANTS. C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). C----------------------------------------------------------------------- K1 = I1MACH(12) K2 = I1MACH(13) R1M5 = R1MACH(5) K = MIN0(IABS(K1),IABS(K2)) ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) K1 = I1MACH(11) - 1 AA = R1M5*FLOAT(K1) DIG = AMIN1(AA,18.0E0) AA = AA*2.303E0 ALIM = ELIM + AMAX1(-AA,-41.45E0) RL = 1.2E0*DIG + 3.0E0 ALAZ=ALOG(AZ) C----------------------------------------------------------------------- C TEST FOR RANGE C----------------------------------------------------------------------- AA=0.5E0/TOL BB=FLOAT(I1MACH(9))*0.5E0 AA=AMIN1(AA,BB) AA=AA**TTH IF (AZ.GT.AA) GO TO 260 AA=SQRT(AA) IF (AZ.GT.AA) IERR=3 CSQ=CSQRT(Z) ZTA=Z*CSQ*CMPLX(TTH,0.0E0) C----------------------------------------------------------------------- C RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL C----------------------------------------------------------------------- IFLAG = 0 SFAC = 1.0E0 ZI = AIMAG(Z) ZR = REAL(Z) AK = AIMAG(ZTA) IF (ZR.GE.0.0E0) GO TO 70 BK = REAL(ZTA) CK = -ABS(BK) ZTA = CMPLX(CK,AK) 70 CONTINUE IF (ZI.NE.0.0E0) GO TO 80 IF (ZR.GT.0.0E0) GO TO 80 ZTA = CMPLX(0.0E0,AK) 80 CONTINUE AA = REAL(ZTA) IF (AA.GE.0.0E0 .AND. ZR.GT.0.0E0) GO TO 100 IF (KODE.EQ.2) GO TO 90 C----------------------------------------------------------------------- C OVERFLOW TEST C----------------------------------------------------------------------- IF (AA.GT.(-ALIM)) GO TO 90 AA = -AA + 0.25E0*ALAZ IFLAG = 1 SFAC = TOL IF (AA.GT.ELIM) GO TO 240 90 CONTINUE C----------------------------------------------------------------------- C CBKNU AND CACAI RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2 C----------------------------------------------------------------------- MR = 1 IF (ZI.LT.0.0E0) MR = -1 CALL CACAI(ZTA, FNU, KODE, MR, 1, CY, NN, RL, TOL, ELIM, ALIM) IF (NN.LT.0) GO TO 250 NZ = NZ + NN GO TO 120 100 CONTINUE IF (KODE.EQ.2) GO TO 110 C----------------------------------------------------------------------- C UNDERFLOW TEST C----------------------------------------------------------------------- IF (AA.LT.ALIM) GO TO 110 AA = -AA - 0.25E0*ALAZ IFLAG = 2 SFAC = 1.0E0/TOL IF (AA.LT.(-ELIM)) GO TO 180 110 CONTINUE CALL CBKNU(ZTA, FNU, KODE, 1, CY, NZ, TOL, ELIM, ALIM) 120 CONTINUE S1 = CY(1)*CMPLX(COEF,0.0E0) IF (IFLAG.NE.0) GO TO 140 IF (ID.EQ.1) GO TO 130 AI = CSQ*S1 RETURN 130 AI = -Z*S1 RETURN 140 CONTINUE S1 = S1*CMPLX(SFAC,0.0E0) IF (ID.EQ.1) GO TO 150 S1 = S1*CSQ AI = S1*CMPLX(1.0E0/SFAC,0.0E0) RETURN 150 CONTINUE S1 = -S1*Z AI = S1*CMPLX(1.0E0/SFAC,0.0E0) RETURN 160 CONTINUE AA = 1.0E+3*R1MACH(1) S1 = CMPLX(0.0E0,0.0E0) IF (ID.EQ.1) GO TO 170 IF (AZ.GT.AA) S1 = CMPLX(C2,0.0E0)*Z AI = CMPLX(C1,0.0E0) - S1 RETURN 170 CONTINUE AI = -CMPLX(C2,0.0E0) AA = SQRT(AA) IF (AZ.GT.AA) S1 = Z*Z*CMPLX(0.5E0,0.0E0) AI = AI + S1*CMPLX(C1,0.0E0) RETURN 180 CONTINUE NZ = 1 AI = CMPLX(0.0E0,0.0E0) RETURN 240 CONTINUE NZ = 0 IERR=2 RETURN 250 CONTINUE IF(NN.EQ.(-1)) GO TO 240 NZ=0 IERR=5 RETURN 260 CONTINUE IERR=4 NZ=0 RETURN END SUBROUTINE CBIRY(Z, ID, KODE, BI, IERR) C***BEGIN PROLOGUE CBIRY C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C***CATEGORY NO. B5K C***KEYWORDS AIRY FUNCTION,BESSEL FUNCTIONS OF ORDER ONE THIRD C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES C***PURPOSE TO COMPUTE AIRY FUNCTIONS BI(Z) AND DBI(Z) FOR COMPLEX Z C***DESCRIPTION C C ON KODE=1, CBIRY COMPUTES THE COMPLEX AIRY FUNCTION BI(Z) OR C ITS DERIVATIVE DBI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON C KODE=2, A SCALING OPTION CEXP(-AXZTA)*BI(Z) OR CEXP(-AXZTA)* C DBI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL BEHAVIOR IN C BOTH THE LEFT AND RIGHT HALF PLANES WHERE C ZTA=(2/3)*Z*CSQRT(Z)=CMPLX(XZTA,YZTA) AND AXZTA=ABS(XZTA). C DEFINITIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF C MATHEMATICAL FUNCTIONS (REF. 1). C C INPUT C Z - Z=CMPLX(X,Y) C ID - ORDER OF DERIVATIVE, ID=0 OR ID=1 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION C KODE= 1 RETURNS C BI=BI(Z) ON ID=0 OR C BI=DBI(Z)/DZ ON ID=1 C = 2 RETURNS C BI=CEXP(-AXZTA)*BI(Z) ON ID=0 OR C BI=CEXP(-AXZTA)*DBI(Z)/DZ ON ID=1 WHERE C ZTA=(2/3)*Z*CSQRT(Z)=CMPLX(XZTA,YZTA) C AND AXZTA=ABS(XZTA) C C OUTPUT C BI - COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND C KODE C IERR - ERROR FLAG C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED C IERR=1, INPUT ERROR - NO COMPUTATION C IERR=2, OVERFLOW - NO COMPUTATION, REAL(Z) C TOO LARGE WITH KODE=1 C IERR=3, CABS(Z) LARGE - COMPUTATION COMPLETED C LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION C PRODUCE LESS THAN HALF OF MACHINE ACCURACY C IERR=4, CABS(Z) TOO LARGE - NO COMPUTATION C COMPLETE LOSS OF ACCURACY BY ARGUMENT C REDUCTION C IERR=5, ERROR - NO COMPUTATION, C ALGORITHM TERMINATION CONDITION NOT MET C C***LONG DESCRIPTION C C BI AND DBI ARE COMPUTED FOR CABS(Z).GT.1.0 FROM THE I BESSEL C FUNCTIONS BY C C BI(Z)=C*SQRT(Z)*( I(-1/3,ZTA) + I(1/3,ZTA) ) C DBI(Z)=C * Z * ( I(-2/3,ZTA) + I(2/3,ZTA) ) C C=1.0/SQRT(3.0) C ZTA=(2/3)*Z**(3/2) C C WITH THE POWER SERIES FOR CABS(Z).LE.1.0. C C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z IS LARGE, LOSSES C OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. CONSEQUENTLY, IF C THE MAGNITUDE OF ZETA=(2/3)*Z**1.5 EXCEEDS U1=SQRT(0.5/UR), C THEN LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR C FLAG IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. C ALSO, IF THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN C ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT C FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE C LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA C MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, C AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE C PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE C PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT- C ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG- C NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN C DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN C EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, C NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE C PRECISION ARITHMETIC. C C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, C OR -PI/2+P. C C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF C COMMERCE, 1955. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 C C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- C 1018, MAY, 1985 C C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, ACM C TRANS. MATH. SOFTWARE, VOL. 12, NO. 3, SEPTEMBER 1986, C PP 265-273. C C***ROUTINES CALLED CBINU,I1MACH,R1MACH C***END PROLOGUE CBIRY COMPLEX BI, CONE, CSQ, CY, S1, S2, TRM1, TRM2, Z, ZTA, Z3 REAL AA, AD, AK, ALIM, ATRM, AZ, AZ3, BB, BK, CK, COEF, C1, C2, * DIG, DK, D1, D2, ELIM, FID, FMR, FNU, FNUL, PI, RL, R1M5, SFAC, * TOL, TTH, ZI, ZR, Z3I, Z3R, R1MACH INTEGER ID, IERR, K, KODE, K1, K2, NZ, I1MACH DIMENSION CY(2) DATA TTH, C1, C2, COEF, PI /6.66666666666666667E-01, * 6.14926627446000736E-01,4.48288357353826359E-01, * 5.77350269189625765E-01,3.14159265358979324E+00/ DATA CONE / (1.0E0,0.0E0) / C***FIRST EXECUTABLE STATEMENT CBIRY IERR = 0 NZ=0 IF (ID.LT.0 .OR. ID.GT.1) IERR=1 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 IF (IERR.NE.0) RETURN AZ = CABS(Z) TOL = AMAX1(R1MACH(4),1.0E-18) FID = FLOAT(ID) IF (AZ.GT.1.0E0) GO TO 60 C----------------------------------------------------------------------- C POWER SERIES FOR CABS(Z).LE.1. C----------------------------------------------------------------------- S1 = CONE S2 = CONE IF (AZ.LT.TOL) GO TO 110 AA = AZ*AZ IF (AA.LT.TOL/AZ) GO TO 40 TRM1 = CONE TRM2 = CONE ATRM = 1.0E0 Z3 = Z*Z*Z AZ3 = AZ*AA AK = 2.0E0 + FID BK = 3.0E0 - FID - FID CK = 4.0E0 - FID DK = 3.0E0 + FID + FID D1 = AK*DK D2 = BK*CK AD = AMIN1(D1,D2) AK = 24.0E0 + 9.0E0*FID BK = 30.0E0 - 9.0E0*FID Z3R = REAL(Z3) Z3I = AIMAG(Z3) DO 30 K=1,25 TRM1 = TRM1*CMPLX(Z3R/D1,Z3I/D1) S1 = S1 + TRM1 TRM2 = TRM2*CMPLX(Z3R/D2,Z3I/D2) S2 = S2 + TRM2 ATRM = ATRM*AZ3/AD D1 = D1 + AK D2 = D2 + BK AD = AMIN1(D1,D2) IF (ATRM.LT.TOL*AD) GO TO 40 AK = AK + 18.0E0 BK = BK + 18.0E0 30 CONTINUE 40 CONTINUE IF (ID.EQ.1) GO TO 50 BI = S1*CMPLX(C1,0.0E0) + Z*S2*CMPLX(C2,0.0E0) IF (KODE.EQ.1) RETURN ZTA = Z*CSQRT(Z)*CMPLX(TTH,0.0E0) AA = REAL(ZTA) AA = -ABS(AA) BI = BI*CMPLX(EXP(AA),0.0E0) RETURN 50 CONTINUE BI = S2*CMPLX(C2,0.0E0) IF (AZ.GT.TOL) BI = BI + Z*Z*S1*CMPLX(C1/(1.0E0+FID),0.0E0) IF (KODE.EQ.1) RETURN ZTA = Z*CSQRT(Z)*CMPLX(TTH,0.0E0) AA = REAL(ZTA) AA = -ABS(AA) BI = BI*CMPLX(EXP(AA),0.0E0) RETURN C----------------------------------------------------------------------- C CASE FOR CABS(Z).GT.1.0 C----------------------------------------------------------------------- 60 CONTINUE FNU = (1.0E0+FID)/3.0E0 C----------------------------------------------------------------------- C SET PARAMETERS RELATED TO MACHINE CONSTANTS. C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU. C----------------------------------------------------------------------- K1 = I1MACH(12) K2 = I1MACH(13) R1M5 = R1MACH(5) K = MIN0(IABS(K1),IABS(K2)) ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) K1 = I1MACH(11) - 1 AA = R1M5*FLOAT(K1) DIG = AMIN1(AA,18.0E0) AA = AA*2.303E0 ALIM = ELIM + AMAX1(-AA,-41.45E0) RL = 1.2E0*DIG + 3.0E0 FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) C----------------------------------------------------------------------- C TEST FOR RANGE C----------------------------------------------------------------------- AA=0.5E0/TOL BB=FLOAT(I1MACH(9))*0.5E0 AA=AMIN1(AA,BB) AA=AA**TTH IF (AZ.GT.AA) GO TO 190 AA=SQRT(AA) IF (AZ.GT.AA) IERR=3 CSQ=CSQRT(Z) ZTA=Z*CSQ*CMPLX(TTH,0.0E0) C----------------------------------------------------------------------- C RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL C----------------------------------------------------------------------- SFAC = 1.0E0 ZI = AIMAG(Z) ZR = REAL(Z) AK = AIMAG(ZTA) IF (ZR.GE.0.0E0) GO TO 70 BK = REAL(ZTA) CK = -ABS(BK) ZTA = CMPLX(CK,AK) 70 CONTINUE IF (ZI.EQ.0.0E0 .AND. ZR.LE.0.0E0) ZTA = CMPLX(0.0E0,AK) AA = REAL(ZTA) IF (KODE.EQ.2) GO TO 80 C----------------------------------------------------------------------- C OVERFLOW TEST C----------------------------------------------------------------------- BB = ABS(AA) IF (BB.LT.ALIM) GO TO 80 BB = BB + 0.25E0*ALOG(AZ) SFAC = TOL IF (BB.GT.ELIM) GO TO 170 80 CONTINUE FMR = 0.0E0 IF (AA.GE.0.0E0 .AND. ZR.GT.0.0E0) GO TO 90 FMR = PI IF (ZI.LT.0.0E0) FMR = -PI ZTA = -ZTA 90 CONTINUE C----------------------------------------------------------------------- C AA=FACTOR FOR ANALYTIC CONTINUATION OF I(FNU,ZTA) C KODE=2 RETURNS EXP(-ABS(XZTA))*I(FNU,ZTA) FROM CBINU C----------------------------------------------------------------------- CALL CBINU(ZTA, FNU, KODE, 1, CY, NZ, RL, FNUL, TOL, ELIM, ALIM) IF (NZ.LT.0) GO TO 180 AA = FMR*FNU Z3 = CMPLX(SFAC,0.0E0) S1 = CY(1)*CMPLX(COS(AA),SIN(AA))*Z3 FNU = (2.0E0-FID)/3.0E0 CALL CBINU(ZTA, FNU, KODE, 2, CY, NZ, RL, FNUL, TOL, ELIM, ALIM) CY(1) = CY(1)*Z3 CY(2) = CY(2)*Z3 C----------------------------------------------------------------------- C BACKWARD RECUR ONE STEP FOR ORDERS -1/3 OR -2/3 C----------------------------------------------------------------------- S2 = CY(1)*CMPLX(FNU+FNU,0.0E0)/ZTA + CY(2) AA = FMR*(FNU-1.0E0) S1 = (S1+S2*CMPLX(COS(AA),SIN(AA)))*CMPLX(COEF,0.0E0) IF (ID.EQ.1) GO TO 100 S1 = CSQ*S1 BI = S1*CMPLX(1.0E0/SFAC,0.0E0) RETURN 100 CONTINUE S1 = Z*S1 BI = S1*CMPLX(1.0E0/SFAC,0.0E0) RETURN 110 CONTINUE AA = C1*(1.0E0-FID) + FID*C2 BI = CMPLX(AA,0.0E0) RETURN 170 CONTINUE NZ=0 IERR=2 RETURN 180 CONTINUE IF(NZ.EQ.(-1)) GO TO 170 NZ=0 IERR=5 RETURN 190 CONTINUE IERR=4 NZ=0 RETURN END SUBROUTINE CUNIK(ZR, FNU, IKFLG, IPMTR, TOL, INIT, PHI, ZETA1, * ZETA2, SUM, CWRK) C***BEGIN PROLOGUE CUNIK C***REFER TO CBESI,CBESK C C CUNIK COMPUTES PARAMETERS FOR THE UNIFORM ASYMPTOTIC C EXPANSIONS OF THE I AND K FUNCTIONS ON IKFLG= 1 OR 2 C RESPECTIVELY BY C C W(FNU,ZR) = PHI*EXP(ZETA)*SUM C C WHERE ZETA=-ZETA1 + ZETA2 OR C ZETA1 - ZETA2 C C THE FIRST CALL MUST HAVE INIT=0. SUBSEQUENT CALLS WITH THE C SAME ZR AND FNU WILL RETURN THE I OR K FUNCTION ON IKFLG= C 1 OR 2 WITH NO CHANGE IN INIT. CWRK IS A COMPLEX WORK C ARRAY. IPMTR=0 COMPUTES ALL PARAMETERS. IPMTR=1 COMPUTES PHI, C ZETA1,ZETA2. C C***ROUTINES CALLED (NONE) C***END PROLOGUE CUNIK COMPLEX CFN, CON, CONE, CRFN, CWRK, CZERO, PHI, S, SR, SUM, T, * T2, ZETA1, ZETA2, ZN, ZR REAL AC, C, FNU, RFN, TEST, TOL, TSTR, TSTI INTEGER I, IKFLG, INIT, IPMTR, J, K, L DIMENSION C(120), CWRK(16), CON(2) DATA CZERO, CONE / (0.0E0,0.0E0), (1.0E0,0.0E0) / DATA CON(1), CON(2) / 1(3.98942280401432678E-01,0.0E0),(1.25331413731550025E+00,0.0E0)/ DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10), 1 C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18), 2 C(19), C(20), C(21), C(22), C(23), C(24)/ 3 1.00000000000000000E+00, -2.08333333333333333E-01, 4 1.25000000000000000E-01, 3.34201388888888889E-01, 5 -4.01041666666666667E-01, 7.03125000000000000E-02, 6 -1.02581259645061728E+00, 1.84646267361111111E+00, 7 -8.91210937500000000E-01, 7.32421875000000000E-02, 8 4.66958442342624743E+00, -1.12070026162229938E+01, 9 8.78912353515625000E+00, -2.36408691406250000E+00, A 1.12152099609375000E-01, -2.82120725582002449E+01, B 8.46362176746007346E+01, -9.18182415432400174E+01, C 4.25349987453884549E+01, -7.36879435947963170E+00, D 2.27108001708984375E-01, 2.12570130039217123E+02, E -7.65252468141181642E+02, 1.05999045252799988E+03/ DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32), 1 C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40), 2 C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/ 3 -6.99579627376132541E+02, 2.18190511744211590E+02, 4 -2.64914304869515555E+01, 5.72501420974731445E-01, 5 -1.91945766231840700E+03, 8.06172218173730938E+03, 6 -1.35865500064341374E+04, 1.16553933368645332E+04, 7 -5.30564697861340311E+03, 1.20090291321635246E+03, 8 -1.08090919788394656E+02, 1.72772750258445740E+00, 9 2.02042913309661486E+04, -9.69805983886375135E+04, A 1.92547001232531532E+05, -2.03400177280415534E+05, B 1.22200464983017460E+05, -4.11926549688975513E+04, C 7.10951430248936372E+03, -4.93915304773088012E+02, D 6.07404200127348304E+00, -2.42919187900551333E+05, E 1.31176361466297720E+06, -2.99801591853810675E+06/ DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56), 1 C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64), 2 C(65), C(66), C(67), C(68), C(69), C(70), C(71), C(72)/ 3 3.76327129765640400E+06, -2.81356322658653411E+06, 4 1.26836527332162478E+06, -3.31645172484563578E+05, 5 4.52187689813627263E+04, -2.49983048181120962E+03, 6 2.43805296995560639E+01, 3.28446985307203782E+06, 7 -1.97068191184322269E+07, 5.09526024926646422E+07, 8 -7.41051482115326577E+07, 6.63445122747290267E+07, 9 -3.75671766607633513E+07, 1.32887671664218183E+07, A -2.78561812808645469E+06, 3.08186404612662398E+05, B -1.38860897537170405E+04, 1.10017140269246738E+02, C -4.93292536645099620E+07, 3.25573074185765749E+08, D -9.39462359681578403E+08, 1.55359689957058006E+09, E -1.62108055210833708E+09, 1.10684281682301447E+09/ DATA C(73), C(74), C(75), C(76), C(77), C(78), C(79), C(80), 1 C(81), C(82), C(83), C(84), C(85), C(86), C(87), C(88), 2 C(89), C(90), C(91), C(92), C(93), C(94), C(95), C(96)/ 3 -4.95889784275030309E+08, 1.42062907797533095E+08, 4 -2.44740627257387285E+07, 2.24376817792244943E+06, 5 -8.40054336030240853E+04, 5.51335896122020586E+02, 6 8.14789096118312115E+08, -5.86648149205184723E+09, 7 1.86882075092958249E+10, -3.46320433881587779E+10, 8 4.12801855797539740E+10, -3.30265997498007231E+10, 9 1.79542137311556001E+10, -6.56329379261928433E+09, A 1.55927986487925751E+09, -2.25105661889415278E+08, B 1.73951075539781645E+07, -5.49842327572288687E+05, C 3.03809051092238427E+03, -1.46792612476956167E+10, D 1.14498237732025810E+11, -3.99096175224466498E+11, E 8.19218669548577329E+11, -1.09837515608122331E+12/ DATA C(97), C(98), C(99), C(100), C(101), C(102), C(103), C(104), 1 C(105), C(106), C(107), C(108), C(109), C(110), C(111), 2 C(112), C(113), C(114), C(115), C(116), C(117), C(118)/ 3 1.00815810686538209E+12, -6.45364869245376503E+11, 4 2.87900649906150589E+11, -8.78670721780232657E+10, 5 1.76347306068349694E+10, -2.16716498322379509E+09, 6 1.43157876718888981E+08, -3.87183344257261262E+06, 7 1.82577554742931747E+04, 2.86464035717679043E+11, 8 -2.40629790002850396E+12, 9.10934118523989896E+12, 9 -2.05168994109344374E+13, 3.05651255199353206E+13, A -3.16670885847851584E+13, 2.33483640445818409E+13, B -1.23204913055982872E+13, 4.61272578084913197E+12, C -1.19655288019618160E+12, 2.05914503232410016E+11, D -2.18229277575292237E+10, 1.24700929351271032E+09/ DATA C(119), C(120)/ 1 -2.91883881222208134E+07, 1.18838426256783253E+05/ C IF (INIT.NE.0) GO TO 40 C----------------------------------------------------------------------- C INITIALIZE ALL VARIABLES C----------------------------------------------------------------------- RFN = 1.0E0/FNU CRFN = CMPLX(RFN,0.0E0) C T = ZR*CRFN C----------------------------------------------------------------------- C OVERFLOW TEST (ZR/FNU TOO SMALL) C----------------------------------------------------------------------- TSTR = REAL(ZR) TSTI = AIMAG(ZR) TEST = R1MACH(1)*1.0E+3 AC = FNU*TEST IF (ABS(TSTR).GT.AC .OR. ABS(TSTI).GT.AC) GO TO 15 AC = 2.0E0*ABS(ALOG(TEST))+FNU ZETA1 = CMPLX(AC,0.0E0) ZETA2 = CMPLX(FNU,0.0E0) PHI=CONE RETURN 15 CONTINUE T=ZR*CRFN S = CONE + T*T SR = CSQRT(S) CFN = CMPLX(FNU,0.0E0) ZN = (CONE+SR)/T ZETA1 = CFN*CLOG(ZN) ZETA2 = CFN*SR T = CONE/SR SR = T*CRFN CWRK(16) = CSQRT(SR) PHI = CWRK(16)*CON(IKFLG) IF (IPMTR.NE.0) RETURN T2 = CONE/S CWRK(1) = CONE CRFN = CONE AC = 1.0E0 L = 1 DO 20 K=2,15 S = CZERO DO 10 J=1,K L = L + 1 S = S*T2 + CMPLX(C(L),0.0E0) 10 CONTINUE CRFN = CRFN*SR CWRK(K) = CRFN*S AC = AC*RFN TSTR = REAL(CWRK(K)) TSTI = AIMAG(CWRK(K)) TEST = ABS(TSTR) + ABS(TSTI) IF (AC.LT.TOL .AND. TEST.LT.TOL) GO TO 30 20 CONTINUE K = 15 30 CONTINUE INIT = K 40 CONTINUE IF (IKFLG.EQ.2) GO TO 60 C----------------------------------------------------------------------- C COMPUTE SUM FOR THE I FUNCTION C----------------------------------------------------------------------- S = CZERO DO 50 I=1,INIT S = S + CWRK(I) 50 CONTINUE SUM = S PHI = CWRK(16)*CON(1) RETURN 60 CONTINUE C----------------------------------------------------------------------- C COMPUTE SUM FOR THE K FUNCTION C----------------------------------------------------------------------- S = CZERO T = CONE DO 70 I=1,INIT S = S + T*CWRK(I) T = -T 70 CONTINUE SUM = S PHI = CWRK(16)*CON(2) RETURN END SUBROUTINE CUOIK(Z, FNU, KODE, IKFLG, N, Y, NUF, TOL, ELIM, ALIM) C***BEGIN PROLOGUE CUOIK C***REFER TO CBESI,CBESK,CBESH C C CUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC C EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM C (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW C WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING C EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN C THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER C MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE C EXP(-ELIM)=SMALLEST MACHINE NUMBER*1.0E+3 AND EXP(-ALIM)= C EXP(-ELIM)/TOL C C IKFLG=1 MEANS THE I SEQUENCE IS TESTED C =2 MEANS THE K SEQUENCE IS TESTED C NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE C =-1 MEANS AN OVERFLOW WOULD OCCUR C IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO C THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE C IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO C IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY C ANOTHER ROUTINE C C***ROUTINES CALLED CUCHK,CUNHJ,CUNIK,R1MACH C***END PROLOGUE CUOIK COMPLEX ARG, ASUM, BSUM, CWRK, CZ, CZERO, PHI, SUM, Y, Z, ZB, * ZETA1, ZETA2, ZN, ZR REAL AARG, AIC, ALIM, APHI, ASCLE, AX, AY, ELIM, FNN, FNU, GNN, * GNU, RCZ, TOL, X, YY INTEGER I, IFORM, IKFLG, INIT, KODE, N, NN, NUF, NW DIMENSION Y(N), CWRK(16) DATA CZERO / (0.0E0,0.0E0) / DATA AIC / 1.265512123484645396E+00 / NUF = 0 NN = N X = REAL(Z) ZR = Z IF (X.LT.0.0E0) ZR = -Z ZB = ZR YY = AIMAG(ZR) AX = ABS(X)*1.7321E0 AY = ABS(YY) IFORM = 1 IF (AY.GT.AX) IFORM = 2 GNU = AMAX1(FNU,1.0E0) IF (IKFLG.EQ.1) GO TO 10 FNN = FLOAT(NN) GNN = FNU + FNN - 1.0E0 GNU = AMAX1(GNN,FNN) 10 CONTINUE C----------------------------------------------------------------------- C ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE C REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET C THE SIGN OF THE IMAGINARY PART CORRECT. C----------------------------------------------------------------------- IF (IFORM.EQ.2) GO TO 20 INIT = 0 CALL CUNIK(ZR, GNU, IKFLG, 1, TOL, INIT, PHI, ZETA1, ZETA2, SUM, * CWRK) CZ = -ZETA1 + ZETA2 GO TO 40 20 CONTINUE ZN = -ZR*CMPLX(0.0E0,1.0E0) IF (YY.GT.0.0E0) GO TO 30 ZN = CONJG(-ZN) 30 CONTINUE CALL CUNHJ(ZN, GNU, 1, TOL, PHI, ARG, ZETA1, ZETA2, ASUM, BSUM) CZ = -ZETA1 + ZETA2 AARG = CABS(ARG) 40 CONTINUE IF (KODE.EQ.2) CZ = CZ - ZB IF (IKFLG.EQ.2) CZ = -CZ APHI = CABS(PHI) RCZ = REAL(CZ) C----------------------------------------------------------------------- C OVERFLOW TEST C----------------------------------------------------------------------- IF (RCZ.GT.ELIM) GO TO 170 IF (RCZ.LT.ALIM) GO TO 50 RCZ = RCZ + ALOG(APHI) IF (IFORM.EQ.2) RCZ = RCZ - 0.25E0*ALOG(AARG) - AIC IF (RCZ.GT.ELIM) GO TO 170 GO TO 100 50 CONTINUE C----------------------------------------------------------------------- C UNDERFLOW TEST C----------------------------------------------------------------------- IF (RCZ.LT.(-ELIM)) GO TO 60 IF (RCZ.GT.(-ALIM)) GO TO 100 RCZ = RCZ + ALOG(APHI) IF (IFORM.EQ.2) RCZ = RCZ - 0.25E0*ALOG(AARG) - AIC IF (RCZ.GT.(-ELIM)) GO TO 80 60 CONTINUE DO 70 I=1,NN Y(I) = CZERO 70 CONTINUE NUF = NN RETURN 80 CONTINUE ASCLE = 1.0E+3*R1MACH(1)/TOL CZ = CZ + CLOG(PHI) IF (IFORM.EQ.1) GO TO 90 CZ = CZ - CMPLX(0.25E0,0.0E0)*CLOG(ARG) - CMPLX(AIC,0.0E0) 90 CONTINUE AX = EXP(RCZ)/TOL AY = AIMAG(CZ) CZ = CMPLX(AX,0.0E0)*CMPLX(COS(AY),SIN(AY)) CALL CUCHK(CZ, NW, ASCLE, TOL) IF (NW.EQ.1) GO TO 60 100 CONTINUE IF (IKFLG.EQ.2) RETURN IF (N.EQ.1) RETURN C----------------------------------------------------------------------- C SET UNDERFLOWS ON I SEQUENCE C----------------------------------------------------------------------- 110 CONTINUE GNU = FNU + FLOAT(NN-1) IF (IFORM.EQ.2) GO TO 120 INIT = 0 CALL CUNIK(ZR, GNU, IKFLG, 1, TOL, INIT, PHI, ZETA1, ZETA2, SUM, * CWRK) CZ = -ZETA1 + ZETA2 GO TO 130 120 CONTINUE CALL CUNHJ(ZN, GNU, 1, TOL, PHI, ARG, ZETA1, ZETA2, ASUM, BSUM) CZ = -ZETA1 + ZETA2 AARG = CABS(ARG) 130 CONTINUE IF (KODE.EQ.2) CZ = CZ - ZB APHI = CABS(PHI) RCZ = REAL(CZ) IF (RCZ.LT.(-ELIM)) GO TO 140 IF (RCZ.GT.(-ALIM)) RETURN RCZ = RCZ + ALOG(APHI) IF (IFORM.EQ.2) RCZ = RCZ - 0.25E0*ALOG(AARG) - AIC IF (RCZ.GT.(-ELIM)) GO TO 150 140 CONTINUE Y(NN) = CZERO NN = NN - 1 NUF = NUF + 1 IF (NN.EQ.0) RETURN GO TO 110 150 CONTINUE ASCLE = 1.0E+3*R1MACH(1)/TOL CZ = CZ + CLOG(PHI) IF (IFORM.EQ.1) GO TO 160 CZ = CZ - CMPLX(0.25E0,0.0E0)*CLOG(ARG) - CMPLX(AIC,0.0E0) 160 CONTINUE AX = EXP(RCZ)/TOL AY = AIMAG(CZ) CZ = CMPLX(AX,0.0E0)*CMPLX(COS(AY),SIN(AY)) CALL CUCHK(CZ, NW, ASCLE, TOL) IF (NW.EQ.1) GO TO 140 RETURN 170 CONTINUE NUF = -1 RETURN END SUBROUTINE CWRSK(ZR, FNU, KODE, N, Y, NZ, CW, TOL, ELIM, ALIM) C***BEGIN PROLOGUE CWRSK C***REFER TO CBESI,CBESK C C CWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY C NORMALIZING THE I FUNCTION RATIOS FROM CRATI BY THE WRONSKIAN C C***ROUTINES CALLED CBKNU,CRATI,R1MACH C***END PROLOGUE CWRSK COMPLEX CINU, CSCL, CT, CW, C1, C2, RCT, ST, Y, ZR REAL ACT, ACW, ALIM, ASCLE, ELIM, FNU, S1, S2, TOL, YY INTEGER I, KODE, N, NW, NZ DIMENSION Y(N), CW(2) C----------------------------------------------------------------------- C I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS C Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE C WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU. C----------------------------------------------------------------------- NZ = 0 CALL CBKNU(ZR, FNU, KODE, 2, CW, NW, TOL, ELIM, ALIM) IF (NW.NE.0) GO TO 50 CALL CRATI(ZR, FNU, N, Y, TOL) C----------------------------------------------------------------------- C RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z), C R(FNU+J-1,Z)=Y(J), J=1,...,N C----------------------------------------------------------------------- CINU = CMPLX(1.0E0,0.0E0) IF (KODE.EQ.1) GO TO 10 YY = AIMAG(ZR) S1 = COS(YY) S2 = SIN(YY) CINU = CMPLX(S1,S2) 10 CONTINUE C----------------------------------------------------------------------- C ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH C THE UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE C SCALED TO PREVENT OVER OR UNDERFLOW. CUOIK HAS DETERMINED THAT C THE RESULT IS ON SCALE. C----------------------------------------------------------------------- ACW = CABS(CW(2)) ASCLE = 1.0E+3*R1MACH(1)/TOL CSCL = CMPLX(1.0E0,0.0E0) IF (ACW.GT.ASCLE) GO TO 20 CSCL = CMPLX(1.0E0/TOL,0.0E0) GO TO 30 20 CONTINUE ASCLE = 1.0E0/ASCLE IF (ACW.LT.ASCLE) GO TO 30 CSCL = CMPLX(TOL,0.0E0) 30 CONTINUE C1 = CW(1)*CSCL C2 = CW(2)*CSCL ST = Y(1) C----------------------------------------------------------------------- C CINU=CINU*(CONJG(CT)/CABS(CT))*(1.0E0/CABS(CT) PREVENTS C UNDER- OR OVERFLOW PREMATURELY BY SQUARING CABS(CT) C----------------------------------------------------------------------- CT = ZR*(C2+ST*C1) ACT = CABS(CT) RCT = CMPLX(1.0E0/ACT,0.0E0) CT = CONJG(CT)*RCT CINU = CINU*RCT*CT Y(1) = CINU*CSCL IF (N.EQ.1) RETURN DO 40 I=2,N CINU = ST*CINU ST = Y(I) Y(I) = CINU*CSCL 40 CONTINUE RETURN 50 CONTINUE NZ = -1 IF(NW.EQ.(-2)) NZ=-2 RETURN END SUBROUTINE CMLRI(Z, FNU, KODE, N, Y, NZ, TOL) C***BEGIN PROLOGUE CMLRI C***REFER TO CBESI,CBESK C C CMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE C MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES. C C***ROUTINES CALLED GAMLN,R1MACH C***END PROLOGUE CMLRI COMPLEX CK, CNORM, CONE, CTWO, CZERO, PT, P1, P2, RZ, SUM, Y, Z REAL ACK, AK, AP, AT, AZ, BK, FKAP, FKK, FLAM, FNF, FNU, RHO, * RHO2, SCLE, TFNF, TOL, TST, X, GAMLN, R1MACH INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N DIMENSION Y(N) DATA CZERO,CONE,CTWO /(0.0E0,0.0E0),(1.0E0,0.0E0),(2.0E0,0.0E0)/ SCLE = 1.0E+3*R1MACH(1)/TOL NZ=0 AZ = CABS(Z) X = REAL(Z) IAZ = INT(AZ) IFNU = INT(FNU) INU = IFNU + N - 1 AT = FLOAT(IAZ) + 1.0E0 CK = CMPLX(AT,0.0E0)/Z RZ = CTWO/Z P1 = CZERO P2 = CONE ACK = (AT+1.0E0)/AZ RHO = ACK + SQRT(ACK*ACK-1.0E0) RHO2 = RHO*RHO TST = (RHO2+RHO2)/((RHO2-1.0E0)*(RHO-1.0E0)) TST = TST/TOL C----------------------------------------------------------------------- C COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES C----------------------------------------------------------------------- AK = AT DO 10 I=1,80 PT = P2 P2 = P1 - CK*P2 P1 = PT CK = CK + RZ AP = CABS(P2) IF (AP.GT.TST*AK*AK) GO TO 20 AK = AK + 1.0E0 10 CONTINUE GO TO 110 20 CONTINUE I = I + 1 K = 0 IF (INU.LT.IAZ) GO TO 40 C----------------------------------------------------------------------- C COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS C----------------------------------------------------------------------- P1 = CZERO P2 = CONE AT = FLOAT(INU) + 1.0E0 CK = CMPLX(AT,0.0E0)/Z ACK = AT/AZ TST = SQRT(ACK/TOL) ITIME = 1 DO 30 K=1,80 PT = P2 P2 = P1 - CK*P2 P1 = PT CK = CK + RZ AP = CABS(P2) IF (AP.LT.TST) GO TO 30 IF (ITIME.EQ.2) GO TO 40 ACK = CABS(CK) FLAM = ACK + SQRT(ACK*ACK-1.0E0) FKAP = AP/CABS(P1) RHO = AMIN1(FLAM,FKAP) TST = TST*SQRT(RHO/(RHO*RHO-1.0E0)) ITIME = 2 30 CONTINUE GO TO 110 40 CONTINUE C----------------------------------------------------------------------- C BACKWARD RECURRENCE AND SUM NORMALIZING RELATION C----------------------------------------------------------------------- K = K + 1 KK = MAX0(I+IAZ,K+INU) FKK = FLOAT(KK) P1 = CZERO C----------------------------------------------------------------------- C SCALE P2 AND SUM BY SCLE C----------------------------------------------------------------------- P2 = CMPLX(SCLE,0.0E0) FNF = FNU - FLOAT(IFNU) TFNF = FNF + FNF BK = GAMLN(FKK+TFNF+1.0E0,IDUM) - GAMLN(FKK+1.0E0,IDUM) * -GAMLN(TFNF+1.0E0,IDUM) BK = EXP(BK) SUM = CZERO KM = KK - INU DO 50 I=1,KM PT = P2 P2 = P1 + CMPLX(FKK+FNF,0.0E0)*RZ*P2 P1 = PT AK = 1.0E0 - TFNF/(FKK+TFNF) ACK = BK*AK SUM = SUM + CMPLX(ACK+BK,0.0E0)*P1 BK = ACK FKK = FKK - 1.0E0 50 CONTINUE Y(N) = P2 IF (N.EQ.1) GO TO 70 DO 60 I=2,N PT = P2 P2 = P1 + CMPLX(FKK+FNF,0.0E0)*RZ*P2 P1 = PT AK = 1.0E0 - TFNF/(FKK+TFNF) ACK = BK*AK SUM = SUM + CMPLX(ACK+BK,0.0E0)*P1 BK = ACK FKK = FKK - 1.0E0 M = N - I + 1 Y(M) = P2 60 CONTINUE 70 CONTINUE IF (IFNU.LE.0) GO TO 90 DO 80 I=1,IFNU PT = P2 P2 = P1 + CMPLX(FKK+FNF,0.0E0)*RZ*P2 P1 = PT AK = 1.0E0 - TFNF/(FKK+TFNF) ACK = BK*AK SUM = SUM + CMPLX(ACK+BK,0.0E0)*P1 BK = ACK FKK = FKK - 1.0E0 80 CONTINUE 90 CONTINUE PT = Z IF (KODE.EQ.2) PT = PT - CMPLX(X,0.0E0) P1 = -CMPLX(FNF,0.0E0)*CLOG(RZ) + PT AP = GAMLN(1.0E0+FNF,IDUM) PT = P1 - CMPLX(AP,0.0E0) C----------------------------------------------------------------------- C THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW C IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES C----------------------------------------------------------------------- P2 = P2 + SUM AP = CABS(P2) P1 = CMPLX(1.0E0/AP,0.0E0) CK = CEXP(PT)*P1 PT = CONJG(P2)*P1 CNORM = CK*PT DO 100 I=1,N Y(I) = Y(I)*CNORM 100 CONTINUE RETURN 110 CONTINUE NZ=-2 RETURN END SUBROUTINE CUNHJ(Z, FNU, IPMTR, TOL, PHI, ARG, ZETA1, ZETA2, * ASUM, BSUM) C***BEGIN PROLOGUE CUNHJ C***REFER TO CBESI,CBESK C C REFERENCES C HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ AND I.A. C STEGUN, AMS55, NATIONAL BUREAU OF STANDARDS, 1965, CHAPTER 9. C C ASYMPTOTICS AND SPECIAL FUNCTIONS BY F.W.J. OLVER, ACADEMIC C PRESS, N.Y., 1974, PAGE 420 C C ABSTRACT C CUNHJ COMPUTES PARAMETERS FOR BESSEL FUNCTIONS C(FNU,Z) = C J(FNU,Z), Y(FNU,Z) OR H(I,FNU,Z) I=1,2 FOR LARGE ORDERS FNU C BY MEANS OF THE UNIFORM ASYMPTOTIC EXPANSION C C C(FNU,Z)=C1*PHI*( ASUM*AIRY(ARG) + C2*BSUM*DAIRY(ARG) ) C C FOR PROPER CHOICES OF C1, C2, AIRY AND DAIRY WHERE AIRY IS C AN AIRY FUNCTION AND DAIRY IS ITS DERIVATIVE. C C (2/3)*FNU*ZETA**1.5 = ZETA1-ZETA2, C C ZETA1=0.5*FNU*CLOG((1+W)/(1-W)), ZETA2=FNU*W FOR SCALING C PURPOSES IN AIRY FUNCTIONS FROM CAIRY OR CBIRY. C C MCONJ=SIGN OF AIMAG(Z), BUT IS AMBIGUOUS WHEN Z IS REAL AND C MUST BE SPECIFIED. IPMTR=0 RETURNS ALL PARAMETERS. IPMTR= C 1 COMPUTES ALL EXCEPT ASUM AND BSUM. C C***ROUTINES CALLED (NONE) C***END PROLOGUE CUNHJ COMPLEX ARG, ASUM, BSUM, CFNU, CONE, CR, CZERO, DR, P, PHI, * PRZTH, PTFN, RFN13, RTZTA, RZTH, SUMA, SUMB, TFN, T2, UP, W, W2, * Z, ZA, ZB, ZC, ZETA, ZETA1, ZETA2, ZTH REAL ALFA, ANG, AP, AR, ATOL, AW2, AZTH, BETA, BR, BTOL, C, EX1, * EX2, FNU, FN13, FN23, GAMA, HPI, PI, PP, RFNU, RFNU2, THPI, TOL, * WI, WR, ZCI, ZCR, ZETAI, ZETAR, ZTHI, ZTHR, ASUMR, ASUMI, BSUMR, * BSUMI, TEST, TSTR, TSTI, AC INTEGER IAS, IBS, IPMTR, IS, J, JR, JU, K, KMAX, KP1, KS, L, LR, * LRP1, L1, L2, M DIMENSION AR(14), BR(14), C(105), ALFA(180), BETA(210), GAMA(30), * AP(30), P(30), UP(14), CR(14), DR(14) DATA AR(1), AR(2), AR(3), AR(4), AR(5), AR(6), AR(7), AR(8), 1 AR(9), AR(10), AR(11), AR(12), AR(13), AR(14)/ 2 1.00000000000000000E+00, 1.04166666666666667E-01, 3 8.35503472222222222E-02, 1.28226574556327160E-01, 4 2.91849026464140464E-01, 8.81627267443757652E-01, 5 3.32140828186276754E+00, 1.49957629868625547E+01, 6 7.89230130115865181E+01, 4.74451538868264323E+02, 7 3.20749009089066193E+03, 2.40865496408740049E+04, 8 1.98923119169509794E+05, 1.79190200777534383E+06/ DATA BR(1), BR(2), BR(3), BR(4), BR(5), BR(6), BR(7), BR(8), 1 BR(9), BR(10), BR(11), BR(12), BR(13), BR(14)/ 2 1.00000000000000000E+00, -1.45833333333333333E-01, 3 -9.87413194444444444E-02, -1.43312053915895062E-01, 4 -3.17227202678413548E-01, -9.42429147957120249E-01, 5 -3.51120304082635426E+00, -1.57272636203680451E+01, 6 -8.22814390971859444E+01, -4.92355370523670524E+02, 7 -3.31621856854797251E+03, -2.48276742452085896E+04, 8 -2.04526587315129788E+05, -1.83844491706820990E+06/ DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10), 1 C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18), 2 C(19), C(20), C(21), C(22), C(23), C(24)/ 3 1.00000000000000000E+00, -2.08333333333333333E-01, 4 1.25000000000000000E-01, 3.34201388888888889E-01, 5 -4.01041666666666667E-01, 7.03125000000000000E-02, 6 -1.02581259645061728E+00, 1.84646267361111111E+00, 7 -8.91210937500000000E-01, 7.32421875000000000E-02, 8 4.66958442342624743E+00, -1.12070026162229938E+01, 9 8.78912353515625000E+00, -2.36408691406250000E+00, A 1.12152099609375000E-01, -2.82120725582002449E+01, B 8.46362176746007346E+01, -9.18182415432400174E+01, C 4.25349987453884549E+01, -7.36879435947963170E+00, D 2.27108001708984375E-01, 2.12570130039217123E+02, E -7.65252468141181642E+02, 1.05999045252799988E+03/ DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32), 1 C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40), 2 C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/ 3 -6.99579627376132541E+02, 2.18190511744211590E+02, 4 -2.64914304869515555E+01, 5.72501420974731445E-01, 5 -1.91945766231840700E+03, 8.06172218173730938E+03, 6 -1.35865500064341374E+04, 1.16553933368645332E+04, 7 -5.30564697861340311E+03, 1.20090291321635246E+03, 8 -1.08090919788394656E+02, 1.72772750258445740E+00, 9 2.02042913309661486E+04, -9.69805983886375135E+04, A 1.92547001232531532E+05, -2.03400177280415534E+05, B 1.22200464983017460E+05, -4.11926549688975513E+04, C 7.10951430248936372E+03, -4.93915304773088012E+02, D 6.07404200127348304E+00, -2.42919187900551333E+05, E 1.31176361466297720E+06, -2.99801591853810675E+06/ DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56), 1 C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64), 2 C(65), C(66), C(67), C(68), C(69), C(70), C(71), C(72)/ 3 3.76327129765640400E+06, -2.81356322658653411E+06, 4 1.26836527332162478E+06, -3.31645172484563578E+05, 5 4.52187689813627263E+04, -2.49983048181120962E+03, 6 2.43805296995560639E+01, 3.28446985307203782E+06, 7 -1.97068191184322269E+07, 5.09526024926646422E+07, 8 -7.41051482115326577E+07, 6.63445122747290267E+07, 9 -3.75671766607633513E+07, 1.32887671664218183E+07, A -2.78561812808645469E+06, 3.08186404612662398E+05, B -1.38860897537170405E+04, 1.10017140269246738E+02, C -4.93292536645099620E+07, 3.25573074185765749E+08, D -9.39462359681578403E+08, 1.55359689957058006E+09, E -1.62108055210833708E+09, 1.10684281682301447E+09/ DATA C(73), C(74), C(75), C(76), C(77), C(78), C(79), C(80), 1 C(81), C(82), C(83), C(84), C(85), C(86), C(87), C(88), 2 C(89), C(90), C(91), C(92), C(93), C(94), C(95), C(96)/ 3 -4.95889784275030309E+08, 1.42062907797533095E+08, 4 -2.44740627257387285E+07, 2.24376817792244943E+06, 5 -8.40054336030240853E+04, 5.51335896122020586E+02, 6 8.14789096118312115E+08, -5.86648149205184723E+09, 7 1.86882075092958249E+10, -3.46320433881587779E+10, 8 4.12801855797539740E+10, -3.30265997498007231E+10, 9 1.79542137311556001E+10, -6.56329379261928433E+09, A 1.55927986487925751E+09, -2.25105661889415278E+08, B 1.73951075539781645E+07, -5.49842327572288687E+05, C 3.03809051092238427E+03, -1.46792612476956167E+10, D 1.14498237732025810E+11, -3.99096175224466498E+11, E 8.19218669548577329E+11, -1.09837515608122331E+12/ DATA C(97), C(98), C(99), C(100), C(101), C(102), C(103), C(104), 1 C(105)/ 2 1.00815810686538209E+12, -6.45364869245376503E+11, 3 2.87900649906150589E+11, -8.78670721780232657E+10, 4 1.76347306068349694E+10, -2.16716498322379509E+09, 5 1.43157876718888981E+08, -3.87183344257261262E+06, 6 1.82577554742931747E+04/ DATA ALFA(1), ALFA(2), ALFA(3), ALFA(4), ALFA(5), ALFA(6), 1 ALFA(7), ALFA(8), ALFA(9), ALFA(10), ALFA(11), ALFA(12), 2 ALFA(13), ALFA(14), ALFA(15), ALFA(16), ALFA(17), ALFA(18), 3 ALFA(19), ALFA(20), ALFA(21), ALFA(22)/ 4 -4.44444444444444444E-03, -9.22077922077922078E-04, 5 -8.84892884892884893E-05, 1.65927687832449737E-04, 6 2.46691372741792910E-04, 2.65995589346254780E-04, 7 2.61824297061500945E-04, 2.48730437344655609E-04, 8 2.32721040083232098E-04, 2.16362485712365082E-04, 9 2.00738858762752355E-04, 1.86267636637545172E-04, A 1.73060775917876493E-04, 1.61091705929015752E-04, B 1.50274774160908134E-04, 1.40503497391269794E-04, C 1.31668816545922806E-04, 1.23667445598253261E-04, D 1.16405271474737902E-04, 1.09798298372713369E-04, E 1.03772410422992823E-04, 9.82626078369363448E-05/ DATA ALFA(23), ALFA(24), ALFA(25), ALFA(26), ALFA(27), ALFA(28), 1 ALFA(29), ALFA(30), ALFA(31), ALFA(32), ALFA(33), ALFA(34), 2 ALFA(35), ALFA(36), ALFA(37), ALFA(38), ALFA(39), ALFA(40), 3 ALFA(41), ALFA(42), ALFA(43), ALFA(44)/ 4 9.32120517249503256E-05, 8.85710852478711718E-05, 5 8.42963105715700223E-05, 8.03497548407791151E-05, 6 7.66981345359207388E-05, 7.33122157481777809E-05, 7 7.01662625163141333E-05, 6.72375633790160292E-05, 8 6.93735541354588974E-04, 2.32241745182921654E-04, 9 -1.41986273556691197E-05, -1.16444931672048640E-04, A -1.50803558053048762E-04, -1.55121924918096223E-04, B -1.46809756646465549E-04, -1.33815503867491367E-04, C -1.19744975684254051E-04, -1.06184319207974020E-04, D -9.37699549891194492E-05, -8.26923045588193274E-05, E -7.29374348155221211E-05, -6.44042357721016283E-05/ DATA ALFA(45), ALFA(46), ALFA(47), ALFA(48), ALFA(49), ALFA(50), 1 ALFA(51), ALFA(52), ALFA(53), ALFA(54), ALFA(55), ALFA(56), 2 ALFA(57), ALFA(58), ALFA(59), ALFA(60), ALFA(61), ALFA(62), 3 ALFA(63), ALFA(64), ALFA(65), ALFA(66)/ 4 -5.69611566009369048E-05, -5.04731044303561628E-05, 5 -4.48134868008882786E-05, -3.98688727717598864E-05, 6 -3.55400532972042498E-05, -3.17414256609022480E-05, 7 -2.83996793904174811E-05, -2.54522720634870566E-05, 8 -2.28459297164724555E-05, -2.05352753106480604E-05, 9 -1.84816217627666085E-05, -1.66519330021393806E-05, A -1.50179412980119482E-05, -1.35554031379040526E-05, B -1.22434746473858131E-05, -1.10641884811308169E-05, C -3.54211971457743841E-04, -1.56161263945159416E-04, D 3.04465503594936410E-05, 1.30198655773242693E-04, E 1.67471106699712269E-04, 1.70222587683592569E-04/ DATA ALFA(67), ALFA(68), ALFA(69), ALFA(70), ALFA(71), ALFA(72), 1 ALFA(73), ALFA(74), ALFA(75), ALFA(76), ALFA(77), ALFA(78), 2 ALFA(79), ALFA(80), ALFA(81), ALFA(82), ALFA(83), ALFA(84), 3 ALFA(85), ALFA(86), ALFA(87), ALFA(88)/ 4 1.56501427608594704E-04, 1.36339170977445120E-04, 5 1.14886692029825128E-04, 9.45869093034688111E-05, 6 7.64498419250898258E-05, 6.07570334965197354E-05, 7 4.74394299290508799E-05, 3.62757512005344297E-05, 8 2.69939714979224901E-05, 1.93210938247939253E-05, 9 1.30056674793963203E-05, 7.82620866744496661E-06, A 3.59257485819351583E-06, 1.44040049814251817E-07, B -2.65396769697939116E-06, -4.91346867098485910E-06, C -6.72739296091248287E-06, -8.17269379678657923E-06, D -9.31304715093561232E-06, -1.02011418798016441E-05, E -1.08805962510592880E-05, -1.13875481509603555E-05/ DATA ALFA(89), ALFA(90), ALFA(91), ALFA(92), ALFA(93), ALFA(94), 1 ALFA(95), ALFA(96), ALFA(97), ALFA(98), ALFA(99), ALFA(100), 2 ALFA(101), ALFA(102), ALFA(103), ALFA(104), ALFA(105), 3 ALFA(106), ALFA(107), ALFA(108), ALFA(109), ALFA(110)/ 4 -1.17519675674556414E-05, -1.19987364870944141E-05, 5 3.78194199201772914E-04, 2.02471952761816167E-04, 6 -6.37938506318862408E-05, -2.38598230603005903E-04, 7 -3.10916256027361568E-04, -3.13680115247576316E-04, 8 -2.78950273791323387E-04, -2.28564082619141374E-04, 9 -1.75245280340846749E-04, -1.25544063060690348E-04, A -8.22982872820208365E-05, -4.62860730588116458E-05, B -1.72334302366962267E-05, 5.60690482304602267E-06, C 2.31395443148286800E-05, 3.62642745856793957E-05, D 4.58006124490188752E-05, 5.24595294959114050E-05, E 5.68396208545815266E-05, 5.94349820393104052E-05/ DATA ALFA(111), ALFA(112), ALFA(113), ALFA(114), ALFA(115), 1 ALFA(116), ALFA(117), ALFA(118), ALFA(119), ALFA(120), 2 ALFA(121), ALFA(122), ALFA(123), ALFA(124), ALFA(125), 3 ALFA(126), ALFA(127), ALFA(128), ALFA(129), ALFA(130)/ 4 6.06478527578421742E-05, 6.08023907788436497E-05, 5 6.01577894539460388E-05, 5.89199657344698500E-05, 6 5.72515823777593053E-05, 5.52804375585852577E-05, 7 5.31063773802880170E-05, 5.08069302012325706E-05, 8 4.84418647620094842E-05, 4.60568581607475370E-05, 9 -6.91141397288294174E-04, -4.29976633058871912E-04, A 1.83067735980039018E-04, 6.60088147542014144E-04, B 8.75964969951185931E-04, 8.77335235958235514E-04, C 7.49369585378990637E-04, 5.63832329756980918E-04, D 3.68059319971443156E-04, 1.88464535514455599E-04/ DATA ALFA(131), ALFA(132), ALFA(133), ALFA(134), ALFA(135), 1 ALFA(136), ALFA(137), ALFA(138), ALFA(139), ALFA(140), 2 ALFA(141), ALFA(142), ALFA(143), ALFA(144), ALFA(145), 3 ALFA(146), ALFA(147), ALFA(148), ALFA(149), ALFA(150)/ 4 3.70663057664904149E-05, -8.28520220232137023E-05, 5 -1.72751952869172998E-04, -2.36314873605872983E-04, 6 -2.77966150694906658E-04, -3.02079514155456919E-04, 7 -3.12594712643820127E-04, -3.12872558758067163E-04, 8 -3.05678038466324377E-04, -2.93226470614557331E-04, 9 -2.77255655582934777E-04, -2.59103928467031709E-04, A -2.39784014396480342E-04, -2.20048260045422848E-04, B -2.00443911094971498E-04, -1.81358692210970687E-04, C -1.63057674478657464E-04, -1.45712672175205844E-04, D -1.29425421983924587E-04, -1.14245691942445952E-04/ DATA ALFA(151), ALFA(152), ALFA(153), ALFA(154), ALFA(155), 1 ALFA(156), ALFA(157), ALFA(158), ALFA(159), ALFA(160), 2 ALFA(161), ALFA(162), ALFA(163), ALFA(164), ALFA(165), 3 ALFA(166), ALFA(167), ALFA(168), ALFA(169), ALFA(170)/ 4 1.92821964248775885E-03, 1.35592576302022234E-03, 5 -7.17858090421302995E-04, -2.58084802575270346E-03, 6 -3.49271130826168475E-03, -3.46986299340960628E-03, 7 -2.82285233351310182E-03, -1.88103076404891354E-03, 8 -8.89531718383947600E-04, 3.87912102631035228E-06, 9 7.28688540119691412E-04, 1.26566373053457758E-03, A 1.62518158372674427E-03, 1.83203153216373172E-03, B 1.91588388990527909E-03, 1.90588846755546138E-03, C 1.82798982421825727E-03, 1.70389506421121530E-03, D 1.55097127171097686E-03, 1.38261421852276159E-03/ DATA ALFA(171), ALFA(172), ALFA(173), ALFA(174), ALFA(175), 1 ALFA(176), ALFA(177), ALFA(178), ALFA(179), ALFA(180)/ 2 1.20881424230064774E-03, 1.03676532638344962E-03, 3 8.71437918068619115E-04, 7.16080155297701002E-04, 4 5.72637002558129372E-04, 4.42089819465802277E-04, 5 3.24724948503090564E-04, 2.20342042730246599E-04, 6 1.28412898401353882E-04, 4.82005924552095464E-05/ DATA BETA(1), BETA(2), BETA(3), BETA(4), BETA(5), BETA(6), 1 BETA(7), BETA(8), BETA(9), BETA(10), BETA(11), BETA(12), 2 BETA(13), BETA(14), BETA(15), BETA(16), BETA(17), BETA(18), 3 BETA(19), BETA(20), BETA(21), BETA(22)/ 4 1.79988721413553309E-02, 5.59964911064388073E-03, 5 2.88501402231132779E-03, 1.80096606761053941E-03, 6 1.24753110589199202E-03, 9.22878876572938311E-04, 7 7.14430421727287357E-04, 5.71787281789704872E-04, 8 4.69431007606481533E-04, 3.93232835462916638E-04, 9 3.34818889318297664E-04, 2.88952148495751517E-04, A 2.52211615549573284E-04, 2.22280580798883327E-04, B 1.97541838033062524E-04, 1.76836855019718004E-04, C 1.59316899661821081E-04, 1.44347930197333986E-04, D 1.31448068119965379E-04, 1.20245444949302884E-04, E 1.10449144504599392E-04, 1.01828770740567258E-04/ DATA BETA(23), BETA(24), BETA(25), BETA(26), BETA(27), BETA(28), 1 BETA(29), BETA(30), BETA(31), BETA(32), BETA(33), BETA(34), 2 BETA(35), BETA(36), BETA(37), BETA(38), BETA(39), BETA(40), 3 BETA(41), BETA(42), BETA(43), BETA(44)/ 4 9.41998224204237509E-05, 8.74130545753834437E-05, 5 8.13466262162801467E-05, 7.59002269646219339E-05, 6 7.09906300634153481E-05, 6.65482874842468183E-05, 7 6.25146958969275078E-05, 5.88403394426251749E-05, 8 -1.49282953213429172E-03, -8.78204709546389328E-04, 9 -5.02916549572034614E-04, -2.94822138512746025E-04, A -1.75463996970782828E-04, -1.04008550460816434E-04, B -5.96141953046457895E-05, -3.12038929076098340E-05, C -1.26089735980230047E-05, -2.42892608575730389E-07, D 8.05996165414273571E-06, 1.36507009262147391E-05, E 1.73964125472926261E-05, 1.98672978842133780E-05/ DATA BETA(45), BETA(46), BETA(47), BETA(48), BETA(49), BETA(50), 1 BETA(51), BETA(52), BETA(53), BETA(54), BETA(55), BETA(56), 2 BETA(57), BETA(58), BETA(59), BETA(60), BETA(61), BETA(62), 3 BETA(6