src/numerical/slatec/fortran/zrati.f

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