amos/zmlri.f

   1       SUBROUTINE ZMLRI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL)
   2 C***BEGIN PROLOGUE  ZMLRI
   3 C***REFER TO  ZBESI,ZBESK
   4 C
   5 C     ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE
   6 C     MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES.
   7 C
   8 C***ROUTINES CALLED  DGAMLN,D1MACH,AZABS,AZEXP,AZLOG,ZMLT
   9 C***END PROLOGUE  ZMLRI
  10 C     COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z
  11       DOUBLE PRECISION ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI,
  12      * CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, FNU, PTI, PTR, P1I,
  13      * P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI,
  14      * SUMR, TFNF, TOL, TST, YI, YR, ZEROI, ZEROR, ZI, ZR, DGAMLN,
  15      * D1MACH, AZABS
  16       INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N, NZ
  17       DIMENSION YR(N), YI(N)
  18       DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
  19       SCLE = D1MACH(1)/TOL
  20       NZ=0
  21       AZ = AZABS(ZR,ZI)
  22       IAZ = INT(SNGL(AZ))
  23       IFNU = INT(SNGL(FNU))
  24       INU = IFNU + N - 1
  25       AT = DBLE(FLOAT(IAZ)) + 1.0D0
  26       RAZ = 1.0D0/AZ
  27       STR = ZR*RAZ
  28       STI = -ZI*RAZ
  29       CKR = STR*AT*RAZ
  30       CKI = STI*AT*RAZ
  31       RZR = (STR+STR)*RAZ
  32       RZI = (STI+STI)*RAZ
  33       P1R = ZEROR
  34       P1I = ZEROI
  35       P2R = CONER
  36       P2I = CONEI
  37       ACK = (AT+1.0D0)*RAZ
  38       RHO = ACK + DSQRT(ACK*ACK-1.0D0)
  39       RHO2 = RHO*RHO
  40       TST = (RHO2+RHO2)/((RHO2-1.0D0)*(RHO-1.0D0))
  41       TST = TST/TOL
  42 C-----------------------------------------------------------------------
  43 C     COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
  44 C-----------------------------------------------------------------------
  45       AK = AT
  46       DO 10 I=1,80
  47         PTR = P2R
  48         PTI = P2I
  49         P2R = P1R - (CKR*PTR-CKI*PTI)
  50         P2I = P1I - (CKI*PTR+CKR*PTI)
  51         P1R = PTR
  52         P1I = PTI
  53         CKR = CKR + RZR
  54         CKI = CKI + RZI
  55         AP = AZABS(P2R,P2I)
  56         IF (AP.GT.TST*AK*AK) GO TO 20
  57         AK = AK + 1.0D0
  58    10 CONTINUE
  59       GO TO 110
  60    20 CONTINUE
  61       I = I + 1
  62       K = 0
  63       IF (INU.LT.IAZ) GO TO 40
  64 C-----------------------------------------------------------------------
  65 C     COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
  66 C-----------------------------------------------------------------------
  67       P1R = ZEROR
  68       P1I = ZEROI
  69       P2R = CONER
  70       P2I = CONEI
  71       AT = DBLE(FLOAT(INU)) + 1.0D0
  72       STR = ZR*RAZ
  73       STI = -ZI*RAZ
  74       CKR = STR*AT*RAZ
  75       CKI = STI*AT*RAZ
  76       ACK = AT*RAZ
  77       TST = DSQRT(ACK/TOL)
  78       ITIME = 1
  79       DO 30 K=1,80
  80         PTR = P2R
  81         PTI = P2I
  82         P2R = P1R - (CKR*PTR-CKI*PTI)
  83         P2I = P1I - (CKR*PTI+CKI*PTR)
  84         P1R = PTR
  85         P1I = PTI
  86         CKR = CKR + RZR
  87         CKI = CKI + RZI
  88         AP = AZABS(P2R,P2I)
  89         IF (AP.LT.TST) GO TO 30
  90         IF (ITIME.EQ.2) GO TO 40
  91         ACK = AZABS(CKR,CKI)
  92         FLAM = ACK + DSQRT(ACK*ACK-1.0D0)
  93         FKAP = AP/AZABS(P1R,P1I)
  94         RHO = DMIN1(FLAM,FKAP)
  95         TST = TST*DSQRT(RHO/(RHO*RHO-1.0D0))
  96         ITIME = 2
  97    30 CONTINUE
  98       GO TO 110
  99    40 CONTINUE
 100 C-----------------------------------------------------------------------
 101 C     BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
 102 C-----------------------------------------------------------------------
 103       K = K + 1
 104       KK = MAX0(I+IAZ,K+INU)
 105       FKK = DBLE(FLOAT(KK))
 106       P1R = ZEROR
 107       P1I = ZEROI
 108 C-----------------------------------------------------------------------
 109 C     SCALE P2 AND SUM BY SCLE
 110 C-----------------------------------------------------------------------
 111       P2R = SCLE
 112       P2I = ZEROI
 113       FNF = FNU - DBLE(FLOAT(IFNU))
 114       TFNF = FNF + FNF
 115       BK = DGAMLN(FKK+TFNF+1.0D0,IDUM) - DGAMLN(FKK+1.0D0,IDUM) -
 116      * DGAMLN(TFNF+1.0D0,IDUM)
 117       BK = DEXP(BK)
 118       SUMR = ZEROR
 119       SUMI = ZEROI
 120       KM = KK - INU
 121       DO 50 I=1,KM
 122         PTR = P2R
 123         PTI = P2I
 124         P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
 125         P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI)
 126         P1R = PTR
 127         P1I = PTI
 128         AK = 1.0D0 - TFNF/(FKK+TFNF)
 129         ACK = BK*AK
 130         SUMR = SUMR + (ACK+BK)*P1R
 131         SUMI = SUMI + (ACK+BK)*P1I
 132         BK = ACK
 133         FKK = FKK - 1.0D0
 134    50 CONTINUE
 135       YR(N) = P2R
 136       YI(N) = P2I
 137       IF (N.EQ.1) GO TO 70
 138       DO 60 I=2,N
 139         PTR = P2R
 140         PTI = P2I
 141         P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
 142         P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI)
 143         P1R = PTR
 144         P1I = PTI
 145         AK = 1.0D0 - TFNF/(FKK+TFNF)
 146         ACK = BK*AK
 147         SUMR = SUMR + (ACK+BK)*P1R
 148         SUMI = SUMI + (ACK+BK)*P1I
 149         BK = ACK
 150         FKK = FKK - 1.0D0
 151         M = N - I + 1
 152         YR(M) = P2R
 153         YI(M) = P2I
 154    60 CONTINUE
 155    70 CONTINUE
 156       IF (IFNU.LE.0) GO TO 90
 157       DO 80 I=1,IFNU
 158         PTR = P2R
 159         PTI = P2I
 160         P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
 161         P2I = P1I + (FKK+FNF)*(RZR*PTI+RZI*PTR)
 162         P1R = PTR
 163         P1I = PTI
 164         AK = 1.0D0 - TFNF/(FKK+TFNF)
 165         ACK = BK*AK
 166         SUMR = SUMR + (ACK+BK)*P1R
 167         SUMI = SUMI + (ACK+BK)*P1I
 168         BK = ACK
 169         FKK = FKK - 1.0D0
 170    80 CONTINUE
 171    90 CONTINUE
 172       PTR = ZR
 173       PTI = ZI
 174       IF (KODE.EQ.2) PTR = ZEROR
 175       CALL AZLOG(RZR, RZI, STR, STI, IDUM)
 176       P1R = -FNF*STR + PTR
 177       P1I = -FNF*STI + PTI
 178       AP = DGAMLN(1.0D0+FNF,IDUM)
 179       PTR = P1R - AP
 180       PTI = P1I
 181 C-----------------------------------------------------------------------
 182 C     THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
 183 C     IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
 184 C-----------------------------------------------------------------------
 185       P2R = P2R + SUMR
 186       P2I = P2I + SUMI
 187       AP = AZABS(P2R,P2I)
 188       P1R = 1.0D0/AP
 189       CALL AZEXP(PTR, PTI, STR, STI)
 190       CKR = STR*P1R
 191       CKI = STI*P1R
 192       PTR = P2R*P1R
 193       PTI = -P2I*P1R
 194       CALL ZMLT(CKR, CKI, PTR, PTI, CNORMR, CNORMI)
 195       DO 100 I=1,N
 196         STR = YR(I)*CNORMR - YI(I)*CNORMI
 197         YI(I) = YR(I)*CNORMI + YI(I)*CNORMR
 198         YR(I) = STR
 199   100 CONTINUE
 200       RETURN
 201   110 CONTINUE
 202       NZ=-2
 203       RETURN
 204       END