amos/zrati.f

   1       SUBROUTINE ZRATI(ZR, ZI, FNU, N, CYR, CYI, TOL)
   2 C***BEGIN PROLOGUE  ZRATI
   3 C***REFER TO  ZBESI,ZBESK,ZBESH
   4 C
   5 C     ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD
   6 C     RECURRENCE.  THE STARTING INDEX IS DETERMINED BY FORWARD
   7 C     RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B,
   8 C     MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973,
   9 C     BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER,
  10 C     BY D. J. SOOKNE.
  11 C
  12 C***ROUTINES CALLED  AZABS,ZDIV
  13 C***END PROLOGUE  ZRATI
  14 C     COMPLEX Z,CY(1),CONE,CZERO,P1,P2,T1,RZ,PT,CDFNU
  15       DOUBLE PRECISION AK, AMAGZ, AP1, AP2, ARG, AZ, CDFNUI, CDFNUR,
  16      * CONEI, CONER, CYI, CYR, CZEROI, CZEROR, DFNU, FDNU, FLAM, FNU,
  17      * FNUP, PTI, PTR, P1I, P1R, P2I, P2R, RAK, RAP1, RHO, RT2, RZI,
  18      * RZR, TEST, TEST1, TOL, TTI, TTR, T1I, T1R, ZI, ZR, AZABS
  19       INTEGER I, ID, IDNU, INU, ITIME, K, KK, MAGZ, N
  20       DIMENSION CYR(N), CYI(N)
  21       DATA CZEROR,CZEROI,CONER,CONEI,RT2/
  22      1 0.0D0, 0.0D0, 1.0D0, 0.0D0, 1.41421356237309505D0 /
  23       AZ = AZABS(ZR,ZI)
  24       INU = INT(SNGL(FNU))
  25       IDNU = INU + N - 1
  26       MAGZ = INT(SNGL(AZ))
  27       AMAGZ = DBLE(FLOAT(MAGZ+1))
  28       FDNU = DBLE(FLOAT(IDNU))
  29       FNUP = DMAX1(AMAGZ,FDNU)
  30       ID = IDNU - MAGZ - 1
  31       ITIME = 1
  32       K = 1
  33       PTR = 1.0D0/AZ
  34       RZR = PTR*(ZR+ZR)*PTR
  35       RZI = -PTR*(ZI+ZI)*PTR
  36       T1R = RZR*FNUP
  37       T1I = RZI*FNUP
  38       P2R = -T1R
  39       P2I = -T1I
  40       P1R = CONER
  41       P1I = CONEI
  42       T1R = T1R + RZR
  43       T1I = T1I + RZI
  44       IF (ID.GT.0) ID = 0
  45       AP2 = AZABS(P2R,P2I)
  46       AP1 = AZABS(P1R,P1I)
  47 C-----------------------------------------------------------------------
  48 C     THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU
  49 C     GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT
  50 C     P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR
  51 C     PREMATURELY.
  52 C-----------------------------------------------------------------------
  53       ARG = (AP2+AP2)/(AP1*TOL)
  54       TEST1 = DSQRT(ARG)
  55       TEST = TEST1
  56       RAP1 = 1.0D0/AP1
  57       P1R = P1R*RAP1
  58       P1I = P1I*RAP1
  59       P2R = P2R*RAP1
  60       P2I = P2I*RAP1
  61       AP2 = AP2*RAP1
  62    10 CONTINUE
  63       K = K + 1
  64       AP1 = AP2
  65       PTR = P2R
  66       PTI = P2I
  67       P2R = P1R - (T1R*PTR-T1I*PTI)
  68       P2I = P1I - (T1R*PTI+T1I*PTR)
  69       P1R = PTR
  70       P1I = PTI
  71       T1R = T1R + RZR
  72       T1I = T1I + RZI
  73       AP2 = AZABS(P2R,P2I)
  74       IF (AP1.LE.TEST) GO TO 10
  75       IF (ITIME.EQ.2) GO TO 20
  76       AK = AZABS(T1R,T1I)*0.5D0
  77       FLAM = AK + DSQRT(AK*AK-1.0D0)
  78       RHO = DMIN1(AP2/AP1,FLAM)
  79       TEST = TEST1*DSQRT(RHO/(RHO*RHO-1.0D0))
  80       ITIME = 2
  81       GO TO 10
  82    20 CONTINUE
  83       KK = K + 1 - ID
  84       AK = DBLE(FLOAT(KK))
  85       T1R = AK
  86       T1I = CZEROI
  87       DFNU = FNU + DBLE(FLOAT(N-1))
  88       P1R = 1.0D0/AP2
  89       P1I = CZEROI
  90       P2R = CZEROR
  91       P2I = CZEROI
  92       DO 30 I=1,KK
  93         PTR = P1R
  94         PTI = P1I
  95         RAP1 = DFNU + T1R
  96         TTR = RZR*RAP1
  97         TTI = RZI*RAP1
  98         P1R = (PTR*TTR-PTI*TTI) + P2R
  99         P1I = (PTR*TTI+PTI*TTR) + P2I
 100         P2R = PTR
 101         P2I = PTI
 102         T1R = T1R - CONER
 103    30 CONTINUE
 104       IF (P1R.NE.CZEROR .OR. P1I.NE.CZEROI) GO TO 40
 105       P1R = TOL
 106       P1I = TOL
 107    40 CONTINUE
 108       CALL ZDIV(P2R, P2I, P1R, P1I, CYR(N), CYI(N))
 109       IF (N.EQ.1) RETURN
 110       K = N - 1
 111       AK = DBLE(FLOAT(K))
 112       T1R = AK
 113       T1I = CZEROI
 114       CDFNUR = FNU*RZR
 115       CDFNUI = FNU*RZI
 116       DO 60 I=2,N
 117         PTR = CDFNUR + (T1R*RZR-T1I*RZI) + CYR(K+1)
 118         PTI = CDFNUI + (T1R*RZI+T1I*RZR) + CYI(K+1)
 119         AK = AZABS(PTR,PTI)
 120         IF (AK.NE.CZEROR) GO TO 50
 121         PTR = TOL
 122         PTI = TOL
 123         AK = TOL*RT2
 124    50   CONTINUE
 125         RAK = CONER/AK
 126         CYR(K) = RAK*PTR*RAK
 127         CYI(K) = -RAK*PTI*RAK
 128         T1R = T1R - CONER
 129         K = K - 1
 130    60 CONTINUE
 131       RETURN
 132       END