qpms/amos/zkscl.f

122 lines
3.0 KiB
FortranFixed
Raw Permalink Normal View History

SUBROUTINE ZKSCL(ZRR,ZRI,FNU,N,YR,YI,NZ,RZR,RZI,ASCLE,TOL,ELIM)
C***BEGIN PROLOGUE ZKSCL
C***REFER TO ZBESK
C
C SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE
C ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN
C RETURN WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL.
C
C***ROUTINES CALLED ZUCHK,AZABS,AZLOG
C***END PROLOGUE ZKSCL
C COMPLEX CK,CS,CY,CZERO,RZ,S1,S2,Y,ZR,ZD,CELM
DOUBLE PRECISION ACS, AS, ASCLE, CKI, CKR, CSI, CSR, CYI,
* CYR, ELIM, FN, FNU, RZI, RZR, STR, S1I, S1R, S2I,
* S2R, TOL, YI, YR, ZEROI, ZEROR, ZRI, ZRR, AZABS,
* ZDR, ZDI, CELMR, ELM, HELIM, ALAS
INTEGER I, IC, IDUM, KK, N, NN, NW, NZ
DIMENSION YR(N), YI(N), CYR(2), CYI(2)
DATA ZEROR,ZEROI / 0.0D0 , 0.0D0 /
C
NZ = 0
IC = 0
NN = MIN0(2,N)
DO 10 I=1,NN
S1R = YR(I)
S1I = YI(I)
CYR(I) = S1R
CYI(I) = S1I
AS = AZABS(S1R,S1I)
ACS = -ZRR + DLOG(AS)
NZ = NZ + 1
YR(I) = ZEROR
YI(I) = ZEROI
IF (ACS.LT.(-ELIM)) GO TO 10
CALL AZLOG(S1R, S1I, CSR, CSI, IDUM)
CSR = CSR - ZRR
CSI = CSI - ZRI
STR = DEXP(CSR)/TOL
CSR = STR*DCOS(CSI)
CSI = STR*DSIN(CSI)
CALL ZUCHK(CSR, CSI, NW, ASCLE, TOL)
IF (NW.NE.0) GO TO 10
YR(I) = CSR
YI(I) = CSI
IC = I
NZ = NZ - 1
10 CONTINUE
IF (N.EQ.1) RETURN
IF (IC.GT.1) GO TO 20
YR(1) = ZEROR
YI(1) = ZEROI
NZ = 2
20 CONTINUE
IF (N.EQ.2) RETURN
IF (NZ.EQ.0) RETURN
FN = FNU + 1.0D0
CKR = FN*RZR
CKI = FN*RZI
S1R = CYR(1)
S1I = CYI(1)
S2R = CYR(2)
S2I = CYI(2)
HELIM = 0.5D0*ELIM
ELM = DEXP(-ELIM)
CELMR = ELM
ZDR = ZRR
ZDI = ZRI
C
C FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF
C S2 GETS LARGER THAN EXP(ELIM/2)
C
DO 30 I=3,N
KK = I
CSR = S2R
CSI = S2I
S2R = CKR*CSR - CKI*CSI + S1R
S2I = CKI*CSR + CKR*CSI + S1I
S1R = CSR
S1I = CSI
CKR = CKR + RZR
CKI = CKI + RZI
AS = AZABS(S2R,S2I)
ALAS = DLOG(AS)
ACS = -ZDR + ALAS
NZ = NZ + 1
YR(I) = ZEROR
YI(I) = ZEROI
IF (ACS.LT.(-ELIM)) GO TO 25
CALL AZLOG(S2R, S2I, CSR, CSI, IDUM)
CSR = CSR - ZDR
CSI = CSI - ZDI
STR = DEXP(CSR)/TOL
CSR = STR*DCOS(CSI)
CSI = STR*DSIN(CSI)
CALL ZUCHK(CSR, CSI, NW, ASCLE, TOL)
IF (NW.NE.0) GO TO 25
YR(I) = CSR
YI(I) = CSI
NZ = NZ - 1
IF (IC.EQ.KK-1) GO TO 40
IC = KK
GO TO 30
25 CONTINUE
IF(ALAS.LT.HELIM) GO TO 30
ZDR = ZDR - ELIM
S1R = S1R*CELMR
S1I = S1I*CELMR
S2R = S2R*CELMR
S2I = S2I*CELMR
30 CONTINUE
NZ = N
IF(IC.EQ.N) NZ=N-1
GO TO 45
40 CONTINUE
NZ = KK - 2
45 CONTINUE
DO 50 I=1,NZ
YR(I) = ZEROR
YI(I) = ZEROI
50 CONTINUE
RETURN
END