QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

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
This commit is contained in:
Marek Nečada 2019-03-20 12:25:42 +02:00
parent e6e7485ebb
commit 522a6d9b85
31 changed files with 1117 additions and 172 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
amos/README.md Normal file
View File

@ -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.
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
```

14
amos/zabs.f
View File

@ -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

6
amos/zacai.f
View File

@ -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

8
amos/zacon.f
View File

@ -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

18
amos/zairy.f
View File

@ -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-----------------------------------------------------------------------

14
amos/zasyi.f
View File

@ -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

6
amos/zbesh.f
View File

@ -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)

269
amos/zbesi.f Normal file
View File

@ -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

4
amos/zbesj.f
View File

@ -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

281
amos/zbesk.f Normal file
View File

@ -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

6
amos/zbinu.f
View File

@ -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

364
amos/zbiry.f Normal file
View File

@ -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

26
amos/zbknu.f
View File

@ -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

8
amos/zbuni.f
View File

@ -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

6
amos/zdiv.f
View File

@ -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

6
amos/zexp.f
View File

@ -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)

12
amos/zkscl.f
View File

@ -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

12
amos/zlog.f
View File

@ -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

20
amos/zmlri.f
View File

@ -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

16
amos/zrati.f
View File

@ -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

14
amos/zs1s2.f
View File

@ -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)

12
amos/zseri.f
View File

@ -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

12
amos/zsqrt.f
View File

@ -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

20
amos/zunhj.f
View File

@ -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

12
amos/zuni1.f
View File

@ -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

14
amos/zuni2.f
View File

@ -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

8
amos/zunik.f
View File

@ -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

18
amos/zunk1.f
View File

@ -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

24
amos/zunk2.f
View File

@ -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

20
amos/zuoik.f
View File

@ -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

8
amos/zwrsk.f
View File

@ -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