123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581 |
- *DECK ZBKNU
- SUBROUTINE ZBKNU (ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM,
- + ALIM)
- C***BEGIN PROLOGUE ZBKNU
- C***SUBSIDIARY
- C***PURPOSE Subsidiary to ZAIRY, ZBESH, ZBESI and ZBESK
- C***LIBRARY SLATEC
- C***TYPE ALL (CBKNU-A, ZBKNU-A)
- C***AUTHOR Amos, D. E., (SNL)
- C***DESCRIPTION
- C
- C ZBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE.
- C
- C***SEE ALSO ZAIRY, ZBESH, ZBESI, ZBESK
- C***ROUTINES CALLED D1MACH, DGAMLN, I1MACH, ZABS, ZDIV, ZEXP, ZKSCL,
- C ZLOG, ZMLT, ZSHCH, ZSQRT, ZUCHK
- C***REVISION HISTORY (YYMMDD)
- C 830501 DATE WRITTEN
- C 910415 Prologue converted to Version 4.0 format. (BAB)
- C 930122 Added ZEXP, ZLOG and ZSQRT to EXTERNAL statement. (RWC)
- C***END PROLOGUE ZBKNU
- C
- DOUBLE PRECISION AA, AK, ALIM, ASCLE, A1, A2, BB, BK, BRY, CAZ,
- * CBI, CBR, CC, CCHI, CCHR, CKI, CKR, COEFI, COEFR, CONEI, CONER,
- * CRSCR, CSCLR, CSHI, CSHR, CSI, CSR, CSRR, CSSR, CTWOR,
- * CZEROI, CZEROR, CZI, CZR, DNU, DNU2, DPI, ELIM, ETEST, FC, FHS,
- * FI, FK, FKS, FMUI, FMUR, FNU, FPI, FR, G1, G2, HPI, PI, PR, PTI,
- * PTR, P1I, P1R, P2I, P2M, P2R, QI, QR, RAK, RCAZ, RTHPI, RZI,
- * RZR, R1, S, SMUI, SMUR, SPI, STI, STR, S1I, S1R, S2I, S2R, TM,
- * TOL, TTH, T1, T2, YI, YR, ZI, ZR, DGAMLN, D1MACH, ZABS, ELM,
- * CELMR, ZDR, ZDI, AS, ALAS, HELIM, CYR, CYI
- INTEGER I, IFLAG, INU, K, KFLAG, KK, KMAX, KODE, KODED, N, NZ,
- * IDUM, I1MACH, J, IC, INUB, NW
- DIMENSION YR(N), YI(N), CC(8), CSSR(3), CSRR(3), BRY(3), CYR(2),
- * CYI(2)
- EXTERNAL ZABS, ZEXP, ZLOG, ZSQRT
- C COMPLEX Z,Y,A,B,RZ,SMU,FU,FMU,F,FLRZ,CZ,S1,S2,CSH,CCH
- C COMPLEX CK,P,Q,COEF,P1,P2,CBK,PT,CZERO,CONE,CTWO,ST,EZ,CS,DK
- C
- DATA KMAX / 30 /
- DATA CZEROR,CZEROI,CONER,CONEI,CTWOR,R1/
- 1 0.0D0 , 0.0D0 , 1.0D0 , 0.0D0 , 2.0D0 , 2.0D0 /
- DATA DPI, RTHPI, SPI ,HPI, FPI, TTH /
- 1 3.14159265358979324D0, 1.25331413731550025D0,
- 2 1.90985931710274403D0, 1.57079632679489662D0,
- 3 1.89769999331517738D0, 6.66666666666666666D-01/
- DATA CC(1), CC(2), CC(3), CC(4), CC(5), CC(6), CC(7), CC(8)/
- 1 5.77215664901532861D-01, -4.20026350340952355D-02,
- 2 -4.21977345555443367D-02, 7.21894324666309954D-03,
- 3 -2.15241674114950973D-04, -2.01348547807882387D-05,
- 4 1.13302723198169588D-06, 6.11609510448141582D-09/
- C***FIRST EXECUTABLE STATEMENT ZBKNU
- CAZ = ZABS(ZR,ZI)
- CSCLR = 1.0D0/TOL
- CRSCR = TOL
- CSSR(1) = CSCLR
- CSSR(2) = 1.0D0
- CSSR(3) = CRSCR
- CSRR(1) = CRSCR
- CSRR(2) = 1.0D0
- CSRR(3) = CSCLR
- BRY(1) = 1.0D+3*D1MACH(1)/TOL
- BRY(2) = 1.0D0/BRY(1)
- BRY(3) = D1MACH(2)
- NZ = 0
- IFLAG = 0
- KODED = KODE
- RCAZ = 1.0D0/CAZ
- STR = ZR*RCAZ
- STI = -ZI*RCAZ
- RZR = (STR+STR)*RCAZ
- RZI = (STI+STI)*RCAZ
- INU = FNU+0.5D0
- DNU = FNU - INU
- IF (ABS(DNU).EQ.0.5D0) GO TO 110
- DNU2 = 0.0D0
- IF (ABS(DNU).GT.TOL) DNU2 = DNU*DNU
- IF (CAZ.GT.R1) GO TO 110
- C-----------------------------------------------------------------------
- C SERIES FOR ABS(Z).LE.R1
- C-----------------------------------------------------------------------
- FC = 1.0D0
- CALL ZLOG(RZR, RZI, SMUR, SMUI, IDUM)
- FMUR = SMUR*DNU
- FMUI = SMUI*DNU
- CALL ZSHCH(FMUR, FMUI, CSHR, CSHI, CCHR, CCHI)
- IF (DNU.EQ.0.0D0) GO TO 10
- FC = DNU*DPI
- FC = FC/SIN(FC)
- SMUR = CSHR/DNU
- SMUI = CSHI/DNU
- 10 CONTINUE
- A2 = 1.0D0 + DNU
- C-----------------------------------------------------------------------
- C GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU)
- C-----------------------------------------------------------------------
- T2 = EXP(-DGAMLN(A2,IDUM))
- T1 = 1.0D0/(T2*FC)
- IF (ABS(DNU).GT.0.1D0) GO TO 40
- C-----------------------------------------------------------------------
- C SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
- C-----------------------------------------------------------------------
- AK = 1.0D0
- S = CC(1)
- DO 20 K=2,8
- AK = AK*DNU2
- TM = CC(K)*AK
- S = S + TM
- IF (ABS(TM).LT.TOL) GO TO 30
- 20 CONTINUE
- 30 G1 = -S
- GO TO 50
- 40 CONTINUE
- G1 = (T1-T2)/(DNU+DNU)
- 50 CONTINUE
- G2 = (T1+T2)*0.5D0
- FR = FC*(CCHR*G1+SMUR*G2)
- FI = FC*(CCHI*G1+SMUI*G2)
- CALL ZEXP(FMUR, FMUI, STR, STI)
- PR = 0.5D0*STR/T2
- PI = 0.5D0*STI/T2
- CALL ZDIV(0.5D0, 0.0D0, STR, STI, PTR, PTI)
- QR = PTR/T1
- QI = PTI/T1
- S1R = FR
- S1I = FI
- S2R = PR
- S2I = PI
- AK = 1.0D0
- A1 = 1.0D0
- CKR = CONER
- CKI = CONEI
- BK = 1.0D0 - DNU2
- IF (INU.GT.0 .OR. N.GT.1) GO TO 80
- C-----------------------------------------------------------------------
- C GENERATE K(FNU,Z), 0.0D0 .LE. FNU .LT. 0.5D0 AND N=1
- C-----------------------------------------------------------------------
- IF (CAZ.LT.TOL) GO TO 70
- CALL ZMLT(ZR, ZI, ZR, ZI, CZR, CZI)
- CZR = 0.25D0*CZR
- CZI = 0.25D0*CZI
- T1 = 0.25D0*CAZ*CAZ
- 60 CONTINUE
- FR = (FR*AK+PR+QR)/BK
- FI = (FI*AK+PI+QI)/BK
- STR = 1.0D0/(AK-DNU)
- PR = PR*STR
- PI = PI*STR
- STR = 1.0D0/(AK+DNU)
- QR = QR*STR
- QI = QI*STR
- STR = CKR*CZR - CKI*CZI
- RAK = 1.0D0/AK
- CKI = (CKR*CZI+CKI*CZR)*RAK
- CKR = STR*RAK
- S1R = CKR*FR - CKI*FI + S1R
- S1I = CKR*FI + CKI*FR + S1I
- A1 = A1*T1*RAK
- BK = BK + AK + AK + 1.0D0
- AK = AK + 1.0D0
- IF (A1.GT.TOL) GO TO 60
- 70 CONTINUE
- YR(1) = S1R
- YI(1) = S1I
- IF (KODED.EQ.1) RETURN
- CALL ZEXP(ZR, ZI, STR, STI)
- CALL ZMLT(S1R, S1I, STR, STI, YR(1), YI(1))
- RETURN
- C-----------------------------------------------------------------------
- C GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE
- C-----------------------------------------------------------------------
- 80 CONTINUE
- IF (CAZ.LT.TOL) GO TO 100
- CALL ZMLT(ZR, ZI, ZR, ZI, CZR, CZI)
- CZR = 0.25D0*CZR
- CZI = 0.25D0*CZI
- T1 = 0.25D0*CAZ*CAZ
- 90 CONTINUE
- FR = (FR*AK+PR+QR)/BK
- FI = (FI*AK+PI+QI)/BK
- STR = 1.0D0/(AK-DNU)
- PR = PR*STR
- PI = PI*STR
- STR = 1.0D0/(AK+DNU)
- QR = QR*STR
- QI = QI*STR
- STR = CKR*CZR - CKI*CZI
- RAK = 1.0D0/AK
- CKI = (CKR*CZI+CKI*CZR)*RAK
- CKR = STR*RAK
- S1R = CKR*FR - CKI*FI + S1R
- S1I = CKR*FI + CKI*FR + S1I
- STR = PR - FR*AK
- STI = PI - FI*AK
- S2R = CKR*STR - CKI*STI + S2R
- S2I = CKR*STI + CKI*STR + S2I
- A1 = A1*T1*RAK
- BK = BK + AK + AK + 1.0D0
- AK = AK + 1.0D0
- IF (A1.GT.TOL) GO TO 90
- 100 CONTINUE
- KFLAG = 2
- A1 = FNU + 1.0D0
- AK = A1*ABS(SMUR)
- IF (AK.GT.ALIM) KFLAG = 3
- STR = CSSR(KFLAG)
- P2R = S2R*STR
- P2I = S2I*STR
- CALL ZMLT(P2R, P2I, RZR, RZI, S2R, S2I)
- S1R = S1R*STR
- S1I = S1I*STR
- IF (KODED.EQ.1) GO TO 210
- CALL ZEXP(ZR, ZI, FR, FI)
- CALL ZMLT(S1R, S1I, FR, FI, S1R, S1I)
- CALL ZMLT(S2R, S2I, FR, FI, S2R, S2I)
- GO TO 210
- C-----------------------------------------------------------------------
- C IFLAG=0 MEANS NO UNDERFLOW OCCURRED
- C IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH
- C KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD
- C RECURSION
- C-----------------------------------------------------------------------
- 110 CONTINUE
- CALL ZSQRT(ZR, ZI, STR, STI)
- CALL ZDIV(RTHPI, CZEROI, STR, STI, COEFR, COEFI)
- KFLAG = 2
- IF (KODED.EQ.2) GO TO 120
- IF (ZR.GT.ALIM) GO TO 290
- C BLANK LINE
- STR = EXP(-ZR)*CSSR(KFLAG)
- STI = -STR*SIN(ZI)
- STR = STR*COS(ZI)
- CALL ZMLT(COEFR, COEFI, STR, STI, COEFR, COEFI)
- 120 CONTINUE
- IF (ABS(DNU).EQ.0.5D0) GO TO 300
- C-----------------------------------------------------------------------
- C MILLER ALGORITHM FOR ABS(Z).GT.R1
- C-----------------------------------------------------------------------
- AK = COS(DPI*DNU)
- AK = ABS(AK)
- IF (AK.EQ.CZEROR) GO TO 300
- FHS = ABS(0.25D0-DNU2)
- IF (FHS.EQ.CZEROR) GO TO 300
- C-----------------------------------------------------------------------
- C COMPUTE R2=F(E). IF ABS(Z).GE.R2, USE FORWARD RECURRENCE TO
- C DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON
- C 12.LE.E.LE.60. E IS COMPUTED FROM 2**(-E)=B**(1-I1MACH(14))=
- C TOL WHERE B IS THE BASE OF THE ARITHMETIC.
- C-----------------------------------------------------------------------
- T1 = I1MACH(14)-1
- T1 = T1*D1MACH(5)*3.321928094D0
- T1 = MAX(T1,12.0D0)
- T1 = MIN(T1,60.0D0)
- T2 = TTH*T1 - 6.0D0
- IF (ZR.NE.0.0D0) GO TO 130
- T1 = HPI
- GO TO 140
- 130 CONTINUE
- T1 = DATAN(ZI/ZR)
- T1 = ABS(T1)
- 140 CONTINUE
- IF (T2.GT.CAZ) GO TO 170
- C-----------------------------------------------------------------------
- C FORWARD RECURRENCE LOOP WHEN ABS(Z).GE.R2
- C-----------------------------------------------------------------------
- ETEST = AK/(DPI*CAZ*TOL)
- FK = CONER
- IF (ETEST.LT.CONER) GO TO 180
- FKS = CTWOR
- CKR = CAZ + CAZ + CTWOR
- P1R = CZEROR
- P2R = CONER
- DO 150 I=1,KMAX
- AK = FHS/FKS
- CBR = CKR/(FK+CONER)
- PTR = P2R
- P2R = CBR*P2R - P1R*AK
- P1R = PTR
- CKR = CKR + CTWOR
- FKS = FKS + FK + FK + CTWOR
- FHS = FHS + FK + FK
- FK = FK + CONER
- STR = ABS(P2R)*FK
- IF (ETEST.LT.STR) GO TO 160
- 150 CONTINUE
- GO TO 310
- 160 CONTINUE
- FK = FK + SPI*T1*SQRT(T2/CAZ)
- FHS = ABS(0.25D0-DNU2)
- GO TO 180
- 170 CONTINUE
- C-----------------------------------------------------------------------
- C COMPUTE BACKWARD INDEX K FOR ABS(Z).LT.R2
- C-----------------------------------------------------------------------
- A2 = SQRT(CAZ)
- AK = FPI*AK/(TOL*SQRT(A2))
- AA = 3.0D0*T1/(1.0D0+CAZ)
- BB = 14.7D0*T1/(28.0D0+CAZ)
- AK = (LOG(AK)+CAZ*COS(AA)/(1.0D0+0.008D0*CAZ))/COS(BB)
- FK = 0.12125D0*AK*AK/CAZ + 1.5D0
- 180 CONTINUE
- C-----------------------------------------------------------------------
- C BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM
- C-----------------------------------------------------------------------
- K = FK
- FK = K
- FKS = FK*FK
- P1R = CZEROR
- P1I = CZEROI
- P2R = TOL
- P2I = CZEROI
- CSR = P2R
- CSI = P2I
- DO 190 I=1,K
- A1 = FKS - FK
- AK = (FKS+FK)/(A1+FHS)
- RAK = 2.0D0/(FK+CONER)
- CBR = (FK+ZR)*RAK
- CBI = ZI*RAK
- PTR = P2R
- PTI = P2I
- P2R = (PTR*CBR-PTI*CBI-P1R)*AK
- P2I = (PTI*CBR+PTR*CBI-P1I)*AK
- P1R = PTR
- P1I = PTI
- CSR = CSR + P2R
- CSI = CSI + P2I
- FKS = A1 - FK + CONER
- FK = FK - CONER
- 190 CONTINUE
- C-----------------------------------------------------------------------
- C COMPUTE (P2/CS)=(P2/ABS(CS))*(CONJG(CS)/ABS(CS)) FOR BETTER
- C SCALING
- C-----------------------------------------------------------------------
- TM = ZABS(CSR,CSI)
- PTR = 1.0D0/TM
- S1R = P2R*PTR
- S1I = P2I*PTR
- CSR = CSR*PTR
- CSI = -CSI*PTR
- CALL ZMLT(COEFR, COEFI, S1R, S1I, STR, STI)
- CALL ZMLT(STR, STI, CSR, CSI, S1R, S1I)
- IF (INU.GT.0 .OR. N.GT.1) GO TO 200
- ZDR = ZR
- ZDI = ZI
- IF(IFLAG.EQ.1) GO TO 270
- GO TO 240
- 200 CONTINUE
- C-----------------------------------------------------------------------
- C COMPUTE P1/P2=(P1/ABS(P2)*CONJG(P2)/ABS(P2) FOR SCALING
- C-----------------------------------------------------------------------
- TM = ZABS(P2R,P2I)
- PTR = 1.0D0/TM
- P1R = P1R*PTR
- P1I = P1I*PTR
- P2R = P2R*PTR
- P2I = -P2I*PTR
- CALL ZMLT(P1R, P1I, P2R, P2I, PTR, PTI)
- STR = DNU + 0.5D0 - PTR
- STI = -PTI
- CALL ZDIV(STR, STI, ZR, ZI, STR, STI)
- STR = STR + 1.0D0
- CALL ZMLT(STR, STI, S1R, S1I, S2R, S2I)
- C-----------------------------------------------------------------------
- C FORWARD RECURSION ON THE THREE TERM RECURSION WITH RELATION WITH
- C SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3
- C-----------------------------------------------------------------------
- 210 CONTINUE
- STR = DNU + 1.0D0
- CKR = STR*RZR
- CKI = STR*RZI
- IF (N.EQ.1) INU = INU - 1
- IF (INU.GT.0) GO TO 220
- IF (N.GT.1) GO TO 215
- S1R = S2R
- S1I = S2I
- 215 CONTINUE
- ZDR = ZR
- ZDI = ZI
- IF(IFLAG.EQ.1) GO TO 270
- GO TO 240
- 220 CONTINUE
- INUB = 1
- IF(IFLAG.EQ.1) GO TO 261
- 225 CONTINUE
- P1R = CSRR(KFLAG)
- ASCLE = BRY(KFLAG)
- DO 230 I=INUB,INU
- STR = S2R
- STI = S2I
- S2R = CKR*STR - CKI*STI + S1R
- S2I = CKR*STI + CKI*STR + S1I
- S1R = STR
- S1I = STI
- CKR = CKR + RZR
- CKI = CKI + RZI
- IF (KFLAG.GE.3) GO TO 230
- P2R = S2R*P1R
- P2I = S2I*P1R
- STR = ABS(P2R)
- STI = ABS(P2I)
- P2M = MAX(STR,STI)
- IF (P2M.LE.ASCLE) GO TO 230
- KFLAG = KFLAG + 1
- ASCLE = BRY(KFLAG)
- S1R = S1R*P1R
- S1I = S1I*P1R
- S2R = P2R
- S2I = P2I
- STR = CSSR(KFLAG)
- S1R = S1R*STR
- S1I = S1I*STR
- S2R = S2R*STR
- S2I = S2I*STR
- P1R = CSRR(KFLAG)
- 230 CONTINUE
- IF (N.NE.1) GO TO 240
- S1R = S2R
- S1I = S2I
- 240 CONTINUE
- STR = CSRR(KFLAG)
- YR(1) = S1R*STR
- YI(1) = S1I*STR
- IF (N.EQ.1) RETURN
- YR(2) = S2R*STR
- YI(2) = S2I*STR
- IF (N.EQ.2) RETURN
- KK = 2
- 250 CONTINUE
- KK = KK + 1
- IF (KK.GT.N) RETURN
- P1R = CSRR(KFLAG)
- ASCLE = BRY(KFLAG)
- DO 260 I=KK,N
- P2R = S2R
- P2I = S2I
- S2R = CKR*P2R - CKI*P2I + S1R
- S2I = CKI*P2R + CKR*P2I + S1I
- S1R = P2R
- S1I = P2I
- CKR = CKR + RZR
- CKI = CKI + RZI
- P2R = S2R*P1R
- P2I = S2I*P1R
- YR(I) = P2R
- YI(I) = P2I
- IF (KFLAG.GE.3) GO TO 260
- STR = ABS(P2R)
- STI = ABS(P2I)
- P2M = MAX(STR,STI)
- IF (P2M.LE.ASCLE) GO TO 260
- KFLAG = KFLAG + 1
- ASCLE = BRY(KFLAG)
- S1R = S1R*P1R
- S1I = S1I*P1R
- S2R = P2R
- S2I = P2I
- STR = CSSR(KFLAG)
- S1R = S1R*STR
- S1I = S1I*STR
- S2R = S2R*STR
- S2I = S2I*STR
- P1R = CSRR(KFLAG)
- 260 CONTINUE
- RETURN
- C-----------------------------------------------------------------------
- C IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW
- C-----------------------------------------------------------------------
- 261 CONTINUE
- HELIM = 0.5D0*ELIM
- ELM = EXP(-ELIM)
- CELMR = ELM
- ASCLE = BRY(1)
- ZDR = ZR
- ZDI = ZI
- IC = -1
- J = 2
- DO 262 I=1,INU
- STR = S2R
- STI = S2I
- S2R = STR*CKR-STI*CKI+S1R
- S2I = STI*CKR+STR*CKI+S1I
- S1R = STR
- S1I = STI
- CKR = CKR+RZR
- CKI = CKI+RZI
- AS = ZABS(S2R,S2I)
- ALAS = LOG(AS)
- P2R = -ZDR+ALAS
- IF(P2R.LT.(-ELIM)) GO TO 263
- CALL ZLOG(S2R,S2I,STR,STI,IDUM)
- P2R = -ZDR+STR
- P2I = -ZDI+STI
- P2M = EXP(P2R)/TOL
- P1R = P2M*COS(P2I)
- P1I = P2M*SIN(P2I)
- CALL ZUCHK(P1R,P1I,NW,ASCLE,TOL)
- IF(NW.NE.0) GO TO 263
- J = 3 - J
- CYR(J) = P1R
- CYI(J) = P1I
- IF(IC.EQ.(I-1)) GO TO 264
- IC = I
- GO TO 262
- 263 CONTINUE
- IF(ALAS.LT.HELIM) GO TO 262
- ZDR = ZDR-ELIM
- S1R = S1R*CELMR
- S1I = S1I*CELMR
- S2R = S2R*CELMR
- S2I = S2I*CELMR
- 262 CONTINUE
- IF(N.NE.1) GO TO 270
- S1R = S2R
- S1I = S2I
- GO TO 270
- 264 CONTINUE
- KFLAG = 1
- INUB = I+1
- S2R = CYR(J)
- S2I = CYI(J)
- J = 3 - J
- S1R = CYR(J)
- S1I = CYI(J)
- IF(INUB.LE.INU) GO TO 225
- IF(N.NE.1) GO TO 240
- S1R = S2R
- S1I = S2I
- GO TO 240
- 270 CONTINUE
- YR(1) = S1R
- YI(1) = S1I
- IF(N.EQ.1) GO TO 280
- YR(2) = S2R
- YI(2) = S2I
- 280 CONTINUE
- ASCLE = BRY(1)
- CALL ZKSCL(ZDR,ZDI,FNU,N,YR,YI,NZ,RZR,RZI,ASCLE,TOL,ELIM)
- INU = N - NZ
- IF (INU.LE.0) RETURN
- KK = NZ + 1
- S1R = YR(KK)
- S1I = YI(KK)
- YR(KK) = S1R*CSRR(1)
- YI(KK) = S1I*CSRR(1)
- IF (INU.EQ.1) RETURN
- KK = NZ + 2
- S2R = YR(KK)
- S2I = YI(KK)
- YR(KK) = S2R*CSRR(1)
- YI(KK) = S2I*CSRR(1)
- IF (INU.EQ.2) RETURN
- T2 = FNU + (KK-1)
- CKR = T2*RZR
- CKI = T2*RZI
- KFLAG = 1
- GO TO 250
- 290 CONTINUE
- C-----------------------------------------------------------------------
- C SCALE BY EXP(Z), IFLAG = 1 CASES
- C-----------------------------------------------------------------------
- KODED = 2
- IFLAG = 1
- KFLAG = 2
- GO TO 120
- C-----------------------------------------------------------------------
- C FNU=HALF ODD INTEGER CASE, DNU=-0.5
- C-----------------------------------------------------------------------
- 300 CONTINUE
- S1R = COEFR
- S1I = COEFI
- S2R = COEFR
- S2I = COEFI
- GO TO 210
- C
- C
- 310 CONTINUE
- NZ=-2
- RETURN
- END
|