123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204 |
- SUBROUTINE ZMLRI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL)
- C***BEGIN PROLOGUE ZMLRI
- C***REFER TO ZBESI,ZBESK
- C
- C ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE
- C MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES.
- C
- C***ROUTINES CALLED DGAMLN,D1MACH,ZABS,ZEXP,ZLOG,ZMLT
- C***END PROLOGUE ZMLRI
- C COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z
- DOUBLE PRECISION ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI,
- * CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, FNU, PTI, PTR, P1I,
- * P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI,
- * SUMR, TFNF, TOL, TST, YI, YR, ZEROI, ZEROR, ZI, ZR, DGAMLN,
- * D1MACH, ZABS
- INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N, NZ
- DIMENSION YR(N), YI(N)
- DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
- SCLE = D1MACH(1)/TOL
- NZ=0
- AZ = ZABS(COMPLEX(ZR,ZI))
- IAZ = INT(SNGL(AZ))
- IFNU = INT(SNGL(FNU))
- INU = IFNU + N - 1
- AT = DBLE(FLOAT(IAZ)) + 1.0D0
- RAZ = 1.0D0/AZ
- STR = ZR*RAZ
- STI = -ZI*RAZ
- CKR = STR*AT*RAZ
- CKI = STI*AT*RAZ
- RZR = (STR+STR)*RAZ
- RZI = (STI+STI)*RAZ
- P1R = ZEROR
- P1I = ZEROI
- P2R = CONER
- P2I = CONEI
- ACK = (AT+1.0D0)*RAZ
- RHO = ACK + DSQRT(ACK*ACK-1.0D0)
- RHO2 = RHO*RHO
- TST = (RHO2+RHO2)/((RHO2-1.0D0)*(RHO-1.0D0))
- TST = TST/TOL
- C-----------------------------------------------------------------------
- C COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
- C-----------------------------------------------------------------------
- AK = AT
- DO 10 I=1,80
- PTR = P2R
- PTI = P2I
- P2R = P1R - (CKR*PTR-CKI*PTI)
- P2I = P1I - (CKI*PTR+CKR*PTI)
- P1R = PTR
- P1I = PTI
- CKR = CKR + RZR
- CKI = CKI + RZI
- AP = ZABS(COMPLEX(P2R,P2I))
- IF (AP.GT.TST*AK*AK) GO TO 20
- AK = AK + 1.0D0
- 10 CONTINUE
- GO TO 110
- 20 CONTINUE
- I = I + 1
- K = 0
- IF (INU.LT.IAZ) GO TO 40
- C-----------------------------------------------------------------------
- C COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
- C-----------------------------------------------------------------------
- P1R = ZEROR
- P1I = ZEROI
- P2R = CONER
- P2I = CONEI
- AT = DBLE(FLOAT(INU)) + 1.0D0
- STR = ZR*RAZ
- STI = -ZI*RAZ
- CKR = STR*AT*RAZ
- CKI = STI*AT*RAZ
- ACK = AT*RAZ
- TST = DSQRT(ACK/TOL)
- ITIME = 1
- DO 30 K=1,80
- PTR = P2R
- PTI = P2I
- P2R = P1R - (CKR*PTR-CKI*PTI)
- P2I = P1I - (CKR*PTI+CKI*PTR)
- P1R = PTR
- P1I = PTI
- CKR = CKR + RZR
- CKI = CKI + RZI
- AP = ZABS(COMPLEX(P2R,P2I))
- IF (AP.LT.TST) GO TO 30
- IF (ITIME.EQ.2) GO TO 40
- ACK = ZABS(COMPLEX(CKR,CKI))
- FLAM = ACK + DSQRT(ACK*ACK-1.0D0)
- FKAP = AP/ZABS(COMPLEX(P1R,P1I))
- RHO = DMIN1(FLAM,FKAP)
- TST = TST*DSQRT(RHO/(RHO*RHO-1.0D0))
- ITIME = 2
- 30 CONTINUE
- GO TO 110
- 40 CONTINUE
- C-----------------------------------------------------------------------
- C BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
- C-----------------------------------------------------------------------
- K = K + 1
- KK = MAX0(I+IAZ,K+INU)
- FKK = DBLE(FLOAT(KK))
- P1R = ZEROR
- P1I = ZEROI
- C-----------------------------------------------------------------------
- C SCALE P2 AND SUM BY SCLE
- C-----------------------------------------------------------------------
- P2R = SCLE
- P2I = ZEROI
- FNF = FNU - DBLE(FLOAT(IFNU))
- TFNF = FNF + FNF
- BK = DGAMLN(FKK+TFNF+1.0D0,IDUM) - DGAMLN(FKK+1.0D0,IDUM) -
- * DGAMLN(TFNF+1.0D0,IDUM)
- BK = DEXP(BK)
- SUMR = ZEROR
- SUMI = ZEROI
- KM = KK - INU
- DO 50 I=1,KM
- PTR = P2R
- PTI = P2I
- P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
- P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI)
- P1R = PTR
- P1I = PTI
- AK = 1.0D0 - TFNF/(FKK+TFNF)
- ACK = BK*AK
- SUMR = SUMR + (ACK+BK)*P1R
- SUMI = SUMI + (ACK+BK)*P1I
- BK = ACK
- FKK = FKK - 1.0D0
- 50 CONTINUE
- YR(N) = P2R
- YI(N) = P2I
- IF (N.EQ.1) GO TO 70
- DO 60 I=2,N
- PTR = P2R
- PTI = P2I
- P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
- P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI)
- P1R = PTR
- P1I = PTI
- AK = 1.0D0 - TFNF/(FKK+TFNF)
- ACK = BK*AK
- SUMR = SUMR + (ACK+BK)*P1R
- SUMI = SUMI + (ACK+BK)*P1I
- BK = ACK
- FKK = FKK - 1.0D0
- M = N - I + 1
- YR(M) = P2R
- YI(M) = P2I
- 60 CONTINUE
- 70 CONTINUE
- IF (IFNU.LE.0) GO TO 90
- DO 80 I=1,IFNU
- PTR = P2R
- PTI = P2I
- P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
- P2I = P1I + (FKK+FNF)*(RZR*PTI+RZI*PTR)
- P1R = PTR
- P1I = PTI
- AK = 1.0D0 - TFNF/(FKK+TFNF)
- ACK = BK*AK
- SUMR = SUMR + (ACK+BK)*P1R
- SUMI = SUMI + (ACK+BK)*P1I
- BK = ACK
- FKK = FKK - 1.0D0
- 80 CONTINUE
- 90 CONTINUE
- PTR = ZR
- PTI = ZI
- IF (KODE.EQ.2) PTR = ZEROR
- CALL ZLOG(RZR, RZI, STR, STI, IDUM)
- P1R = -FNF*STR + PTR
- P1I = -FNF*STI + PTI
- AP = DGAMLN(1.0D0+FNF,IDUM)
- PTR = P1R - AP
- PTI = P1I
- C-----------------------------------------------------------------------
- C THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
- C IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
- C-----------------------------------------------------------------------
- P2R = P2R + SUMR
- P2I = P2I + SUMI
- AP = ZABS(COMPLEX(P2R,P2I))
- P1R = 1.0D0/AP
- CALL ZEXP(PTR, PTI, STR, STI)
- CKR = STR*P1R
- CKI = STI*P1R
- PTR = P2R*P1R
- PTI = -P2I*P1R
- CALL ZMLT(CKR, CKI, PTR, PTI, CNORMR, CNORMI)
- DO 100 I=1,N
- STR = YR(I)*CNORMR - YI(I)*CNORMI
- YI(I) = YR(I)*CNORMI + YI(I)*CNORMR
- YR(I) = STR
- 100 CONTINUE
- RETURN
- 110 CONTINUE
- NZ=-2
- RETURN
- END
|