133 lines
3.8 KiB
Fortran
133 lines
3.8 KiB
Fortran
SUBROUTINE ZRATI(ZR, ZI, FNU, N, CYR, CYI, TOL)
|
|
C***BEGIN PROLOGUE ZRATI
|
|
C***REFER TO ZBESI,ZBESK,ZBESH
|
|
C
|
|
C ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD
|
|
C RECURRENCE. THE STARTING INDEX IS DETERMINED BY FORWARD
|
|
C RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B,
|
|
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 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, 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 = AZABS(ZR,ZI)
|
|
INU = INT(SNGL(FNU))
|
|
IDNU = INU + N - 1
|
|
MAGZ = INT(SNGL(AZ))
|
|
AMAGZ = DBLE(FLOAT(MAGZ+1))
|
|
FDNU = DBLE(FLOAT(IDNU))
|
|
FNUP = DMAX1(AMAGZ,FDNU)
|
|
ID = IDNU - MAGZ - 1
|
|
ITIME = 1
|
|
K = 1
|
|
PTR = 1.0D0/AZ
|
|
RZR = PTR*(ZR+ZR)*PTR
|
|
RZI = -PTR*(ZI+ZI)*PTR
|
|
T1R = RZR*FNUP
|
|
T1I = RZI*FNUP
|
|
P2R = -T1R
|
|
P2I = -T1I
|
|
P1R = CONER
|
|
P1I = CONEI
|
|
T1R = T1R + RZR
|
|
T1I = T1I + RZI
|
|
IF (ID.GT.0) ID = 0
|
|
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
|
|
C P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR
|
|
C PREMATURELY.
|
|
C-----------------------------------------------------------------------
|
|
ARG = (AP2+AP2)/(AP1*TOL)
|
|
TEST1 = DSQRT(ARG)
|
|
TEST = TEST1
|
|
RAP1 = 1.0D0/AP1
|
|
P1R = P1R*RAP1
|
|
P1I = P1I*RAP1
|
|
P2R = P2R*RAP1
|
|
P2I = P2I*RAP1
|
|
AP2 = AP2*RAP1
|
|
10 CONTINUE
|
|
K = K + 1
|
|
AP1 = AP2
|
|
PTR = P2R
|
|
PTI = P2I
|
|
P2R = P1R - (T1R*PTR-T1I*PTI)
|
|
P2I = P1I - (T1R*PTI+T1I*PTR)
|
|
P1R = PTR
|
|
P1I = PTI
|
|
T1R = T1R + RZR
|
|
T1I = T1I + RZI
|
|
AP2 = AZABS(P2R,P2I)
|
|
IF (AP1.LE.TEST) GO TO 10
|
|
IF (ITIME.EQ.2) GO TO 20
|
|
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))
|
|
ITIME = 2
|
|
GO TO 10
|
|
20 CONTINUE
|
|
KK = K + 1 - ID
|
|
AK = DBLE(FLOAT(KK))
|
|
T1R = AK
|
|
T1I = CZEROI
|
|
DFNU = FNU + DBLE(FLOAT(N-1))
|
|
P1R = 1.0D0/AP2
|
|
P1I = CZEROI
|
|
P2R = CZEROR
|
|
P2I = CZEROI
|
|
DO 30 I=1,KK
|
|
PTR = P1R
|
|
PTI = P1I
|
|
RAP1 = DFNU + T1R
|
|
TTR = RZR*RAP1
|
|
TTI = RZI*RAP1
|
|
P1R = (PTR*TTR-PTI*TTI) + P2R
|
|
P1I = (PTR*TTI+PTI*TTR) + P2I
|
|
P2R = PTR
|
|
P2I = PTI
|
|
T1R = T1R - CONER
|
|
30 CONTINUE
|
|
IF (P1R.NE.CZEROR .OR. P1I.NE.CZEROI) GO TO 40
|
|
P1R = TOL
|
|
P1I = TOL
|
|
40 CONTINUE
|
|
CALL ZDIV(P2R, P2I, P1R, P1I, CYR(N), CYI(N))
|
|
IF (N.EQ.1) RETURN
|
|
K = N - 1
|
|
AK = DBLE(FLOAT(K))
|
|
T1R = AK
|
|
T1I = CZEROI
|
|
CDFNUR = FNU*RZR
|
|
CDFNUI = FNU*RZI
|
|
DO 60 I=2,N
|
|
PTR = CDFNUR + (T1R*RZR-T1I*RZI) + CYR(K+1)
|
|
PTI = CDFNUI + (T1R*RZI+T1I*RZR) + CYI(K+1)
|
|
AK = AZABS(PTR,PTI)
|
|
IF (AK.NE.CZEROR) GO TO 50
|
|
PTR = TOL
|
|
PTI = TOL
|
|
AK = TOL*RT2
|
|
50 CONTINUE
|
|
RAK = CONER/AK
|
|
CYR(K) = RAK*PTR*RAK
|
|
CYI(K) = -RAK*PTI*RAK
|
|
T1R = T1R - CONER
|
|
K = K - 1
|
|
60 CONTINUE
|
|
RETURN
|
|
END
|