From 522a6d9b85328d1acb3d9b514268191142fe236f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 20 Mar 2019 12:25:42 +0200 Subject: [PATCH] Replace original netlib amos routines with those in scipy. Note: apparently, one could also use amos from openspecfun https://github.com/JuliaMath/openspecfun Former-commit-id: 77d61fe6486803684201480410021fd173172cc0 --- amos/README.md | 31 +++++ amos/zabs.f | 14 +- amos/zacai.f | 6 +- amos/zacon.f | 8 +- amos/zairy.f | 18 +-- amos/zasyi.f | 14 +- amos/zbesh.f | 6 +- amos/zbesi.f | 269 ++++++++++++++++++++++++++++++++++++ amos/zbesj.f | 4 +- amos/zbesk.f | 281 ++++++++++++++++++++++++++++++++++++++ amos/zbinu.f | 6 +- amos/zbiry.f | 364 +++++++++++++++++++++++++++++++++++++++++++++++++ amos/zbknu.f | 26 ++-- amos/zbuni.f | 8 +- amos/zdiv.f | 6 +- amos/zexp.f | 6 +- amos/zkscl.f | 12 +- amos/zlog.f | 12 +- amos/zmlri.f | 20 +-- amos/zrati.f | 16 +-- amos/zs1s2.f | 14 +- amos/zseri.f | 12 +- amos/zsqrt.f | 12 +- amos/zunhj.f | 20 +-- amos/zuni1.f | 12 +- amos/zuni2.f | 14 +- amos/zunik.f | 8 +- amos/zunk1.f | 18 +-- amos/zunk2.f | 24 ++-- amos/zuoik.f | 20 +-- amos/zwrsk.f | 8 +- 31 files changed, 1117 insertions(+), 172 deletions(-) create mode 100644 amos/README.md create mode 100644 amos/zbesi.f create mode 100644 amos/zbesk.f create mode 100644 amos/zbiry.f diff --git a/amos/README.md b/amos/README.md new file mode 100644 index 0000000..95d7029 --- /dev/null +++ b/amos/README.md @@ -0,0 +1,31 @@ +# AMOS + +A Portable Package for Bessel Functions of a Complex Argument +and Nonnegative Order + +This algorithm is a package of subroutines for computing Bessel +functions and Airy functions. The routines are updated +versions of those routines found in TOMS algorithm 644. + +## Disclaimer + +``` +* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * +* ISSUED BY SANDIA LABORATORIES, +* A PRIME CONTRACTOR TO THE +* UNITED STATES DEPARTMENT OF ENERGY +* * * * * * * * * * * * * * NOTICE * * * * * * * * * * * * * * * +* THIS REPORT WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED BY THE +* UNITED STATES GOVERNMENT. NEITHER THE UNITED STATES NOR THE +* UNITED STATES DEPARTMENT OF ENERGY, NOR ANY OF THEIR +* EMPLOYEES, NOR ANY OF THEIR CONTRACTORS, SUBCONTRACTORS, OR THEIR +* EMPLOYEES, MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY +* LEGAL LIABILITY OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS +* OR USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT OR PROCESS +* DISCLOSED, OR REPRESENTS THAT ITS USE WOULD NOT INFRINGE +* PRIVATELY OWNED RIGHTS. +* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * +* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * +* THIS CODE HAS BEEN APPROVED FOR UNLIMITED RELEASE. +* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * +``` diff --git a/amos/zabs.f b/amos/zabs.f index b25a7ad..31514a2 100644 --- a/amos/zabs.f +++ b/amos/zabs.f @@ -1,12 +1,12 @@ - DOUBLE PRECISION FUNCTION ZABS(ZR, ZI) -C***BEGIN PROLOGUE ZABS + DOUBLE PRECISION FUNCTION AZABS(ZR, ZI) +C***BEGIN PROLOGUE AZABS C***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY C -C ZABS COMPUTES THE ABSOLUTE VALUE OR MAGNITUDE OF A DOUBLE +C AZABS COMPUTES THE ABSOLUTE VALUE OR MAGNITUDE OF A DOUBLE C PRECISION COMPLEX VARIABLE CMPLX(ZR,ZI) C C***ROUTINES CALLED (NONE) -C***END PROLOGUE ZABS +C***END PROLOGUE AZABS DOUBLE PRECISION ZR, ZI, U, V, Q, S U = DABS(ZR) V = DABS(ZI) @@ -19,11 +19,11 @@ C----------------------------------------------------------------------- IF (S.EQ.0.0D+0) GO TO 20 IF (U.GT.V) GO TO 10 Q = U/V - ZABS = V*DSQRT(1.D+0+Q*Q) + AZABS = V*DSQRT(1.D+0+Q*Q) RETURN 10 Q = V/U - ZABS = U*DSQRT(1.D+0+Q*Q) + AZABS = U*DSQRT(1.D+0+Q*Q) RETURN - 20 ZABS = 0.0D+0 + 20 AZABS = 0.0D+0 RETURN END diff --git a/amos/zacai.f b/amos/zacai.f index 345e5eb..87eba88 100644 --- a/amos/zacai.f +++ b/amos/zacai.f @@ -14,19 +14,19 @@ C ZACAI IS THE SAME AS ZACON WITH THE PARTS FOR LARGER ORDERS AND C RECURRENCE REMOVED. A RECURSIVE CALL TO ZACON CAN RESULT IF ZACON C IS CALLED FROM ZAIRY. C -C***ROUTINES CALLED ZASYI,ZBKNU,ZMLRI,ZSERI,ZS1S2,D1MACH,ZABS +C***ROUTINES CALLED ZASYI,ZBKNU,ZMLRI,ZSERI,ZS1S2,D1MACH,AZABS C***END PROLOGUE ZACAI C COMPLEX CSGN,CSPN,C1,C2,Y,Z,ZN,CY DOUBLE PRECISION ALIM, ARG, ASCLE, AZ, CSGNR, CSGNI, CSPNR, * CSPNI, C1R, C1I, C2R, C2I, CYR, CYI, DFNU, ELIM, FMR, FNU, PI, - * RL, SGN, TOL, YY, YR, YI, ZR, ZI, ZNR, ZNI, D1MACH, ZABS + * RL, SGN, TOL, YY, YR, YI, ZR, ZI, ZNR, ZNI, D1MACH, AZABS INTEGER INU, IUF, KODE, MR, N, NN, NW, NZ DIMENSION YR(N), YI(N), CYR(2), CYI(2) DATA PI / 3.14159265358979324D0 / NZ = 0 ZNR = -ZR ZNI = -ZI - AZ = ZABS(ZR,ZI) + AZ = AZABS(ZR,ZI) NN = N DFNU = FNU + DBLE(FLOAT(N-1)) IF (AZ.LE.2.0D0) GO TO 10 diff --git a/amos/zacon.f b/amos/zacon.f index c3964e1..b0dbb91 100644 --- a/amos/zacon.f +++ b/amos/zacon.f @@ -11,7 +11,7 @@ C C TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT C HALF Z PLANE C -C***ROUTINES CALLED ZBINU,ZBKNU,ZS1S2,D1MACH,ZABS,ZMLT +C***ROUTINES CALLED ZBINU,ZBKNU,ZS1S2,D1MACH,AZABS,ZMLT C***END PROLOGUE ZACON C COMPLEX CK,CONE,CSCL,CSCR,CSGN,CSPN,CY,CZERO,C1,C2,RZ,SC1,SC2,ST, C *S1,S2,Y,Z,ZN @@ -20,7 +20,7 @@ C *S1,S2,Y,Z,ZN * CSR, CSRR, CSSR, CYI, CYR, C1I, C1M, C1R, C2I, C2R, ELIM, FMR, * FN, FNU, FNUL, PI, PTI, PTR, RAZN, RL, RZI, RZR, SC1I, SC1R, * SC2I, SC2R, SGN, SPN, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, - * YY, ZEROR, ZI, ZNI, ZNR, ZR, D1MACH, ZABS + * YY, ZEROR, ZI, ZNI, ZNR, ZR, D1MACH, AZABS INTEGER I, INU, IUF, KFLAG, KODE, MR, N, NN, NW, NZ DIMENSION YR(N), YI(N), CYR(2), CYI(2), CSSR(3), CSRR(3), BRY(3) DATA PI / 3.14159265358979324D0 / @@ -102,7 +102,7 @@ C----------------------------------------------------------------------- IF (N.EQ.2) RETURN CSPNR = -CSPNR CSPNI = -CSPNI - AZN = ZABS(ZNR,ZNI) + AZN = AZABS(ZNR,ZNI) RAZN = 1.0D0/AZN STR = ZNR*RAZN STI = -ZNI*RAZN @@ -125,7 +125,7 @@ C----------------------------------------------------------------------- BRY(1) = ASCLE BRY(2) = 1.0D0/ASCLE BRY(3) = D1MACH(2) - AS2 = ZABS(S2R,S2I) + AS2 = AZABS(S2R,S2I) KFLAG = 2 IF (AS2.GT.BRY(1)) GO TO 50 KFLAG = 1 diff --git a/amos/zairy.f b/amos/zairy.f index d6ee764..1563d1d 100644 --- a/amos/zairy.f +++ b/amos/zairy.f @@ -19,7 +19,7 @@ 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 DEFINTIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF +C DEFINITIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF C MATHEMATICAL FUNCTIONS (REF. 1). C C INPUT ZR,ZI ARE DOUBLE PRECISION @@ -124,14 +124,14 @@ C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. C MATH. SOFTWARE, 1986 C -C***ROUTINES CALLED ZACAI,ZBKNU,ZEXP,ZSQRT,I1MACH,D1MACH +C***ROUTINES CALLED ZACAI,ZBKNU,AZEXP,AZSQRT,I1MACH,D1MACH C***END PROLOGUE ZAIRY C COMPLEX AI,CONE,CSQ,CY,S1,S2,TRM1,TRM2,Z,ZTA,Z3 DOUBLE PRECISION AA, AD, AII, AIR, AK, ALIM, ATRM, AZ, AZ3, BK, * CC, CK, COEF, CONEI, CONER, CSQI, CSQR, CYI, CYR, C1, C2, DIG, * DK, D1, D2, ELIM, FID, FNU, PTR, RL, R1M5, SFAC, STI, STR, * S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I, TRM2R, TTH, ZEROI, - * ZEROR, ZI, ZR, ZTAI, ZTAR, Z3I, Z3R, D1MACH, ZABS, ALAZ, BB + * ZEROR, ZI, ZR, ZTAI, ZTAR, Z3I, Z3R, D1MACH, AZABS, ALAZ, BB INTEGER ID, IERR, IFLAG, K, KODE, K1, K2, MR, NN, NZ, I1MACH DIMENSION CYR(1), CYI(1) DATA TTH, C1, C2, COEF /6.66666666666666667D-01, @@ -144,7 +144,7 @@ C***FIRST EXECUTABLE STATEMENT ZAIRY 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 = ZABS(ZR,ZI) + AZ = AZABS(ZR,ZI) TOL = DMAX1(D1MACH(4),1.0D-18) FID = DBLE(FLOAT(ID)) IF (AZ.GT.1.0D0) GO TO 70 @@ -201,10 +201,10 @@ C----------------------------------------------------------------------- AIR = S1R*C1 - C2*(ZR*S2R-ZI*S2I) AII = S1I*C1 - C2*(ZR*S2I+ZI*S2R) IF (KODE.EQ.1) RETURN - CALL ZSQRT(ZR, ZI, STR, STI) + CALL AZSQRT(ZR, ZI, STR, STI) ZTAR = TTH*(ZR*STR-ZI*STI) ZTAI = TTH*(ZR*STI+ZI*STR) - CALL ZEXP(ZTAR, ZTAI, STR, STI) + CALL AZEXP(ZTAR, ZTAI, STR, STI) PTR = AIR*STR - AII*STI AII = AIR*STI + AII*STR AIR = PTR @@ -220,10 +220,10 @@ C----------------------------------------------------------------------- AII = AII + CC*(STR*ZI+STI*ZR) 60 CONTINUE IF (KODE.EQ.1) RETURN - CALL ZSQRT(ZR, ZI, STR, STI) + CALL AZSQRT(ZR, ZI, STR, STI) ZTAR = TTH*(ZR*STR-ZI*STI) ZTAI = TTH*(ZR*STI+ZI*STR) - CALL ZEXP(ZTAR, ZTAI, STR, STI) + CALL AZEXP(ZTAR, ZTAI, STR, STI) PTR = STR*AIR - STI*AII AII = STR*AII + STI*AIR AIR = PTR @@ -265,7 +265,7 @@ C----------------------------------------------------------------------- IF (AZ.GT.AA) GO TO 260 AA=DSQRT(AA) IF (AZ.GT.AA) IERR=3 - CALL ZSQRT(ZR, ZI, CSQR, CSQI) + CALL AZSQRT(ZR, ZI, CSQR, CSQI) ZTAR = TTH*(ZR*CSQR-ZI*CSQI) ZTAI = TTH*(ZR*CSQI+ZI*CSQR) C----------------------------------------------------------------------- diff --git a/amos/zasyi.f b/amos/zasyi.f index 5cd09fa..578136f 100644 --- a/amos/zasyi.f +++ b/amos/zasyi.f @@ -8,21 +8,21 @@ C MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE C REGION CABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN. C NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1. C -C***ROUTINES CALLED D1MACH,ZABS,ZDIV,ZEXP,ZMLT,ZSQRT +C***ROUTINES CALLED D1MACH,AZABS,ZDIV,AZEXP,ZMLT,AZSQRT C***END PROLOGUE ZASYI C COMPLEX AK1,CK,CONE,CS1,CS2,CZ,CZERO,DK,EZ,P1,RZ,S2,Y,Z DOUBLE PRECISION AA, AEZ, AK, AK1I, AK1R, ALIM, ARG, ARM, ATOL, * AZ, BB, BK, CKI, CKR, CONEI, CONER, CS1I, CS1R, CS2I, CS2R, CZI, * CZR, DFNU, DKI, DKR, DNU2, ELIM, EZI, EZR, FDN, FNU, PI, P1I, * P1R, RAZ, RL, RTPI, RTR1, RZI, RZR, S, SGN, SQK, STI, STR, S2I, - * S2R, TOL, TZI, TZR, YI, YR, ZEROI, ZEROR, ZI, ZR, D1MACH, ZABS + * S2R, TOL, TZI, TZR, YI, YR, ZEROI, ZEROR, ZI, ZR, D1MACH, AZABS INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ DIMENSION YR(N), YI(N) DATA PI, RTPI /3.14159265358979324D0 , 0.159154943091895336D0 / DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / C NZ = 0 - AZ = ZABS(ZR,ZI) + AZ = AZABS(ZR,ZI) ARM = 1.0D+3*D1MACH(1) RTR1 = DSQRT(ARM) IL = MIN0(2,N) @@ -35,7 +35,7 @@ C----------------------------------------------------------------------- STI = -ZI*RAZ AK1R = RTPI*STR*RAZ AK1I = RTPI*STI*RAZ - CALL ZSQRT(AK1R, AK1I, AK1R, AK1I) + CALL AZSQRT(AK1R, AK1I, AK1R, AK1I) CZR = ZR CZI = ZI IF (KODE.NE.2) GO TO 10 @@ -47,7 +47,7 @@ C----------------------------------------------------------------------- KODED = 1 IF ((DABS(CZR).GT.ALIM) .AND. (N.GT.2)) GO TO 20 KODED = 0 - CALL ZEXP(CZR, CZI, STR, STI) + CALL AZEXP(CZR, CZI, STR, STI) CALL ZMLT(AK1R, AK1I, STR, STI, AK1R, AK1I) 20 CONTINUE FDN = 0.0D0 @@ -120,7 +120,7 @@ C----------------------------------------------------------------------- IF (ZR+ZR.GE.ELIM) GO TO 60 TZR = ZR + ZR TZI = ZI + ZI - CALL ZEXP(-TZR, -TZI, STR, STI) + CALL AZEXP(-TZR, -TZI, STR, STI) CALL ZMLT(STR, STI, P1R, P1I, STR, STI) CALL ZMLT(STR, STI, CS2R, CS2I, STR, STI) S2R = S2R + STR @@ -149,7 +149,7 @@ C----------------------------------------------------------------------- K = K - 1 80 CONTINUE IF (KODED.EQ.0) RETURN - CALL ZEXP(CZR, CZI, CKR, CKI) + CALL AZEXP(CZR, CZI, CKR, CKI) DO 90 I=1,NN STR = YR(I)*CKR - YI(I)*CKI YI(I) = YR(I)*CKI + YI(I)*CKR diff --git a/amos/zbesh.f b/amos/zbesh.f index cedeb8b..5aacecf 100644 --- a/amos/zbesh.f +++ b/amos/zbesh.f @@ -152,13 +152,13 @@ C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. C MATH. SOFTWARE, 1986 C -C***ROUTINES CALLED ZACON,ZBKNU,ZBUNK,ZUOIK,ZABS,I1MACH,D1MACH +C***ROUTINES CALLED ZACON,ZBKNU,ZBUNK,ZUOIK,AZABS,I1MACH,D1MACH C***END PROLOGUE ZBESH C C COMPLEX CY,Z,ZN,ZT,CSGN DOUBLE PRECISION AA, ALIM, ALN, ARG, AZ, CYI, CYR, DIG, ELIM, * FMM, FN, FNU, FNUL, HPI, RHPI, RL, R1M5, SGN, STR, TOL, UFL, ZI, - * ZNI, ZNR, ZR, ZTI, D1MACH, ZABS, BB, ASCLE, RTOL, ATOL, STI, + * ZNI, ZNR, ZR, ZTI, D1MACH, AZABS, BB, ASCLE, RTOL, ATOL, STI, * CSGNR, CSGNI INTEGER I, IERR, INU, INUH, IR, K, KODE, K1, K2, M, * MM, MR, N, NN, NUF, NW, NZ, I1MACH @@ -208,7 +208,7 @@ C----------------------------------------------------------------------- C----------------------------------------------------------------------- C TEST FOR PROPER RANGE C----------------------------------------------------------------------- - AZ = ZABS(ZR,ZI) + AZ = AZABS(ZR,ZI) AA = 0.5D0/TOL BB=DBLE(FLOAT(I1MACH(9)))*0.5D0 AA = DMIN1(AA,BB) diff --git a/amos/zbesi.f b/amos/zbesi.f new file mode 100644 index 0000000..a2ddd8c --- /dev/null +++ b/amos/zbesi.f @@ -0,0 +1,269 @@ + SUBROUTINE ZBESI(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR) +C***BEGIN PROLOGUE ZBESI +C***DATE WRITTEN 830501 (YYMMDD) +C***REVISION DATE 890801 (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 ***A DOUBLE PRECISION ROUTINE*** +C ON KODE=1, ZBESI 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, ZBESI 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 NOTATION +C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS +C (REF. 1). +C +C INPUT ZR,ZI,FNU ARE DOUBLE PRECISION +C ZR,ZI - Z=CMPLX(ZR,ZI), -PI.LT.ARG(Z).LE.PI +C FNU - ORDER OF INITIAL I FUNCTION, FNU.GE.0.0D0 +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 CYR,CYI ARE DOUBLE PRECISION +C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS +C CONTAIN REAL AND IMAGINARY PARTS 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 TO UNDERFLOW, CY(J)=CMPLX(0.0D0,0.0D0) +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=DMAX1(D1MACH(4),1.0D-18) IS +C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. +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, TRANS. +C MATH. SOFTWARE, 1986 +C +C***ROUTINES CALLED ZBINU,I1MACH,D1MACH +C***END PROLOGUE ZBESI +C COMPLEX CONE,CSGN,CW,CY,CZERO,Z,ZN + DOUBLE PRECISION AA, ALIM, ARG, CONEI, CONER, CSGNI, CSGNR, CYI, + * CYR, DIG, ELIM, FNU, FNUL, PI, RL, R1M5, STR, TOL, ZI, ZNI, ZNR, + * ZR, D1MACH, AZ, BB, FN, AZABS, ASCLE, RTOL, ATOL, STI + INTEGER I, IERR, INU, K, KODE, K1,K2,N,NZ,NN, I1MACH + DIMENSION CYR(N), CYI(N) + DATA PI /3.14159265358979324D0/ + DATA CONER, CONEI /1.0D0,0.0D0/ +C +C***FIRST EXECUTABLE STATEMENT ZBESI + IERR = 0 + NZ=0 + IF (FNU.LT.0.0D0) 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 = DMAX1(D1MACH(4),1.0D-18) + K1 = I1MACH(15) + K2 = I1MACH(16) + R1M5 = D1MACH(5) + K = MIN0(IABS(K1),IABS(K2)) + ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) + K1 = I1MACH(14) - 1 + AA = R1M5*DBLE(FLOAT(K1)) + DIG = DMIN1(AA,18.0D0) + AA = AA*2.303D0 + ALIM = ELIM + DMAX1(-AA,-41.45D0) + RL = 1.2D0*DIG + 3.0D0 + FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) +C----------------------------------------------------------------------------- +C TEST FOR PROPER RANGE +C----------------------------------------------------------------------- + AZ = AZABS(ZR,ZI) + FN = FNU+DBLE(FLOAT(N-1)) + AA = 0.5D0/TOL + BB=DBLE(FLOAT(I1MACH(9)))*0.5D0 + AA = DMIN1(AA,BB) + IF (AZ.GT.AA) GO TO 260 + IF (FN.GT.AA) GO TO 260 + AA = DSQRT(AA) + IF (AZ.GT.AA) IERR=3 + IF (FN.GT.AA) IERR=3 + ZNR = ZR + ZNI = ZI + CSGNR = CONER + CSGNI = CONEI + IF (ZR.GE.0.0D0) GO TO 40 + ZNR = -ZR + ZNI = -ZI +C----------------------------------------------------------------------- +C CALCULATE CSGN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE +C WHEN FNU IS LARGE +C----------------------------------------------------------------------- + INU = INT(SNGL(FNU)) + ARG = (FNU-DBLE(FLOAT(INU)))*PI + IF (ZI.LT.0.0D0) ARG = -ARG + CSGNR = DCOS(ARG) + CSGNI = DSIN(ARG) + IF (MOD(INU,2).EQ.0) GO TO 40 + CSGNR = -CSGNR + CSGNI = -CSGNI + 40 CONTINUE + CALL ZBINU(ZNR, ZNI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL, TOL, + * ELIM, ALIM) + IF (NZ.LT.0) GO TO 120 + IF (ZR.GE.0.0D0) RETURN +C----------------------------------------------------------------------- +C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE +C----------------------------------------------------------------------- + NN = N - NZ + IF (NN.EQ.0) RETURN + RTOL = 1.0D0/TOL + ASCLE = D1MACH(1)*RTOL*1.0D+3 + DO 50 I=1,NN +C STR = CYR(I)*CSGNR - CYI(I)*CSGNI +C CYI(I) = CYR(I)*CSGNI + CYI(I)*CSGNR +C CYR(I) = STR + AA = CYR(I) + BB = CYI(I) + ATOL = 1.0D0 + IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 55 + AA = AA*RTOL + BB = BB*RTOL + ATOL = TOL + 55 CONTINUE + STR = AA*CSGNR - BB*CSGNI + STI = AA*CSGNI + BB*CSGNR + CYR(I) = STR*ATOL + CYI(I) = STI*ATOL + CSGNR = -CSGNR + CSGNI = -CSGNI + 50 CONTINUE + RETURN + 120 CONTINUE + IF(NZ.EQ.(-2)) GO TO 130 + NZ = 0 + IERR=2 + RETURN + 130 CONTINUE + NZ=0 + IERR=5 + RETURN + 260 CONTINUE + NZ=0 + IERR=4 + RETURN + END diff --git a/amos/zbesj.f b/amos/zbesj.f index f3a73de..afe588f 100644 --- a/amos/zbesj.f +++ b/amos/zbesj.f @@ -147,7 +147,7 @@ C C COMPLEX CI,CSGN,CY,Z,ZN DOUBLE PRECISION AA, ALIM, ARG, CII, CSGNI, CSGNR, CYI, CYR, DIG, * ELIM, FNU, FNUL, HPI, RL, R1M5, STR, TOL, ZI, ZNI, ZNR, ZR, - * D1MACH, BB, FN, AZ, ZABS, ASCLE, RTOL, ATOL, STI + * D1MACH, BB, FN, AZ, AZABS, ASCLE, RTOL, ATOL, STI INTEGER I, IERR, INU, INUH, IR, K, KODE, K1, K2, N, NL, NZ, I1MACH DIMENSION CYR(N), CYI(N) DATA HPI /1.57079632679489662D0/ @@ -186,7 +186,7 @@ C----------------------------------------------------------------------- C----------------------------------------------------------------------- C TEST FOR PROPER RANGE C----------------------------------------------------------------------- - AZ = ZABS(ZR,ZI) + AZ = AZABS(ZR,ZI) FN = FNU+DBLE(FLOAT(N-1)) AA = 0.5D0/TOL BB=DBLE(FLOAT(I1MACH(9)))*0.5D0 diff --git a/amos/zbesk.f b/amos/zbesk.f new file mode 100644 index 0000000..cd8eeda --- /dev/null +++ b/amos/zbesk.f @@ -0,0 +1,281 @@ + SUBROUTINE ZBESK(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR) +C***BEGIN PROLOGUE ZBESK +C***DATE WRITTEN 830501 (YYMMDD) +C***REVISION DATE 890801 (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 ***A DOUBLE PRECISION ROUTINE*** +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 ZR,ZI,FNU ARE DOUBLE PRECISION +C ZR,ZI - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0), +C -PI.LT.ARG(Z).LE.PI +C FNU - ORDER OF INITIAL K FUNCTION, FNU.GE.0.0D0 +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 CYR,CYI ARE DOUBLE PRECISION +C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS +C CONTAIN REAL AND IMAGINARY PARTS 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 DUE +C TO UNDERFLOW, CY(I)=CMPLX(0.0D0,0.0D0), +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 +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 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=DMAX1(D1MACH(4),1.0D-18) IS +C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. +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, TRANS. +C MATH. SOFTWARE, 1986 +C +C***ROUTINES CALLED ZACON,ZBKNU,ZBUNK,ZUOIK,AZABS,I1MACH,D1MACH +C***END PROLOGUE ZBESK +C +C COMPLEX CY,Z + DOUBLE PRECISION AA, ALIM, ALN, ARG, AZ, CYI, CYR, DIG, ELIM, FN, + * FNU, FNUL, RL, R1M5, TOL, UFL, ZI, ZR, D1MACH, AZABS, BB + INTEGER IERR, K, KODE, K1, K2, MR, N, NN, NUF, NW, NZ, I1MACH + DIMENSION CYR(N), CYI(N) +C***FIRST EXECUTABLE STATEMENT ZBESK + IERR = 0 + NZ=0 + IF (ZI.EQ.0.0E0 .AND. ZR.EQ.0.0E0) IERR=1 + IF (FNU.LT.0.0D0) 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 = DMAX1(D1MACH(4),1.0D-18) + K1 = I1MACH(15) + K2 = I1MACH(16) + R1M5 = D1MACH(5) + K = MIN0(IABS(K1),IABS(K2)) + ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) + K1 = I1MACH(14) - 1 + AA = R1M5*DBLE(FLOAT(K1)) + DIG = DMIN1(AA,18.0D0) + AA = AA*2.303D0 + ALIM = ELIM + DMAX1(-AA,-41.45D0) + FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) + RL = 1.2D0*DIG + 3.0D0 +C----------------------------------------------------------------------------- +C TEST FOR PROPER RANGE +C----------------------------------------------------------------------- + AZ = AZABS(ZR,ZI) + FN = FNU + DBLE(FLOAT(NN-1)) + AA = 0.5D0/TOL + BB=DBLE(FLOAT(I1MACH(9)))*0.5D0 + AA = DMIN1(AA,BB) + IF (AZ.GT.AA) GO TO 260 + IF (FN.GT.AA) GO TO 260 + AA = DSQRT(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 = DEXP(-ELIM) + UFL = D1MACH(1)*1.0D+3 + IF (AZ.LT.UFL) GO TO 180 + IF (FNU.GT.FNUL) GO TO 80 + IF (FN.LE.1.0D0) GO TO 60 + IF (FN.GT.2.0D0) GO TO 50 + IF (AZ.GT.TOL) GO TO 60 + ARG = 0.5D0*AZ + ALN = -FN*DLOG(ARG) + IF (ALN.GT.ELIM) GO TO 180 + GO TO 60 + 50 CONTINUE + CALL ZUOIK(ZR, ZI, FNU, KODE, 2, NN, CYR, CYI, 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 (ZR.LT.0.0D0) GO TO 70 +C----------------------------------------------------------------------- +C RIGHT HALF PLANE COMPUTATION, REAL(Z).GE.0. +C----------------------------------------------------------------------- + CALL ZBKNU(ZR, ZI, FNU, KODE, NN, CYR, CYI, 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 (ZI.LT.0.0D0) MR = -1 + CALL ZACON(ZR, ZI, FNU, KODE, MR, NN, CYR, CYI, 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 (ZR.GE.0.0D0) GO TO 90 + MR = 1 + IF (ZI.LT.0.0D0) MR = -1 + 90 CONTINUE + CALL ZBUNK(ZR, ZI, FNU, KODE, MR, NN, CYR, CYI, NW, TOL, ELIM, + * ALIM) + IF (NW.LT.0) GO TO 200 + NZ = NZ + NW + RETURN + 100 CONTINUE + IF (ZR.LT.0.0D0) 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 + 260 CONTINUE + NZ=0 + IERR=4 + RETURN + END diff --git a/amos/zbinu.f b/amos/zbinu.f index 7294f41..c76846a 100644 --- a/amos/zbinu.f +++ b/amos/zbinu.f @@ -5,16 +5,16 @@ C***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZAIRY,ZBIRY C C ZBINU COMPUTES THE I FUNCTION IN THE RIGHT HALF Z PLANE C -C***ROUTINES CALLED ZABS,ZASYI,ZBUNI,ZMLRI,ZSERI,ZUOIK,ZWRSK +C***ROUTINES CALLED AZABS,ZASYI,ZBUNI,ZMLRI,ZSERI,ZUOIK,ZWRSK C***END PROLOGUE ZBINU DOUBLE PRECISION ALIM, AZ, CWI, CWR, CYI, CYR, DFNU, ELIM, FNU, - * FNUL, RL, TOL, ZEROI, ZEROR, ZI, ZR, ZABS + * FNUL, RL, TOL, ZEROI, ZEROR, ZI, ZR, AZABS INTEGER I, INW, KODE, N, NLAST, NN, NUI, NW, NZ DIMENSION CYR(N), CYI(N), CWR(2), CWI(2) DATA ZEROR,ZEROI / 0.0D0, 0.0D0 / C NZ = 0 - AZ = ZABS(ZR,ZI) + AZ = AZABS(ZR,ZI) NN = N DFNU = FNU + DBLE(FLOAT(N-1)) IF (AZ.LE.2.0D0) GO TO 10 diff --git a/amos/zbiry.f b/amos/zbiry.f new file mode 100644 index 0000000..e0b6b77 --- /dev/null +++ b/amos/zbiry.f @@ -0,0 +1,364 @@ + SUBROUTINE ZBIRY(ZR, ZI, ID, KODE, BIR, BII, IERR) +C***BEGIN PROLOGUE ZBIRY +C***DATE WRITTEN 830501 (YYMMDD) +C***REVISION DATE 890801 (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 ***A DOUBLE PRECISION ROUTINE*** +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 ZR,ZI ARE DOUBLE PRECISION +C ZR,ZI - Z=CMPLX(ZR,ZI) +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 BIR,BII ARE DOUBLE PRECISION +C BIR,BII- 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 ON 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=DMAX1(D1MACH(4),1.0D-18) IS +C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. +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, TRANS. +C MATH. SOFTWARE, 1986 +C +C***ROUTINES CALLED ZBINU,AZABS,ZDIV,AZSQRT,D1MACH,I1MACH +C***END PROLOGUE ZBIRY +C COMPLEX BI,CONE,CSQ,CY,S1,S2,TRM1,TRM2,Z,ZTA,Z3 + DOUBLE PRECISION AA, AD, AK, ALIM, ATRM, AZ, AZ3, BB, BII, BIR, + * BK, CC, CK, COEF, CONEI, CONER, CSQI, CSQR, CYI, CYR, C1, C2, + * DIG, DK, D1, D2, EAA, ELIM, FID, FMR, FNU, FNUL, PI, RL, R1M5, + * SFAC, STI, STR, S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I, + * TRM2R, TTH, ZI, ZR, ZTAI, ZTAR, Z3I, Z3R, D1MACH, AZABS + INTEGER ID, IERR, K, KODE, K1, K2, NZ, I1MACH + DIMENSION CYR(2), CYI(2) + DATA TTH, C1, C2, COEF, PI /6.66666666666666667D-01, + * 6.14926627446000736D-01,4.48288357353826359D-01, + * 5.77350269189625765D-01,3.14159265358979324D+00/ + DATA CONER, CONEI /1.0D0,0.0D0/ +C***FIRST EXECUTABLE STATEMENT ZBIRY + 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 = AZABS(ZR,ZI) + TOL = DMAX1(D1MACH(4),1.0D-18) + FID = DBLE(FLOAT(ID)) + IF (AZ.GT.1.0E0) GO TO 70 +C----------------------------------------------------------------------- +C POWER SERIES FOR CABS(Z).LE.1. +C----------------------------------------------------------------------- + S1R = CONER + S1I = CONEI + S2R = CONER + S2I = CONEI + IF (AZ.LT.TOL) GO TO 130 + AA = AZ*AZ + IF (AA.LT.TOL/AZ) GO TO 40 + TRM1R = CONER + TRM1I = CONEI + TRM2R = CONER + TRM2I = CONEI + ATRM = 1.0D0 + STR = ZR*ZR - ZI*ZI + STI = ZR*ZI + ZI*ZR + Z3R = STR*ZR - STI*ZI + Z3I = STR*ZI + STI*ZR + AZ3 = AZ*AA + AK = 2.0D0 + FID + BK = 3.0D0 - FID - FID + CK = 4.0D0 - FID + DK = 3.0D0 + FID + FID + D1 = AK*DK + D2 = BK*CK + AD = DMIN1(D1,D2) + AK = 24.0D0 + 9.0D0*FID + BK = 30.0D0 - 9.0D0*FID + DO 30 K=1,25 + STR = (TRM1R*Z3R-TRM1I*Z3I)/D1 + TRM1I = (TRM1R*Z3I+TRM1I*Z3R)/D1 + TRM1R = STR + S1R = S1R + TRM1R + S1I = S1I + TRM1I + STR = (TRM2R*Z3R-TRM2I*Z3I)/D2 + TRM2I = (TRM2R*Z3I+TRM2I*Z3R)/D2 + TRM2R = STR + S2R = S2R + TRM2R + S2I = S2I + TRM2I + ATRM = ATRM*AZ3/AD + D1 = D1 + AK + D2 = D2 + BK + AD = DMIN1(D1,D2) + IF (ATRM.LT.TOL*AD) GO TO 40 + AK = AK + 18.0D0 + BK = BK + 18.0D0 + 30 CONTINUE + 40 CONTINUE + IF (ID.EQ.1) GO TO 50 + BIR = C1*S1R + C2*(ZR*S2R-ZI*S2I) + BII = C1*S1I + C2*(ZR*S2I+ZI*S2R) + IF (KODE.EQ.1) RETURN + CALL AZSQRT(ZR, ZI, STR, STI) + ZTAR = TTH*(ZR*STR-ZI*STI) + ZTAI = TTH*(ZR*STI+ZI*STR) + AA = ZTAR + AA = -DABS(AA) + EAA = DEXP(AA) + BIR = BIR*EAA + BII = BII*EAA + RETURN + 50 CONTINUE + BIR = S2R*C2 + BII = S2I*C2 + IF (AZ.LE.TOL) GO TO 60 + CC = C1/(1.0D0+FID) + STR = S1R*ZR - S1I*ZI + STI = S1R*ZI + S1I*ZR + BIR = BIR + CC*(STR*ZR-STI*ZI) + BII = BII + CC*(STR*ZI+STI*ZR) + 60 CONTINUE + IF (KODE.EQ.1) RETURN + CALL AZSQRT(ZR, ZI, STR, STI) + ZTAR = TTH*(ZR*STR-ZI*STI) + ZTAI = TTH*(ZR*STI+ZI*STR) + AA = ZTAR + AA = -DABS(AA) + EAA = DEXP(AA) + BIR = BIR*EAA + BII = BII*EAA + RETURN +C----------------------------------------------------------------------- +C CASE FOR CABS(Z).GT.1.0 +C----------------------------------------------------------------------- + 70 CONTINUE + FNU = (1.0D0+FID)/3.0D0 +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(15) + K2 = I1MACH(16) + R1M5 = D1MACH(5) + K = MIN0(IABS(K1),IABS(K2)) + ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) + K1 = I1MACH(14) - 1 + AA = R1M5*DBLE(FLOAT(K1)) + DIG = DMIN1(AA,18.0D0) + AA = AA*2.303D0 + ALIM = ELIM + DMAX1(-AA,-41.45D0) + RL = 1.2D0*DIG + 3.0D0 + FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) +C----------------------------------------------------------------------- +C TEST FOR RANGE +C----------------------------------------------------------------------- + AA=0.5D0/TOL + BB=DBLE(FLOAT(I1MACH(9)))*0.5D0 + AA=DMIN1(AA,BB) + AA=AA**TTH + IF (AZ.GT.AA) GO TO 260 + AA=DSQRT(AA) + IF (AZ.GT.AA) IERR=3 + CALL AZSQRT(ZR, ZI, CSQR, CSQI) + ZTAR = TTH*(ZR*CSQR-ZI*CSQI) + ZTAI = TTH*(ZR*CSQI+ZI*CSQR) +C----------------------------------------------------------------------- +C RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL +C----------------------------------------------------------------------- + SFAC = 1.0D0 + AK = ZTAI + IF (ZR.GE.0.0D0) GO TO 80 + BK = ZTAR + CK = -DABS(BK) + ZTAR = CK + ZTAI = AK + 80 CONTINUE + IF (ZI.NE.0.0D0 .OR. ZR.GT.0.0D0) GO TO 90 + ZTAR = 0.0D0 + ZTAI = AK + 90 CONTINUE + AA = ZTAR + IF (KODE.EQ.2) GO TO 100 +C----------------------------------------------------------------------- +C OVERFLOW TEST +C----------------------------------------------------------------------- + BB = DABS(AA) + IF (BB.LT.ALIM) GO TO 100 + BB = BB + 0.25D0*DLOG(AZ) + SFAC = TOL + IF (BB.GT.ELIM) GO TO 190 + 100 CONTINUE + FMR = 0.0D0 + IF (AA.GE.0.0D0 .AND. ZR.GT.0.0D0) GO TO 110 + FMR = PI + IF (ZI.LT.0.0D0) FMR = -PI + ZTAR = -ZTAR + ZTAI = -ZTAI + 110 CONTINUE +C----------------------------------------------------------------------- +C AA=FACTOR FOR ANALYTIC CONTINUATION OF I(FNU,ZTA) +C KODE=2 RETURNS EXP(-ABS(XZTA))*I(FNU,ZTA) FROM CBESI +C----------------------------------------------------------------------- + CALL ZBINU(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, NZ, RL, FNUL, TOL, + * ELIM, ALIM) + IF (NZ.LT.0) GO TO 200 + AA = FMR*FNU + Z3R = SFAC + STR = DCOS(AA) + STI = DSIN(AA) + S1R = (STR*CYR(1)-STI*CYI(1))*Z3R + S1I = (STR*CYI(1)+STI*CYR(1))*Z3R + FNU = (2.0D0-FID)/3.0D0 + CALL ZBINU(ZTAR, ZTAI, FNU, KODE, 2, CYR, CYI, NZ, RL, FNUL, TOL, + * ELIM, ALIM) + CYR(1) = CYR(1)*Z3R + CYI(1) = CYI(1)*Z3R + CYR(2) = CYR(2)*Z3R + CYI(2) = CYI(2)*Z3R +C----------------------------------------------------------------------- +C BACKWARD RECUR ONE STEP FOR ORDERS -1/3 OR -2/3 +C----------------------------------------------------------------------- + CALL ZDIV(CYR(1), CYI(1), ZTAR, ZTAI, STR, STI) + S2R = (FNU+FNU)*STR + CYR(2) + S2I = (FNU+FNU)*STI + CYI(2) + AA = FMR*(FNU-1.0D0) + STR = DCOS(AA) + STI = DSIN(AA) + S1R = COEF*(S1R+S2R*STR-S2I*STI) + S1I = COEF*(S1I+S2R*STI+S2I*STR) + IF (ID.EQ.1) GO TO 120 + STR = CSQR*S1R - CSQI*S1I + S1I = CSQR*S1I + CSQI*S1R + S1R = STR + BIR = S1R/SFAC + BII = S1I/SFAC + RETURN + 120 CONTINUE + STR = ZR*S1R - ZI*S1I + S1I = ZR*S1I + ZI*S1R + S1R = STR + BIR = S1R/SFAC + BII = S1I/SFAC + RETURN + 130 CONTINUE + AA = C1*(1.0D0-FID) + FID*C2 + BIR = AA + BII = 0.0D0 + RETURN + 190 CONTINUE + IERR=2 + NZ=0 + RETURN + 200 CONTINUE + IF(NZ.EQ.(-1)) GO TO 190 + NZ=0 + IERR=5 + RETURN + 260 CONTINUE + IERR=4 + NZ=0 + RETURN + END diff --git a/amos/zbknu.f b/amos/zbknu.f index 4d78d06..a8eb50d 100644 --- a/amos/zbknu.f +++ b/amos/zbknu.f @@ -5,8 +5,8 @@ C***REFER TO ZBESI,ZBESK,ZAIRY,ZBESH C C ZBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE. C -C***ROUTINES CALLED DGAMLN,I1MACH,D1MACH,ZKSCL,ZSHCH,ZUCHK,ZABS,ZDIV, -C ZEXP,ZLOG,ZMLT,ZSQRT +C***ROUTINES CALLED DGAMLN,I1MACH,D1MACH,ZKSCL,ZSHCH,ZUCHK,AZABS,ZDIV, +C AZEXP,AZLOG,ZMLT,AZSQRT C***END PROLOGUE ZBKNU C DOUBLE PRECISION AA, AK, ALIM, ASCLE, A1, A2, BB, BK, BRY, CAZ, @@ -16,7 +16,7 @@ C * FI, FK, FKS, FMUI, FMUR, FNU, FPI, FR, G1, G2, HPI, PI, PR, PTI, * PTR, P1I, P1R, P2I, P2M, P2R, QI, QR, RAK, RCAZ, RTHPI, RZI, * RZR, R1, S, SMUI, SMUR, SPI, STI, STR, S1I, S1R, S2I, S2R, TM, - * TOL, TTH, T1, T2, YI, YR, ZI, ZR, DGAMLN, D1MACH, ZABS, ELM, + * TOL, TTH, T1, T2, YI, YR, ZI, ZR, DGAMLN, D1MACH, AZABS, ELM, * CELMR, ZDR, ZDI, AS, ALAS, HELIM, CYR, CYI INTEGER I, IFLAG, INU, K, KFLAG, KK, KMAX, KODE, KODED, N, NZ, * IDUM, I1MACH, J, IC, INUB, NW @@ -38,7 +38,7 @@ C 3 -2.15241674114950973D-04, -2.01348547807882387D-05, 4 1.13302723198169588D-06, 6.11609510448141582D-09/ C - CAZ = ZABS(ZR,ZI) + CAZ = AZABS(ZR,ZI) CSCLR = 1.0D0/TOL CRSCR = TOL CSSR(1) = CSCLR @@ -68,7 +68,7 @@ C----------------------------------------------------------------------- C SERIES FOR CABS(Z).LE.R1 C----------------------------------------------------------------------- FC = 1.0D0 - CALL ZLOG(RZR, RZI, SMUR, SMUI, IDUM) + CALL AZLOG(RZR, RZI, SMUR, SMUI, IDUM) FMUR = SMUR*DNU FMUI = SMUI*DNU CALL ZSHCH(FMUR, FMUI, CSHR, CSHI, CCHR, CCHI) @@ -104,7 +104,7 @@ C----------------------------------------------------------------------- G2 = (T1+T2)*0.5D0 FR = FC*(CCHR*G1+SMUR*G2) FI = FC*(CCHI*G1+SMUI*G2) - CALL ZEXP(FMUR, FMUI, STR, STI) + CALL AZEXP(FMUR, FMUI, STR, STI) PR = 0.5D0*STR/T2 PI = 0.5D0*STI/T2 CALL ZDIV(0.5D0, 0.0D0, STR, STI, PTR, PTI) @@ -151,7 +151,7 @@ C----------------------------------------------------------------------- YR(1) = S1R YI(1) = S1I IF (KODED.EQ.1) RETURN - CALL ZEXP(ZR, ZI, STR, STI) + CALL AZEXP(ZR, ZI, STR, STI) CALL ZMLT(S1R, S1I, STR, STI, YR(1), YI(1)) RETURN C----------------------------------------------------------------------- @@ -198,7 +198,7 @@ C----------------------------------------------------------------------- S1R = S1R*STR S1I = S1I*STR IF (KODED.EQ.1) GO TO 210 - CALL ZEXP(ZR, ZI, FR, FI) + CALL AZEXP(ZR, ZI, FR, FI) CALL ZMLT(S1R, S1I, FR, FI, S1R, S1I) CALL ZMLT(S2R, S2I, FR, FI, S2R, S2I) GO TO 210 @@ -209,7 +209,7 @@ C KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD C RECURSION C----------------------------------------------------------------------- 110 CONTINUE - CALL ZSQRT(ZR, ZI, STR, STI) + CALL AZSQRT(ZR, ZI, STR, STI) CALL ZDIV(RTHPI, CZEROI, STR, STI, COEFR, COEFI) KFLAG = 2 IF (KODED.EQ.2) GO TO 120 @@ -320,7 +320,7 @@ C----------------------------------------------------------------------- C COMPUTE (P2/CS)=(P2/CABS(CS))*(CONJG(CS)/CABS(CS)) FOR BETTER C SCALING C----------------------------------------------------------------------- - TM = ZABS(CSR,CSI) + TM = AZABS(CSR,CSI) PTR = 1.0D0/TM S1R = P2R*PTR S1I = P2I*PTR @@ -337,7 +337,7 @@ C----------------------------------------------------------------------- C----------------------------------------------------------------------- C COMPUTE P1/P2=(P1/CABS(P2)*CONJG(P2)/CABS(P2) FOR SCALING C----------------------------------------------------------------------- - TM = ZABS(P2R,P2I) + TM = AZABS(P2R,P2I) PTR = 1.0D0/TM P1R = P1R*PTR P1I = P1I*PTR @@ -472,11 +472,11 @@ C----------------------------------------------------------------------- S1I = STI CKR = CKR+RZR CKI = CKI+RZI - AS = ZABS(S2R,S2I) + AS = AZABS(S2R,S2I) ALAS = DLOG(AS) P2R = -ZDR+ALAS IF(P2R.LT.(-ELIM)) GO TO 263 - CALL ZLOG(S2R,S2I,STR,STI,IDUM) + CALL AZLOG(S2R,S2I,STR,STI,IDUM) P2R = -ZDR+STR P2I = -ZDI+STI P2M = DEXP(P2R)/TOL diff --git a/amos/zbuni.f b/amos/zbuni.f index 51cd547..965eddf 100644 --- a/amos/zbuni.f +++ b/amos/zbuni.f @@ -9,12 +9,12 @@ C FNU+N-1 GREATER THAN FNUL BY ADDING NUI AND COMPUTING C ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR I(FNU,Z) C ON IFORM=1 AND THE EXPANSION FOR J(FNU,Z) ON IFORM=2 C -C***ROUTINES CALLED ZUNI1,ZUNI2,ZABS,D1MACH +C***ROUTINES CALLED ZUNI1,ZUNI2,AZABS,D1MACH C***END PROLOGUE ZBUNI C COMPLEX CSCL,CSCR,CY,RZ,ST,S1,S2,Y,Z DOUBLE PRECISION ALIM, AX, AY, CSCLR, CSCRR, CYI, CYR, DFNU, * ELIM, FNU, FNUI, FNUL, GNU, RAZ, RZI, RZR, STI, STR, S1I, S1R, - * S2I, S2R, TOL, YI, YR, ZI, ZR, ZABS, ASCLE, BRY, C1R, C1I, C1M, + * S2I, S2R, TOL, YI, YR, ZI, ZR, AZABS, ASCLE, BRY, C1R, C1I, C1M, * D1MACH INTEGER I, IFLAG, IFORM, K, KODE, N, NL, NLAST, NUI, NW, NZ DIMENSION YR(N), YI(N), CYR(2), CYI(2), BRY(3) @@ -46,7 +46,7 @@ C----------------------------------------------------------------------- 20 CONTINUE IF (NW.LT.0) GO TO 50 IF (NW.NE.0) GO TO 90 - STR = ZABS(CYR(1),CYI(1)) + STR = AZABS(CYR(1),CYI(1)) C---------------------------------------------------------------------- C SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED C---------------------------------------------------------------------- @@ -72,7 +72,7 @@ C---------------------------------------------------------------------- S1I = CYI(2)*CSCLR S2R = CYR(1)*CSCLR S2I = CYI(1)*CSCLR - RAZ = 1.0D0/ZABS(ZR,ZI) + RAZ = 1.0D0/AZABS(ZR,ZI) STR = ZR*RAZ STI = -ZI*RAZ RZR = (STR+STR)*RAZ diff --git a/amos/zdiv.f b/amos/zdiv.f index 1384d8e..f897f4e 100644 --- a/amos/zdiv.f +++ b/amos/zdiv.f @@ -4,11 +4,11 @@ C***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY C C DOUBLE PRECISION COMPLEX DIVIDE C=A/B. C -C***ROUTINES CALLED ZABS +C***ROUTINES CALLED AZABS C***END PROLOGUE ZDIV DOUBLE PRECISION AR, AI, BR, BI, CR, CI, BM, CA, CB, CC, CD - DOUBLE PRECISION ZABS - BM = 1.0D0/ZABS(BR,BI) + DOUBLE PRECISION AZABS + BM = 1.0D0/AZABS(BR,BI) CC = BR*BM CD = BI*BM CA = (AR*CC+AI*CD)*BM diff --git a/amos/zexp.f b/amos/zexp.f index fcb553c..8a34276 100644 --- a/amos/zexp.f +++ b/amos/zexp.f @@ -1,11 +1,11 @@ - SUBROUTINE ZEXP(AR, AI, BR, BI) -C***BEGIN PROLOGUE ZEXP + SUBROUTINE AZEXP(AR, AI, BR, BI) +C***BEGIN PROLOGUE AZEXP C***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY C C DOUBLE PRECISION COMPLEX EXPONENTIAL FUNCTION B=EXP(A) C C***ROUTINES CALLED (NONE) -C***END PROLOGUE ZEXP +C***END PROLOGUE AZEXP DOUBLE PRECISION AR, AI, BR, BI, ZM, CA, CB ZM = DEXP(AR) CA = ZM*DCOS(AI) diff --git a/amos/zkscl.f b/amos/zkscl.f index e051745..90df795 100644 --- a/amos/zkscl.f +++ b/amos/zkscl.f @@ -6,12 +6,12 @@ C SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE C ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN C RETURN WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL. C -C***ROUTINES CALLED ZUCHK,ZABS,ZLOG +C***ROUTINES CALLED ZUCHK,AZABS,AZLOG C***END PROLOGUE ZKSCL C COMPLEX CK,CS,CY,CZERO,RZ,S1,S2,Y,ZR,ZD,CELM DOUBLE PRECISION ACS, AS, ASCLE, CKI, CKR, CSI, CSR, CYI, * CYR, ELIM, FN, FNU, RZI, RZR, STR, S1I, S1R, S2I, - * S2R, TOL, YI, YR, ZEROI, ZEROR, ZRI, ZRR, ZABS, + * S2R, TOL, YI, YR, ZEROI, ZEROR, ZRI, ZRR, AZABS, * ZDR, ZDI, CELMR, ELM, HELIM, ALAS INTEGER I, IC, IDUM, KK, N, NN, NW, NZ DIMENSION YR(N), YI(N), CYR(2), CYI(2) @@ -25,13 +25,13 @@ C S1I = YI(I) CYR(I) = S1R CYI(I) = S1I - AS = ZABS(S1R,S1I) + AS = AZABS(S1R,S1I) ACS = -ZRR + DLOG(AS) NZ = NZ + 1 YR(I) = ZEROR YI(I) = ZEROI IF (ACS.LT.(-ELIM)) GO TO 10 - CALL ZLOG(S1R, S1I, CSR, CSI, IDUM) + CALL AZLOG(S1R, S1I, CSR, CSI, IDUM) CSR = CSR - ZRR CSI = CSI - ZRI STR = DEXP(CSR)/TOL @@ -78,14 +78,14 @@ C S1I = CSI CKR = CKR + RZR CKI = CKI + RZI - AS = ZABS(S2R,S2I) + AS = AZABS(S2R,S2I) ALAS = DLOG(AS) ACS = -ZDR + ALAS NZ = NZ + 1 YR(I) = ZEROR YI(I) = ZEROI IF (ACS.LT.(-ELIM)) GO TO 25 - CALL ZLOG(S2R, S2I, CSR, CSI, IDUM) + CALL AZLOG(S2R, S2I, CSR, CSI, IDUM) CSR = CSR - ZDR CSI = CSI - ZDI STR = DEXP(CSR)/TOL diff --git a/amos/zlog.f b/amos/zlog.f index a9368d5..403d8be 100644 --- a/amos/zlog.f +++ b/amos/zlog.f @@ -1,13 +1,13 @@ - SUBROUTINE ZLOG(AR, AI, BR, BI, IERR) -C***BEGIN PROLOGUE ZLOG + SUBROUTINE AZLOG(AR, AI, BR, BI, IERR) +C***BEGIN PROLOGUE AZLOG C***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY C C DOUBLE PRECISION COMPLEX LOGARITHM B=CLOG(A) C IERR=0,NORMAL RETURN IERR=1, Z=CMPLX(0.0,0.0) -C***ROUTINES CALLED ZABS -C***END PROLOGUE ZLOG +C***ROUTINES CALLED AZABS +C***END PROLOGUE AZLOG DOUBLE PRECISION AR, AI, BR, BI, ZM, DTHETA, DPI, DHPI - DOUBLE PRECISION ZABS + DOUBLE PRECISION AZABS DATA DPI , DHPI / 3.141592653589793238462643383D+0, 1 1.570796326794896619231321696D+0/ C @@ -31,7 +31,7 @@ C BI = 0.0D+0 RETURN 40 IF (AR.LT.0.0D+0) DTHETA = DTHETA + DPI - 50 ZM = ZABS(AR,AI) + 50 ZM = AZABS(AR,AI) BR = DLOG(ZM) BI = DTHETA RETURN diff --git a/amos/zmlri.f b/amos/zmlri.f index 27a5547..5bc8d77 100644 --- a/amos/zmlri.f +++ b/amos/zmlri.f @@ -5,20 +5,20 @@ C C ZMLRI 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 DGAMLN,D1MACH,ZABS,ZEXP,ZLOG,ZMLT +C***ROUTINES CALLED DGAMLN,D1MACH,AZABS,AZEXP,AZLOG,ZMLT C***END PROLOGUE ZMLRI C COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z DOUBLE PRECISION ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI, * CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, FNU, PTI, PTR, P1I, * P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI, * SUMR, TFNF, TOL, TST, YI, YR, ZEROI, ZEROR, ZI, ZR, DGAMLN, - * D1MACH, ZABS + * D1MACH, AZABS INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N, NZ DIMENSION YR(N), YI(N) DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / SCLE = D1MACH(1)/TOL NZ=0 - AZ = ZABS(ZR,ZI) + AZ = AZABS(ZR,ZI) IAZ = INT(SNGL(AZ)) IFNU = INT(SNGL(FNU)) INU = IFNU + N - 1 @@ -52,7 +52,7 @@ C----------------------------------------------------------------------- P1I = PTI CKR = CKR + RZR CKI = CKI + RZI - AP = ZABS(P2R,P2I) + AP = AZABS(P2R,P2I) IF (AP.GT.TST*AK*AK) GO TO 20 AK = AK + 1.0D0 10 CONTINUE @@ -85,12 +85,12 @@ C----------------------------------------------------------------------- P1I = PTI CKR = CKR + RZR CKI = CKI + RZI - AP = ZABS(P2R,P2I) + AP = AZABS(P2R,P2I) IF (AP.LT.TST) GO TO 30 IF (ITIME.EQ.2) GO TO 40 - ACK = ZABS(CKR,CKI) + ACK = AZABS(CKR,CKI) FLAM = ACK + DSQRT(ACK*ACK-1.0D0) - FKAP = AP/ZABS(P1R,P1I) + FKAP = AP/AZABS(P1R,P1I) RHO = DMIN1(FLAM,FKAP) TST = TST*DSQRT(RHO/(RHO*RHO-1.0D0)) ITIME = 2 @@ -172,7 +172,7 @@ C----------------------------------------------------------------------- PTR = ZR PTI = ZI IF (KODE.EQ.2) PTR = ZEROR - CALL ZLOG(RZR, RZI, STR, STI, IDUM) + CALL AZLOG(RZR, RZI, STR, STI, IDUM) P1R = -FNF*STR + PTR P1I = -FNF*STI + PTI AP = DGAMLN(1.0D0+FNF,IDUM) @@ -184,9 +184,9 @@ C IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES C----------------------------------------------------------------------- P2R = P2R + SUMR P2I = P2I + SUMI - AP = ZABS(P2R,P2I) + AP = AZABS(P2R,P2I) P1R = 1.0D0/AP - CALL ZEXP(PTR, PTI, STR, STI) + CALL AZEXP(PTR, PTI, STR, STI) CKR = STR*P1R CKI = STI*P1R PTR = P2R*P1R diff --git a/amos/zrati.f b/amos/zrati.f index 0458918..d8ab777 100644 --- a/amos/zrati.f +++ b/amos/zrati.f @@ -9,18 +9,18 @@ C MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973, C BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER, C BY D. J. SOOKNE. C -C***ROUTINES CALLED ZABS,ZDIV +C***ROUTINES CALLED AZABS,ZDIV C***END PROLOGUE ZRATI C COMPLEX Z,CY(1),CONE,CZERO,P1,P2,T1,RZ,PT,CDFNU DOUBLE PRECISION AK, AMAGZ, AP1, AP2, ARG, AZ, CDFNUI, CDFNUR, * CONEI, CONER, CYI, CYR, CZEROI, CZEROR, DFNU, FDNU, FLAM, FNU, * FNUP, PTI, PTR, P1I, P1R, P2I, P2R, RAK, RAP1, RHO, RT2, RZI, - * RZR, TEST, TEST1, TOL, TTI, TTR, T1I, T1R, ZI, ZR, ZABS + * RZR, TEST, TEST1, TOL, TTI, TTR, T1I, T1R, ZI, ZR, AZABS INTEGER I, ID, IDNU, INU, ITIME, K, KK, MAGZ, N DIMENSION CYR(N), CYI(N) DATA CZEROR,CZEROI,CONER,CONEI,RT2/ 1 0.0D0, 0.0D0, 1.0D0, 0.0D0, 1.41421356237309505D0 / - AZ = ZABS(ZR,ZI) + AZ = AZABS(ZR,ZI) INU = INT(SNGL(FNU)) IDNU = INU + N - 1 MAGZ = INT(SNGL(AZ)) @@ -42,8 +42,8 @@ C COMPLEX Z,CY(1),CONE,CZERO,P1,P2,T1,RZ,PT,CDFNU T1R = T1R + RZR T1I = T1I + RZI IF (ID.GT.0) ID = 0 - AP2 = ZABS(P2R,P2I) - AP1 = ZABS(P1R,P1I) + AP2 = AZABS(P2R,P2I) + AP1 = AZABS(P1R,P1I) C----------------------------------------------------------------------- C THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU C GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT @@ -70,10 +70,10 @@ C----------------------------------------------------------------------- P1I = PTI T1R = T1R + RZR T1I = T1I + RZI - AP2 = ZABS(P2R,P2I) + AP2 = AZABS(P2R,P2I) IF (AP1.LE.TEST) GO TO 10 IF (ITIME.EQ.2) GO TO 20 - AK = ZABS(T1R,T1I)*0.5D0 + AK = AZABS(T1R,T1I)*0.5D0 FLAM = AK + DSQRT(AK*AK-1.0D0) RHO = DMIN1(AP2/AP1,FLAM) TEST = TEST1*DSQRT(RHO/(RHO*RHO-1.0D0)) @@ -116,7 +116,7 @@ C----------------------------------------------------------------------- DO 60 I=2,N PTR = CDFNUR + (T1R*RZR-T1I*RZI) + CYR(K+1) PTI = CDFNUI + (T1R*RZI+T1I*RZR) + CYI(K+1) - AK = ZABS(PTR,PTI) + AK = AZABS(PTR,PTI) IF (AK.NE.CZEROR) GO TO 50 PTR = TOL PTI = TOL diff --git a/amos/zs1s2.f b/amos/zs1s2.f index 749e57e..194be44 100644 --- a/amos/zs1s2.f +++ b/amos/zs1s2.f @@ -11,16 +11,16 @@ C MAGNITUDE, BUT FOR KODE=2 THEY CAN BE OF THE SAME ORDER C OF MAGNITUDE AND THE MAXIMUM MUST BE AT LEAST ONE C PRECISION ABOVE THE UNDERFLOW LIMIT. C -C***ROUTINES CALLED ZABS,ZEXP,ZLOG +C***ROUTINES CALLED AZABS,AZEXP,AZLOG C***END PROLOGUE ZS1S2 C COMPLEX CZERO,C1,S1,S1D,S2,ZR DOUBLE PRECISION AA, ALIM, ALN, ASCLE, AS1, AS2, C1I, C1R, S1DI, - * S1DR, S1I, S1R, S2I, S2R, ZEROI, ZEROR, ZRI, ZRR, ZABS + * S1DR, S1I, S1R, S2I, S2R, ZEROI, ZEROR, ZRI, ZRR, AZABS INTEGER IUF, IDUM, NZ DATA ZEROR,ZEROI / 0.0D0 , 0.0D0 / NZ = 0 - AS1 = ZABS(S1R,S1I) - AS2 = ZABS(S2R,S2I) + AS1 = AZABS(S1R,S1I) + AS2 = AZABS(S2R,S2I) IF (S1R.EQ.0.0D0 .AND. S1I.EQ.0.0D0) GO TO 10 IF (AS1.EQ.0.0D0) GO TO 10 ALN = -ZRR - ZRR + DLOG(AS1) @@ -30,11 +30,11 @@ C COMPLEX CZERO,C1,S1,S1D,S2,ZR S1I = ZEROI AS1 = ZEROR IF (ALN.LT.(-ALIM)) GO TO 10 - CALL ZLOG(S1DR, S1DI, C1R, C1I, IDUM) + CALL AZLOG(S1DR, S1DI, C1R, C1I, IDUM) C1R = C1R - ZRR - ZRR C1I = C1I - ZRI - ZRI - CALL ZEXP(C1R, C1I, S1R, S1I) - AS1 = ZABS(S1R,S1I) + CALL AZEXP(C1R, C1I, S1R, S1I) + AS1 = AZABS(S1R,S1I) IUF = IUF + 1 10 CONTINUE AA = DMAX1(AS1,AS2) diff --git a/amos/zseri.f b/amos/zseri.f index 2e8e14a..4c487f0 100644 --- a/amos/zseri.f +++ b/amos/zseri.f @@ -11,20 +11,20 @@ C DUE TO UNDERFLOW. NZ.LT.0 MEANS UNDERFLOW OCCURRED, BUT THE C CONDITION CABS(Z).LE.2*SQRT(FNU+1) WAS VIOLATED AND THE C COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ). C -C***ROUTINES CALLED DGAMLN,D1MACH,ZUCHK,ZABS,ZDIV,ZLOG,ZMLT +C***ROUTINES CALLED DGAMLN,D1MACH,ZUCHK,AZABS,ZDIV,AZLOG,ZMLT C***END PROLOGUE ZSERI C COMPLEX AK1,CK,COEF,CONE,CRSC,CSCL,CZ,CZERO,HZ,RZ,S1,S2,Y,Z DOUBLE PRECISION AA, ACZ, AK, AK1I, AK1R, ALIM, ARM, ASCLE, ATOL, * AZ, CKI, CKR, COEFI, COEFR, CONEI, CONER, CRSCR, CZI, CZR, DFNU, * ELIM, FNU, FNUP, HZI, HZR, RAZ, RS, RTR1, RZI, RZR, S, SS, STI, * STR, S1I, S1R, S2I, S2R, TOL, YI, YR, WI, WR, ZEROI, ZEROR, ZI, - * ZR, DGAMLN, D1MACH, ZABS + * ZR, DGAMLN, D1MACH, AZABS INTEGER I, IB, IDUM, IFLAG, IL, K, KODE, L, M, N, NN, NZ, NW DIMENSION YR(N), YI(N), WR(2), WI(2) DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / C NZ = 0 - AZ = ZABS(ZR,ZI) + AZ = AZABS(ZR,ZI) IF (AZ.EQ.0.0D0) GO TO 160 ARM = 1.0D+3*D1MACH(1) RTR1 = DSQRT(ARM) @@ -38,9 +38,9 @@ C IF (AZ.LE.RTR1) GO TO 10 CALL ZMLT(HZR, HZI, HZR, HZI, CZR, CZI) 10 CONTINUE - ACZ = ZABS(CZR,CZI) + ACZ = AZABS(CZR,CZI) NN = N - CALL ZLOG(HZR, HZI, CKR, CKI, IDUM) + CALL AZLOG(HZR, HZI, CKR, CKI, IDUM) 20 CONTINUE DFNU = FNU + DBLE(FLOAT(NN-1)) FNUP = DFNU + 1.0D0 @@ -157,7 +157,7 @@ C----------------------------------------------------------------------- YI(K) = CKI AK = AK - 1.0D0 K = K - 1 - IF (ZABS(CKR,CKI).GT.ASCLE) GO TO 140 + IF (AZABS(CKR,CKI).GT.ASCLE) GO TO 140 130 CONTINUE RETURN 140 CONTINUE diff --git a/amos/zsqrt.f b/amos/zsqrt.f index 08c17cc..43ba617 100644 --- a/amos/zsqrt.f +++ b/amos/zsqrt.f @@ -1,16 +1,16 @@ - SUBROUTINE ZSQRT(AR, AI, BR, BI) -C***BEGIN PROLOGUE ZSQRT + SUBROUTINE AZSQRT(AR, AI, BR, BI) +C***BEGIN PROLOGUE AZSQRT C***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY C C DOUBLE PRECISION COMPLEX SQUARE ROOT, B=CSQRT(A) C -C***ROUTINES CALLED ZABS -C***END PROLOGUE ZSQRT +C***ROUTINES CALLED AZABS +C***END PROLOGUE AZSQRT DOUBLE PRECISION AR, AI, BR, BI, ZM, DTHETA, DPI, DRT - DOUBLE PRECISION ZABS + DOUBLE PRECISION AZABS DATA DRT , DPI / 7.071067811865475244008443621D-1, 1 3.141592653589793238462643383D+0/ - ZM = ZABS(AR,AI) + ZM = AZABS(AR,AI) ZM = DSQRT(ZM) IF (AR.EQ.0.0D+0) GO TO 10 IF (AI.EQ.0.0D+0) GO TO 20 diff --git a/amos/zunhj.f b/amos/zunhj.f index 41ecc24..e56e43c 100644 --- a/amos/zunhj.f +++ b/amos/zunhj.f @@ -29,7 +29,7 @@ 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 ZABS,ZDIV,ZLOG,ZSQRT,D1MACH +C***ROUTINES CALLED AZABS,ZDIV,AZLOG,AZSQRT,D1MACH C***END PROLOGUE ZUNHJ C COMPLEX ARG,ASUM,BSUM,CFNU,CONE,CR,CZERO,DR,P,PHI,PRZTH,PTFN, C *RFN13,RTZTA,RZTH,SUMA,SUMB,TFN,T2,UP,W,W2,Z,ZA,ZB,ZC,ZETA,ZETA1, @@ -42,7 +42,7 @@ C *ZETA2,ZTH * SUMAI, SUMAR, SUMBI, SUMBR, TEST, TFNI, TFNR, THPI, TOL, TZAI, * TZAR, T2I, T2R, UPI, UPR, WI, WR, W2I, W2R, ZAI, ZAR, ZBI, ZBR, * ZCI, ZCR, ZEROI, ZEROR, ZETAI, ZETAR, ZETA1I, ZETA1R, ZETA2I, - * ZETA2R, ZI, ZR, ZTHI, ZTHR, ZABS, AC, D1MACH + * ZETA2R, ZI, ZR, ZTHI, ZTHR, AZABS, AC, D1MACH INTEGER IAS, IBS, IPMTR, IS, J, JR, JU, K, KMAX, KP1, KS, L, LR, * LRP1, L1, L2, M, IDUM DIMENSION AR(14), BR(14), C(105), ALFA(180), BETA(210), GAMA(30), @@ -457,7 +457,7 @@ C----------------------------------------------------------------------- RFN13 = 1.0D0/FN13 W2R = CONER - ZBR*ZBR + ZBI*ZBI W2I = CONEI - ZBR*ZBI - ZBR*ZBI - AW2 = ZABS(W2R,W2I) + AW2 = AZABS(W2R,W2I) IF (AW2.GT.0.25D0) GO TO 130 C----------------------------------------------------------------------- C POWER SERIES FOR CABS(W2).LE.0.25D0 @@ -484,8 +484,8 @@ C----------------------------------------------------------------------- ZETAI = W2R*SUMAI + W2I*SUMAR ARGR = ZETAR*FN23 ARGI = ZETAI*FN23 - CALL ZSQRT(SUMAR, SUMAI, ZAR, ZAI) - CALL ZSQRT(W2R, W2I, STR, STI) + CALL AZSQRT(SUMAR, SUMAI, ZAR, ZAI) + CALL AZSQRT(W2R, W2I, STR, STI) ZETA2R = STR*FNU ZETA2I = STI*FNU STR = CONER + EX2*(ZETAR*ZAR-ZETAI*ZAI) @@ -494,7 +494,7 @@ C----------------------------------------------------------------------- ZETA1I = STR*ZETA2I + STI*ZETA2R ZAR = ZAR + ZAR ZAI = ZAI + ZAI - CALL ZSQRT(ZAR, ZAI, STR, STI) + CALL AZSQRT(ZAR, ZAI, STR, STI) PHIR = STR*RFN13 PHII = STI*RFN13 IF (IPMTR.EQ.1) GO TO 120 @@ -565,13 +565,13 @@ C----------------------------------------------------------------------- C CABS(W2).GT.0.25D0 C----------------------------------------------------------------------- 130 CONTINUE - CALL ZSQRT(W2R, W2I, WR, WI) + CALL AZSQRT(W2R, W2I, WR, WI) IF (WR.LT.0.0D0) WR = 0.0D0 IF (WI.LT.0.0D0) WI = 0.0D0 STR = CONER + WR STI = WI CALL ZDIV(STR, STI, ZBR, ZBI, ZAR, ZAI) - CALL ZLOG(ZAR, ZAI, ZCR, ZCI, IDUM) + CALL AZLOG(ZAR, ZAI, ZCR, ZCI, IDUM) IF (ZCI.LT.0.0D0) ZCI = 0.0D0 IF (ZCI.GT.HPI) ZCI = HPI IF (ZCR.LT.0.0D0) ZCR = 0.0D0 @@ -581,7 +581,7 @@ C----------------------------------------------------------------------- ZETA1I = ZCI*FNU ZETA2R = WR*FNU ZETA2I = WI*FNU - AZTH = ZABS(ZTHR,ZTHI) + AZTH = AZABS(ZTHR,ZTHI) ANG = THPI IF (ZTHR.GE.0.0D0 .AND. ZTHI.LT.0.0D0) GO TO 140 ANG = HPI @@ -600,7 +600,7 @@ C----------------------------------------------------------------------- CALL ZDIV(RTZTR, RTZTI, WR, WI, ZAR, ZAI) TZAR = ZAR + ZAR TZAI = ZAI + ZAI - CALL ZSQRT(TZAR, TZAI, STR, STI) + CALL AZSQRT(TZAR, TZAI, STR, STI) PHIR = STR*RFN13 PHII = STI*RFN13 IF (IPMTR.EQ.1) GO TO 120 diff --git a/amos/zuni1.f b/amos/zuni1.f index 0e6ab5d..c7173b3 100644 --- a/amos/zuni1.f +++ b/amos/zuni1.f @@ -12,7 +12,7 @@ C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL. C Y(I)=CZERO FOR I=NLAST+1,N C -C***ROUTINES CALLED ZUCHK,ZUNIK,ZUOIK,D1MACH,ZABS +C***ROUTINES CALLED ZUCHK,ZUNIK,ZUOIK,D1MACH,AZABS C***END PROLOGUE ZUNI1 C COMPLEX CFN,CONE,CRSC,CSCL,CSR,CSS,CWRK,CZERO,C1,C2,PHI,RZ,SUM,S1, C *S2,Y,Z,ZETA1,ZETA2 @@ -20,7 +20,7 @@ C *S2,Y,Z,ZETA1,ZETA2 * CSCL, CSRR, CSSR, CWRKI, CWRKR, C1R, C2I, C2M, C2R, ELIM, FN, * FNU, FNUL, PHII, PHIR, RAST, RS1, RZI, RZR, STI, STR, SUMI, * SUMR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZEROI, ZEROR, ZETA1I, - * ZETA1R, ZETA2I, ZETA2R, ZI, ZR, CYR, CYI, D1MACH, ZABS + * ZETA1R, ZETA2I, ZETA2R, ZI, ZR, CYR, CYI, D1MACH, AZABS INTEGER I, IFLAG, INIT, K, KODE, M, N, ND, NLAST, NN, NUF, NW, NZ DIMENSION BRY(3), YR(N), YI(N), CWRKR(16), CWRKI(16), CSSR(3), * CSRR(3), CYR(2), CYI(2) @@ -53,7 +53,7 @@ C----------------------------------------------------------------------- IF (KODE.EQ.1) GO TO 10 STR = ZR + ZETA2R STI = ZI + ZETA2I - RAST = FN/ZABS(STR,STI) + RAST = FN/AZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = -ZETA1R + STR @@ -75,7 +75,7 @@ C----------------------------------------------------------------------- IF (KODE.EQ.1) GO TO 40 STR = ZR + ZETA2R STI = ZI + ZETA2I - RAST = FN/ZABS(STR,STI) + RAST = FN/AZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = -ZETA1R + STR @@ -95,7 +95,7 @@ C----------------------------------------------------------------------- C----------------------------------------------------------------------- C REFINE TEST AND SCALE C----------------------------------------------------------------------- - APHI = ZABS(PHIR,PHII) + APHI = AZABS(PHIR,PHII) RS1 = RS1 + DLOG(APHI) IF (DABS(RS1).GT.ELIM) GO TO 110 IF (I.EQ.1) IFLAG = 1 @@ -124,7 +124,7 @@ C----------------------------------------------------------------------- YI(M) = S2I*CSRR(IFLAG) 80 CONTINUE IF (ND.LE.2) GO TO 100 - RAST = 1.0D0/ZABS(ZR,ZI) + RAST = 1.0D0/AZABS(ZR,ZI) STR = ZR*RAST STI = -ZI*RAST RZR = (STR+STR)*RAST diff --git a/amos/zuni2.f b/amos/zuni2.f index 6f6b579..49061cb 100644 --- a/amos/zuni2.f +++ b/amos/zuni2.f @@ -13,7 +13,7 @@ C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL. C Y(I)=CZERO FOR I=NLAST+1,N C -C***ROUTINES CALLED ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,ZABS +C***ROUTINES CALLED ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,AZABS C***END PROLOGUE ZUNI2 C COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS, C *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN @@ -23,7 +23,7 @@ C *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN * DAIR, ELIM, FN, FNU, FNUL, HPI, PHII, PHIR, RAST, RAZ, RS1, RZI, * RZR, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZBI, ZBR, ZEROI, * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI, ZNI, ZNR, ZR, CYR, - * CYI, D1MACH, ZABS, CAR, SAR + * CYI, D1MACH, AZABS, CAR, SAR INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST, * NN, NUF, NW, NZ, IDUM DIMENSION BRY(3), YR(N), YI(N), CIPR(4), CIPI(4), CSSR(3), @@ -85,7 +85,7 @@ C----------------------------------------------------------------------- IF (KODE.EQ.1) GO TO 20 STR = ZBR + ZETA2R STI = ZBI + ZETA2I - RAST = FN/ZABS(STR,STI) + RAST = FN/AZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = -ZETA1R + STR @@ -106,7 +106,7 @@ C----------------------------------------------------------------------- IF (KODE.EQ.1) GO TO 50 STR = ZBR + ZETA2R STI = ZBI + ZETA2I - RAST = FN/ZABS(STR,STI) + RAST = FN/AZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = -ZETA1R + STR @@ -127,8 +127,8 @@ C----------------------------------------------------------------------- C REFINE TEST AND SCALE C----------------------------------------------------------------------- C----------------------------------------------------------------------- - APHI = ZABS(PHIR,PHII) - AARG = ZABS(ARGR,ARGI) + APHI = AZABS(PHIR,PHII) + AARG = AZABS(ARGR,ARGI) RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC IF (DABS(RS1).GT.ELIM) GO TO 120 IF (I.EQ.1) IFLAG = 1 @@ -171,7 +171,7 @@ C----------------------------------------------------------------------- C2R = STR 90 CONTINUE IF (ND.LE.2) GO TO 110 - RAZ = 1.0D0/ZABS(ZR,ZI) + RAZ = 1.0D0/AZABS(ZR,ZI) STR = ZR*RAZ STI = -ZI*RAZ RZR = (STR+STR)*RAZ diff --git a/amos/zunik.f b/amos/zunik.f index 3e8293e..7f297c3 100644 --- a/amos/zunik.f +++ b/amos/zunik.f @@ -18,7 +18,7 @@ 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 ZDIV,ZLOG,ZSQRT,D1MACH +C***ROUTINES CALLED ZDIV,AZLOG,AZSQRT,D1MACH C***END PROLOGUE ZUNIK C COMPLEX CFN,CON,CONE,CRFN,CWRK,CZERO,PHI,S,SR,SUM,T,T2,ZETA1, C *ZETA2,ZN,ZR @@ -131,11 +131,11 @@ C----------------------------------------------------------------------- TI = ZRI*RFN SR = CONER + (TR*TR-TI*TI) SI = CONEI + (TR*TI+TI*TR) - CALL ZSQRT(SR, SI, SRR, SRI) + CALL AZSQRT(SR, SI, SRR, SRI) STR = CONER + SRR STI = CONEI + SRI CALL ZDIV(STR, STI, TR, TI, ZNR, ZNI) - CALL ZLOG(ZNR, ZNI, STR, STI, IDUM) + CALL AZLOG(ZNR, ZNI, STR, STI, IDUM) ZETA1R = FNU*STR ZETA1I = FNU*STI ZETA2R = FNU*SRR @@ -143,7 +143,7 @@ C----------------------------------------------------------------------- CALL ZDIV(CONER, CONEI, SRR, SRI, TR, TI) SRR = TR*RFN SRI = TI*RFN - CALL ZSQRT(SRR, SRI, CWRKR(16), CWRKI(16)) + CALL AZSQRT(SRR, SRI, CWRKR(16), CWRKI(16)) PHIR = CWRKR(16)*CON(IKFLG) PHII = CWRKI(16)*CON(IKFLG) IF (IPMTR.NE.0) RETURN diff --git a/amos/zunk1.f b/amos/zunk1.f index ea9ce9e..5457d07 100644 --- a/amos/zunk1.f +++ b/amos/zunk1.f @@ -9,7 +9,7 @@ C UNIFORM ASYMPTOTIC EXPANSION. C MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION. C NZ=-1 MEANS AN OVERFLOW WILL OCCUR C -C***ROUTINES CALLED ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,ZABS +C***ROUTINES CALLED ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,AZABS C***END PROLOGUE ZUNK1 C COMPLEX CFN,CK,CONE,CRSC,CS,CSCL,CSGN,CSPN,CSR,CSS,CWRK,CY,CZERO, C *C1,C2,PHI,PHID,RZ,SUM,SUMD,S1,S2,Y,Z,ZETA1,ZETA1D,ZETA2,ZETA2D,ZR @@ -19,7 +19,7 @@ C *C1,C2,PHI,PHID,RZ,SUM,SUMD,S1,S2,Y,Z,ZETA1,ZETA1D,ZETA2,ZETA2D,ZR * FNF, FNU, PHIDI, PHIDR, PHII, PHIR, PI, RAST, RAZR, RS1, RZI, * RZR, SGN, STI, STR, SUMDI, SUMDR, SUMI, SUMR, S1I, S1R, S2I, * S2R, TOL, YI, YR, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, - * ZET1DI, ZET1DR, ZET2DI, ZET2DR, ZI, ZR, ZRI, ZRR, D1MACH, ZABS + * ZET1DI, ZET1DR, ZET2DI, ZET2DR, ZI, ZR, ZRI, ZRR, D1MACH, AZABS INTEGER I, IB, IFLAG, IFN, IL, INIT, INU, IUF, K, KDFLG, KFLAG, * KK, KODE, MR, N, NW, NZ, INITD, IC, IPARD, J DIMENSION BRY(3), INIT(2), YR(N), YI(N), SUMR(2), SUMI(2), @@ -65,7 +65,7 @@ C----------------------------------------------------------------------- IF (KODE.EQ.1) GO TO 20 STR = ZRR + ZETA2R(J) STI = ZRI + ZETA2I(J) - RAST = FN/ZABS(STR,STI) + RAST = FN/AZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = ZETA1R(J) - STR @@ -85,7 +85,7 @@ C----------------------------------------------------------------------- C----------------------------------------------------------------------- C REFINE TEST AND SCALE C----------------------------------------------------------------------- - APHI = ZABS(PHIR(J),PHII(J)) + APHI = AZABS(PHIR(J),PHII(J)) RS1 = RS1 + DLOG(APHI) IF (DABS(RS1).GT.ELIM) GO TO 60 IF (KDFLG.EQ.1) KFLAG = 1 @@ -133,7 +133,7 @@ C----------------------------------------------------------------------- 70 CONTINUE I = N 75 CONTINUE - RAZR = 1.0D0/ZABS(ZRR,ZRI) + RAZR = 1.0D0/AZABS(ZRR,ZRI) STR = ZRR*RAZR STI = -ZRI*RAZR RZR = (STR+STR)*RAZR @@ -156,7 +156,7 @@ C----------------------------------------------------------------------- IF (KODE.EQ.1) GO TO 80 STR = ZRR + ZET2DR STI = ZRI + ZET2DI - RAST = FN/ZABS(STR,STI) + RAST = FN/AZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = ZET1DR - STR @@ -172,7 +172,7 @@ C----------------------------------------------------------------------- C---------------------------------------------------------------------------- C REFINE ESTIMATE AND TEST C------------------------------------------------------------------------- - APHI = ZABS(PHIDR,PHIDI) + APHI = AZABS(PHIDR,PHIDI) RS1 = RS1+DLOG(APHI) IF (DABS(RS1).LT.ELIM) GO TO 100 95 CONTINUE @@ -287,7 +287,7 @@ C----------------------------------------------------------------------- IF (KODE.EQ.1) GO TO 200 STR = ZRR + ZET2DR STI = ZRI + ZET2DI - RAST = FN/ZABS(STR,STI) + RAST = FN/AZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = -ZET1DR + STR @@ -307,7 +307,7 @@ C----------------------------------------------------------------------- C----------------------------------------------------------------------- C REFINE TEST AND SCALE C----------------------------------------------------------------------- - APHI = ZABS(PHIDR,PHIDI) + APHI = AZABS(PHIDR,PHIDI) RS1 = RS1 + DLOG(APHI) IF (DABS(RS1).GT.ELIM) GO TO 260 IF (KDFLG.EQ.1) IFLAG = 1 diff --git a/amos/zunk2.f b/amos/zunk2.f index 6fa8f33..932fe8e 100644 --- a/amos/zunk2.f +++ b/amos/zunk2.f @@ -12,7 +12,7 @@ C HALF PLANE OR ZR=-Z IF Z IS IN THE LEFT HALF PLANE. MR INDIC- C ATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION. C NZ=-1 MEANS AN OVERFLOW WILL OCCUR C -C***ROUTINES CALLED ZAIRY,ZKSCL,ZS1S2,ZUCHK,ZUNHJ,D1MACH,ZABS +C***ROUTINES CALLED ZAIRY,ZKSCL,ZS1S2,ZUCHK,ZUNHJ,D1MACH,AZABS C***END PROLOGUE ZUNK2 C COMPLEX AI,ARG,ARGD,ASUM,ASUMD,BSUM,BSUMD,CFN,CI,CIP,CK,CONE,CRSC, C *CR1,CR2,CS,CSCL,CSGN,CSPN,CSR,CSS,CY,CZERO,C1,C2,DAI,PHI,PHID,RZ, @@ -26,7 +26,7 @@ C *S1,S2,Y,Z,ZB,ZETA1,ZETA1D,ZETA2,ZETA2D,ZN,ZR * PHII, PHIR, PI, PTI, PTR, RAST, RAZR, RS1, RZI, RZR, SAR, SGN, * STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, YY, ZBI, ZBR, ZEROI, * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZET1DI, ZET1DR, ZET2DI, - * ZET2DR, ZI, ZNI, ZNR, ZR, ZRI, ZRR, D1MACH, ZABS + * ZET2DR, ZI, ZNI, ZNR, ZR, ZRI, ZRR, D1MACH, AZABS INTEGER I, IB, IFLAG, IFN, IL, IN, INU, IUF, K, KDFLG, KFLAG, KK, * KODE, MR, N, NAI, NDAI, NW, NZ, IDUM, J, IPARD, IC DIMENSION BRY(3), YR(N), YI(N), ASUMR(2), ASUMI(2), BSUMR(2), @@ -105,7 +105,7 @@ C----------------------------------------------------------------------- IF (KODE.EQ.1) GO TO 30 STR = ZBR + ZETA2R(J) STI = ZBI + ZETA2I(J) - RAST = FN/ZABS(STR,STI) + RAST = FN/AZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = ZETA1R(J) - STR @@ -125,8 +125,8 @@ C----------------------------------------------------------------------- C----------------------------------------------------------------------- C REFINE TEST AND SCALE C----------------------------------------------------------------------- - APHI = ZABS(PHIR(J),PHII(J)) - AARG = ZABS(ARGR(J),ARGI(J)) + APHI = AZABS(PHIR(J),PHII(J)) + AARG = AZABS(ARGR(J),ARGI(J)) RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC IF (DABS(RS1).GT.ELIM) GO TO 70 IF (KDFLG.EQ.1) KFLAG = 1 @@ -193,7 +193,7 @@ C----------------------------------------------------------------------- 80 CONTINUE I = N 85 CONTINUE - RAZR = 1.0D0/ZABS(ZRR,ZRI) + RAZR = 1.0D0/AZABS(ZRR,ZRI) STR = ZRR*RAZR STI = -ZRI*RAZR RZR = (STR+STR)*RAZR @@ -214,7 +214,7 @@ C----------------------------------------------------------------------- IF (KODE.EQ.1) GO TO 90 STR = ZBR + ZET2DR STI = ZBI + ZET2DI - RAST = FN/ZABS(STR,STI) + RAST = FN/AZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = ZET1DR - STR @@ -230,7 +230,7 @@ C----------------------------------------------------------------------- C---------------------------------------------------------------------------- C REFINE ESTIMATE AND TEST C------------------------------------------------------------------------- - APHI = ZABS(PHIDR,PHIDI) + APHI = AZABS(PHIDR,PHIDI) RS1 = RS1+DLOG(APHI) IF (DABS(RS1).LT.ELIM) GO TO 120 105 CONTINUE @@ -291,7 +291,7 @@ C----------------------------------------------------------------------- FMR = DBLE(FLOAT(MR)) SGN = -DSIGN(PI,FMR) C----------------------------------------------------------------------- -C CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP. +C CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP. C----------------------------------------------------------------------- CSGNI = SGN IF (YY.LE.0.0D0) CSGNI = -CSGNI @@ -355,7 +355,7 @@ C----------------------------------------------------------------------- IF (KODE.EQ.1) GO TO 220 STR = ZBR + ZET2DR STI = ZBI + ZET2DI - RAST = FN/ZABS(STR,STI) + RAST = FN/AZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = -ZET1DR + STR @@ -375,8 +375,8 @@ C----------------------------------------------------------------------- C----------------------------------------------------------------------- C REFINE TEST AND SCALE C----------------------------------------------------------------------- - APHI = ZABS(PHIDR,PHIDI) - AARG = ZABS(ARGDR,ARGDI) + APHI = AZABS(PHIDR,PHIDI) + AARG = AZABS(ARGDR,ARGDI) RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC IF (DABS(RS1).GT.ELIM) GO TO 280 IF (KDFLG.EQ.1) IFLAG = 1 diff --git a/amos/zuoik.f b/amos/zuoik.f index 4c1259c..5b05f96 100644 --- a/amos/zuoik.f +++ b/amos/zuoik.f @@ -23,7 +23,7 @@ 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 ZUCHK,ZUNHJ,ZUNIK,D1MACH,ZABS,ZLOG +C***ROUTINES CALLED ZUCHK,ZUNHJ,ZUNIK,D1MACH,AZABS,AZLOG C***END PROLOGUE ZUOIK C COMPLEX ARG,ASUM,BSUM,CWRK,CZ,CZERO,PHI,SUM,Y,Z,ZB,ZETA1,ZETA2,ZN, C *ZR @@ -31,7 +31,7 @@ C *ZR * ASCLE, AX, AY, BSUMI, BSUMR, CWRKI, CWRKR, CZI, CZR, ELIM, FNN, * FNU, GNN, GNU, PHII, PHIR, RCZ, STR, STI, SUMI, SUMR, TOL, YI, * YR, ZBI, ZBR, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI, - * ZNI, ZNR, ZR, ZRI, ZRR, D1MACH, ZABS + * ZNI, ZNR, ZR, ZRI, ZRR, D1MACH, AZABS INTEGER I, IDUM, IFORM, IKFLG, INIT, KODE, N, NN, NUF, NW DIMENSION YR(N), YI(N), CWRKR(16), CWRKI(16) DATA ZEROR,ZEROI / 0.0D0, 0.0D0 / @@ -78,7 +78,7 @@ C----------------------------------------------------------------------- * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) CZR = -ZETA1R + ZETA2R CZI = -ZETA1I + ZETA2I - AARG = ZABS(ARGR,ARGI) + AARG = AZABS(ARGR,ARGI) 50 CONTINUE IF (KODE.EQ.1) GO TO 60 CZR = CZR - ZBR @@ -88,7 +88,7 @@ C----------------------------------------------------------------------- CZR = -CZR CZI = -CZI 70 CONTINUE - APHI = ZABS(PHIR,PHII) + APHI = AZABS(PHIR,PHII) RCZ = CZR C----------------------------------------------------------------------- C OVERFLOW TEST @@ -117,11 +117,11 @@ C----------------------------------------------------------------------- RETURN 110 CONTINUE ASCLE = 1.0D+3*D1MACH(1)/TOL - CALL ZLOG(PHIR, PHII, STR, STI, IDUM) + CALL AZLOG(PHIR, PHII, STR, STI, IDUM) CZR = CZR + STR CZI = CZI + STI IF (IFORM.EQ.1) GO TO 120 - CALL ZLOG(ARGR, ARGI, STR, STI, IDUM) + CALL AZLOG(ARGR, ARGI, STR, STI, IDUM) CZR = CZR - 0.25D0*STR - AIC CZI = CZI - 0.25D0*STI 120 CONTINUE @@ -151,13 +151,13 @@ C----------------------------------------------------------------------- * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) CZR = -ZETA1R + ZETA2R CZI = -ZETA1I + ZETA2I - AARG = ZABS(ARGR,ARGI) + AARG = AZABS(ARGR,ARGI) 160 CONTINUE IF (KODE.EQ.1) GO TO 170 CZR = CZR - ZBR CZI = CZI - ZBI 170 CONTINUE - APHI = ZABS(PHIR,PHII) + APHI = AZABS(PHIR,PHII) RCZ = CZR IF (RCZ.LT.(-ELIM)) GO TO 180 IF (RCZ.GT.(-ALIM)) RETURN @@ -173,11 +173,11 @@ C----------------------------------------------------------------------- GO TO 140 190 CONTINUE ASCLE = 1.0D+3*D1MACH(1)/TOL - CALL ZLOG(PHIR, PHII, STR, STI, IDUM) + CALL AZLOG(PHIR, PHII, STR, STI, IDUM) CZR = CZR + STR CZI = CZI + STI IF (IFORM.EQ.1) GO TO 200 - CALL ZLOG(ARGR, ARGI, STR, STI, IDUM) + CALL AZLOG(ARGR, ARGI, STR, STI, IDUM) CZR = CZR - 0.25D0*STR - AIC CZI = CZI - 0.25D0*STI 200 CONTINUE diff --git a/amos/zwrsk.f b/amos/zwrsk.f index 5f0672f..397340f 100644 --- a/amos/zwrsk.f +++ b/amos/zwrsk.f @@ -6,12 +6,12 @@ C C ZWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY C NORMALIZING THE I FUNCTION RATIOS FROM ZRATI BY THE WRONSKIAN C -C***ROUTINES CALLED D1MACH,ZBKNU,ZRATI,ZABS +C***ROUTINES CALLED D1MACH,ZBKNU,ZRATI,AZABS C***END PROLOGUE ZWRSK C COMPLEX CINU,CSCL,CT,CW,C1,C2,RCT,ST,Y,ZR DOUBLE PRECISION ACT, ACW, ALIM, ASCLE, CINUI, CINUR, CSCLR, CTI, * CTR, CWI, CWR, C1I, C1R, C2I, C2R, ELIM, FNU, PTI, PTR, RACT, - * STI, STR, TOL, YI, YR, ZRI, ZRR, ZABS, D1MACH + * STI, STR, TOL, YI, YR, ZRI, ZRR, AZABS, D1MACH INTEGER I, KODE, N, NW, NZ DIMENSION YR(N), YI(N), CWR(2), CWI(2) C----------------------------------------------------------------------- @@ -39,7 +39,7 @@ 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 = ZABS(CWR(2),CWI(2)) + ACW = AZABS(CWR(2),CWI(2)) ASCLE = 1.0D+3*D1MACH(1)/TOL CSCLR = 1.0D0 IF (ACW.GT.ASCLE) GO TO 20 @@ -66,7 +66,7 @@ C----------------------------------------------------------------------- PTI = PTI + C2I CTR = ZRR*PTR - ZRI*PTI CTI = ZRR*PTI + ZRI*PTR - ACT = ZABS(CTR,CTI) + ACT = AZABS(CTR,CTI) RACT = 1.0D0/ACT CTR = CTR*RACT CTI = -CTI*RACT