cmlri.f

Go to the documentation of this file.
00001       SUBROUTINE CMLRI(Z, FNU, KODE, N, Y, NZ, TOL)
00002 C***BEGIN PROLOGUE  CMLRI
00003 C***REFER TO  CBESI,CBESK
00004 C
00005 C     CMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE
00006 C     MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES.
00007 C
00008 C***ROUTINES CALLED  GAMLN,R1MACH
00009 C***END PROLOGUE  CMLRI
00010       COMPLEX CK, CNORM, CONE, CTWO, CZERO, PT, P1, P2, RZ, SUM, Y, Z
00011       REAL ACK, AK, AP, AT, AZ, BK, FKAP, FKK, FLAM, FNF, FNU, RHO,
00012      * RHO2, SCLE, TFNF, TOL, TST, X, GAMLN, R1MACH
00013       INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N
00014       DIMENSION Y(N)
00015       DATA CZERO,CONE,CTWO /(0.0E0,0.0E0),(1.0E0,0.0E0),(2.0E0,0.0E0)/
00016       SCLE = 1.0E+3*R1MACH(1)/TOL
00017       NZ=0
00018       AZ = CABS(Z)
00019       X = REAL(Z)
00020       IAZ = INT(AZ)
00021       IFNU = INT(FNU)
00022       INU = IFNU + N - 1
00023       AT = FLOAT(IAZ) + 1.0E0
00024       CK = CMPLX(AT,0.0E0)/Z
00025       RZ = CTWO/Z
00026       P1 = CZERO
00027       P2 = CONE
00028       ACK = (AT+1.0E0)/AZ
00029       RHO = ACK + SQRT(ACK*ACK-1.0E0)
00030       RHO2 = RHO*RHO
00031       TST = (RHO2+RHO2)/((RHO2-1.0E0)*(RHO-1.0E0))
00032       TST = TST/TOL
00033 C-----------------------------------------------------------------------
00034 C     COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
00035 C-----------------------------------------------------------------------
00036       AK = AT
00037       DO 10 I=1,80
00038         PT = P2
00039         P2 = P1 - CK*P2
00040         P1 = PT
00041         CK = CK + RZ
00042         AP = CABS(P2)
00043         IF (AP.GT.TST*AK*AK) GO TO 20
00044         AK = AK + 1.0E0
00045    10 CONTINUE
00046       GO TO 110
00047    20 CONTINUE
00048       I = I + 1
00049       K = 0
00050       IF (INU.LT.IAZ) GO TO 40
00051 C-----------------------------------------------------------------------
00052 C     COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
00053 C-----------------------------------------------------------------------
00054       P1 = CZERO
00055       P2 = CONE
00056       AT = FLOAT(INU) + 1.0E0
00057       CK = CMPLX(AT,0.0E0)/Z
00058       ACK = AT/AZ
00059       TST = SQRT(ACK/TOL)
00060       ITIME = 1
00061       DO 30 K=1,80
00062         PT = P2
00063         P2 = P1 - CK*P2
00064         P1 = PT
00065         CK = CK + RZ
00066         AP = CABS(P2)
00067         IF (AP.LT.TST) GO TO 30
00068         IF (ITIME.EQ.2) GO TO 40
00069         ACK = CABS(CK)
00070         FLAM = ACK + SQRT(ACK*ACK-1.0E0)
00071         FKAP = AP/CABS(P1)
00072         RHO = AMIN1(FLAM,FKAP)
00073         TST = TST*SQRT(RHO/(RHO*RHO-1.0E0))
00074         ITIME = 2
00075    30 CONTINUE
00076       GO TO 110
00077    40 CONTINUE
00078 C-----------------------------------------------------------------------
00079 C     BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
00080 C-----------------------------------------------------------------------
00081       K = K + 1
00082       KK = MAX0(I+IAZ,K+INU)
00083       FKK = FLOAT(KK)
00084       P1 = CZERO
00085 C-----------------------------------------------------------------------
00086 C     SCALE P2 AND SUM BY SCLE
00087 C-----------------------------------------------------------------------
00088       P2 = CMPLX(SCLE,0.0E0)
00089       FNF = FNU - FLOAT(IFNU)
00090       TFNF = FNF + FNF
00091       BK = GAMLN(FKK+TFNF+1.0E0,IDUM) - GAMLN(FKK+1.0E0,IDUM)
00092      *     -GAMLN(TFNF+1.0E0,IDUM)
00093       BK = EXP(BK)
00094       SUM = CZERO
00095       KM = KK - INU
00096       DO 50 I=1,KM
00097         PT = P2
00098         P2 = P1 + CMPLX(FKK+FNF,0.0E0)*RZ*P2
00099         P1 = PT
00100         AK = 1.0E0 - TFNF/(FKK+TFNF)
00101         ACK = BK*AK
00102         SUM = SUM + CMPLX(ACK+BK,0.0E0)*P1
00103         BK = ACK
00104         FKK = FKK - 1.0E0
00105    50 CONTINUE
00106       Y(N) = P2
00107       IF (N.EQ.1) GO TO 70
00108       DO 60 I=2,N
00109         PT = P2
00110         P2 = P1 + CMPLX(FKK+FNF,0.0E0)*RZ*P2
00111         P1 = PT
00112         AK = 1.0E0 - TFNF/(FKK+TFNF)
00113         ACK = BK*AK
00114         SUM = SUM + CMPLX(ACK+BK,0.0E0)*P1
00115         BK = ACK
00116         FKK = FKK - 1.0E0
00117         M = N - I + 1
00118         Y(M) = P2
00119    60 CONTINUE
00120    70 CONTINUE
00121       IF (IFNU.LE.0) GO TO 90
00122       DO 80 I=1,IFNU
00123         PT = P2
00124         P2 = P1 + CMPLX(FKK+FNF,0.0E0)*RZ*P2
00125         P1 = PT
00126         AK = 1.0E0 - TFNF/(FKK+TFNF)
00127         ACK = BK*AK
00128         SUM = SUM + CMPLX(ACK+BK,0.0E0)*P1
00129         BK = ACK
00130         FKK = FKK - 1.0E0
00131    80 CONTINUE
00132    90 CONTINUE
00133       PT = Z
00134       IF (KODE.EQ.2) PT = PT - CMPLX(X,0.0E0)
00135       P1 = -CMPLX(FNF,0.0E0)*CLOG(RZ) + PT
00136       AP = GAMLN(1.0E0+FNF,IDUM)
00137       PT = P1 - CMPLX(AP,0.0E0)
00138 C-----------------------------------------------------------------------
00139 C     THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
00140 C     IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
00141 C-----------------------------------------------------------------------
00142       P2 = P2 + SUM
00143       AP = CABS(P2)
00144       P1 = CMPLX(1.0E0/AP,0.0E0)
00145       CK = CEXP(PT)*P1
00146       PT = CONJG(P2)*P1
00147       CNORM = CK*PT
00148       DO 100 I=1,N
00149         Y(I) = Y(I)*CNORM
00150   100 CONTINUE
00151       RETURN
00152   110 CONTINUE
00153       NZ=-2
00154       RETURN
00155       END
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines