zmlri.f

Go to the documentation of this file.
00001       SUBROUTINE ZMLRI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL)
00002 C***BEGIN PROLOGUE  ZMLRI
00003 C***REFER TO  ZBESI,ZBESK
00004 C
00005 C     ZMLRI 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  DGAMLN,D1MACH,XZABS,XZEXP,XZLOG,ZMLT
00009 C***END PROLOGUE  ZMLRI
00010 C     COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z
00011       DOUBLE PRECISION ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI,
00012      * CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, FNU, PTI, PTR, P1I,
00013      * P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI,
00014      * SUMR, TFNF, TOL, TST, YI, YR, ZEROI, ZEROR, ZI, ZR, DGAMLN,
00015      * D1MACH, XZABS
00016       INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N, NZ
00017       DIMENSION YR(N), YI(N)
00018       DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
00019       SCLE = D1MACH(1)/TOL
00020       NZ=0
00021       AZ = XZABS(ZR,ZI)
00022       IAZ = INT(SNGL(AZ))
00023       IFNU = INT(SNGL(FNU))
00024       INU = IFNU + N - 1
00025       AT = DBLE(FLOAT(IAZ)) + 1.0D0
00026       RAZ = 1.0D0/AZ
00027       STR = ZR*RAZ
00028       STI = -ZI*RAZ
00029       CKR = STR*AT*RAZ
00030       CKI = STI*AT*RAZ
00031       RZR = (STR+STR)*RAZ
00032       RZI = (STI+STI)*RAZ
00033       P1R = ZEROR
00034       P1I = ZEROI
00035       P2R = CONER
00036       P2I = CONEI
00037       ACK = (AT+1.0D0)*RAZ
00038       RHO = ACK + DSQRT(ACK*ACK-1.0D0)
00039       RHO2 = RHO*RHO
00040       TST = (RHO2+RHO2)/((RHO2-1.0D0)*(RHO-1.0D0))
00041       TST = TST/TOL
00042 C-----------------------------------------------------------------------
00043 C     COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
00044 C-----------------------------------------------------------------------
00045       AK = AT
00046       DO 10 I=1,80
00047         PTR = P2R
00048         PTI = P2I
00049         P2R = P1R - (CKR*PTR-CKI*PTI)
00050         P2I = P1I - (CKI*PTR+CKR*PTI)
00051         P1R = PTR
00052         P1I = PTI
00053         CKR = CKR + RZR
00054         CKI = CKI + RZI
00055         AP = XZABS(P2R,P2I)
00056         IF (AP.GT.TST*AK*AK) GO TO 20
00057         AK = AK + 1.0D0
00058    10 CONTINUE
00059       GO TO 110
00060    20 CONTINUE
00061       I = I + 1
00062       K = 0
00063       IF (INU.LT.IAZ) GO TO 40
00064 C-----------------------------------------------------------------------
00065 C     COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
00066 C-----------------------------------------------------------------------
00067       P1R = ZEROR
00068       P1I = ZEROI
00069       P2R = CONER
00070       P2I = CONEI
00071       AT = DBLE(FLOAT(INU)) + 1.0D0
00072       STR = ZR*RAZ
00073       STI = -ZI*RAZ
00074       CKR = STR*AT*RAZ
00075       CKI = STI*AT*RAZ
00076       ACK = AT*RAZ
00077       TST = DSQRT(ACK/TOL)
00078       ITIME = 1
00079       DO 30 K=1,80
00080         PTR = P2R
00081         PTI = P2I
00082         P2R = P1R - (CKR*PTR-CKI*PTI)
00083         P2I = P1I - (CKR*PTI+CKI*PTR)
00084         P1R = PTR
00085         P1I = PTI
00086         CKR = CKR + RZR
00087         CKI = CKI + RZI
00088         AP = XZABS(P2R,P2I)
00089         IF (AP.LT.TST) GO TO 30
00090         IF (ITIME.EQ.2) GO TO 40
00091         ACK = XZABS(CKR,CKI)
00092         FLAM = ACK + DSQRT(ACK*ACK-1.0D0)
00093         FKAP = AP/XZABS(P1R,P1I)
00094         RHO = DMIN1(FLAM,FKAP)
00095         TST = TST*DSQRT(RHO/(RHO*RHO-1.0D0))
00096         ITIME = 2
00097    30 CONTINUE
00098       GO TO 110
00099    40 CONTINUE
00100 C-----------------------------------------------------------------------
00101 C     BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
00102 C-----------------------------------------------------------------------
00103       K = K + 1
00104       KK = MAX0(I+IAZ,K+INU)
00105       FKK = DBLE(FLOAT(KK))
00106       P1R = ZEROR
00107       P1I = ZEROI
00108 C-----------------------------------------------------------------------
00109 C     SCALE P2 AND SUM BY SCLE
00110 C-----------------------------------------------------------------------
00111       P2R = SCLE
00112       P2I = ZEROI
00113       FNF = FNU - DBLE(FLOAT(IFNU))
00114       TFNF = FNF + FNF
00115       BK = DGAMLN(FKK+TFNF+1.0D0,IDUM) - DGAMLN(FKK+1.0D0,IDUM) -
00116      * DGAMLN(TFNF+1.0D0,IDUM)
00117       BK = DEXP(BK)
00118       SUMR = ZEROR
00119       SUMI = ZEROI
00120       KM = KK - INU
00121       DO 50 I=1,KM
00122         PTR = P2R
00123         PTI = P2I
00124         P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
00125         P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI)
00126         P1R = PTR
00127         P1I = PTI
00128         AK = 1.0D0 - TFNF/(FKK+TFNF)
00129         ACK = BK*AK
00130         SUMR = SUMR + (ACK+BK)*P1R
00131         SUMI = SUMI + (ACK+BK)*P1I
00132         BK = ACK
00133         FKK = FKK - 1.0D0
00134    50 CONTINUE
00135       YR(N) = P2R
00136       YI(N) = P2I
00137       IF (N.EQ.1) GO TO 70
00138       DO 60 I=2,N
00139         PTR = P2R
00140         PTI = P2I
00141         P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
00142         P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI)
00143         P1R = PTR
00144         P1I = PTI
00145         AK = 1.0D0 - TFNF/(FKK+TFNF)
00146         ACK = BK*AK
00147         SUMR = SUMR + (ACK+BK)*P1R
00148         SUMI = SUMI + (ACK+BK)*P1I
00149         BK = ACK
00150         FKK = FKK - 1.0D0
00151         M = N - I + 1
00152         YR(M) = P2R
00153         YI(M) = P2I
00154    60 CONTINUE
00155    70 CONTINUE
00156       IF (IFNU.LE.0) GO TO 90
00157       DO 80 I=1,IFNU
00158         PTR = P2R
00159         PTI = P2I
00160         P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
00161         P2I = P1I + (FKK+FNF)*(RZR*PTI+RZI*PTR)
00162         P1R = PTR
00163         P1I = PTI
00164         AK = 1.0D0 - TFNF/(FKK+TFNF)
00165         ACK = BK*AK
00166         SUMR = SUMR + (ACK+BK)*P1R
00167         SUMI = SUMI + (ACK+BK)*P1I
00168         BK = ACK
00169         FKK = FKK - 1.0D0
00170    80 CONTINUE
00171    90 CONTINUE
00172       PTR = ZR
00173       PTI = ZI
00174       IF (KODE.EQ.2) PTR = ZEROR
00175       CALL XZLOG(RZR, RZI, STR, STI, IDUM)
00176       P1R = -FNF*STR + PTR
00177       P1I = -FNF*STI + PTI
00178       AP = DGAMLN(1.0D0+FNF,IDUM)
00179       PTR = P1R - AP
00180       PTI = P1I
00181 C-----------------------------------------------------------------------
00182 C     THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
00183 C     IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
00184 C-----------------------------------------------------------------------
00185       P2R = P2R + SUMR
00186       P2I = P2I + SUMI
00187       AP = XZABS(P2R,P2I)
00188       P1R = 1.0D0/AP
00189       CALL XZEXP(PTR, PTI, STR, STI)
00190       CKR = STR*P1R
00191       CKI = STI*P1R
00192       PTR = P2R*P1R
00193       PTI = -P2I*P1R
00194       CALL ZMLT(CKR, CKI, PTR, PTI, CNORMR, CNORMI)
00195       DO 100 I=1,N
00196         STR = YR(I)*CNORMR - YI(I)*CNORMI
00197         YI(I) = YR(I)*CNORMI + YI(I)*CNORMR
00198         YR(I) = STR
00199   100 CONTINUE
00200       RETURN
00201   110 CONTINUE
00202       NZ=-2
00203       RETURN
00204       END
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines