GNU Octave: libcruft/amos/ckscl.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CKSCL(ZR, FNU, N, Y, NZ, RZ, ASCLE, TOL, ELIM)
00002 C***BEGIN PROLOGUE  CKSCL
00003 C***REFER TO  CBKNU,CUNK1,CUNK2
00004 C
00005 C     SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE
00006 C     ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN
00007 C     RETURN WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL.
00008 C
00009 C***ROUTINES CALLED  CUCHK
00010 C***END PROLOGUE  CKSCL
00011       COMPLEX CK, CS, CY, CZERO, RZ, S1, S2, Y, ZR, ZD, CELM
00012       REAL AA, ASCLE, ACS, AS, CSI, CSR, ELIM, FN, FNU, TOL, XX, ZRI,
00013      * ELM, ALAS, HELIM
00014       INTEGER I, IC, K, KK, N, NN, NW, NZ
00015       DIMENSION Y(N), CY(2)
00016       DATA CZERO / (0.0E0,0.0E0) /
00017 C
00018       NZ = 0
00019       IC = 0
00020       XX = REAL(ZR)
00021       NN = MIN0(2,N)
00022       DO 10 I=1,NN
00023         S1 = Y(I)
00024         CY(I) = S1
00025         AS = CABS(S1)
00026         ACS = -XX + ALOG(AS)
00027         NZ = NZ + 1
00028         Y(I) = CZERO
00029         IF (ACS.LT.(-ELIM)) GO TO 10
00030         CS = -ZR + CLOG(S1)
00031         CSR = REAL(CS)
00032         CSI = AIMAG(CS)
00033         AA = EXP(CSR)/TOL
00034         CS = CMPLX(AA,0.0E0)*CMPLX(COS(CSI),SIN(CSI))
00035         CALL CUCHK(CS, NW, ASCLE, TOL)
00036         IF (NW.NE.0) GO TO 10
00037         Y(I) = CS
00038         NZ = NZ - 1
00039         IC = I
00040    10 CONTINUE
00041       IF (N.EQ.1) RETURN
00042       IF (IC.GT.1) GO TO 20
00043       Y(1) = CZERO
00044       NZ = 2
00045    20 CONTINUE
00046       IF (N.EQ.2) RETURN
00047       IF (NZ.EQ.0) RETURN
00048       FN = FNU + 1.0E0
00049       CK = CMPLX(FN,0.0E0)*RZ
00050       S1 = CY(1)
00051       S2 = CY(2)
00052       HELIM = 0.5E0*ELIM
00053       ELM = EXP(-ELIM)
00054       CELM = CMPLX(ELM,0.0E0)
00055       ZRI =AIMAG(ZR)
00056       ZD = ZR
00057 C
00058 C     FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF
00059 C     S2 GETS LARGER THAN EXP(ELIM/2)
00060 C
00061       DO 30 I=3,N
00062         KK = I
00063         CS = S2
00064         S2 = CK*S2 + S1
00065         S1 = CS
00066         CK = CK + RZ
00067         AS = CABS(S2)
00068         ALAS = ALOG(AS)
00069         ACS = -XX + ALAS
00070         NZ = NZ + 1
00071         Y(I) = CZERO
00072         IF (ACS.LT.(-ELIM)) GO TO 25
00073         CS = -ZD + CLOG(S2)
00074         CSR = REAL(CS)
00075         CSI = AIMAG(CS)
00076         AA = EXP(CSR)/TOL
00077         CS = CMPLX(AA,0.0E0)*CMPLX(COS(CSI),SIN(CSI))
00078         CALL CUCHK(CS, NW, ASCLE, TOL)
00079         IF (NW.NE.0) GO TO 25
00080         Y(I) = CS
00081         NZ = NZ - 1
00082         IF (IC.EQ.(KK-1)) GO TO 40
00083         IC = KK
00084         GO TO 30
00085    25   CONTINUE
00086         IF(ALAS.LT.HELIM) GO TO 30
00087         XX = XX-ELIM
00088         S1 = S1*CELM
00089         S2 = S2*CELM
00090         ZD = CMPLX(XX,ZRI)
00091    30 CONTINUE
00092       NZ = N
00093       IF(IC.EQ.N) NZ=N-1
00094       GO TO 45
00095    40 CONTINUE
00096       NZ = KK - 2
00097    45 CONTINUE
00098       DO 50 K=1,NZ
00099         Y(K) = CZERO
00100    50 CONTINUE
00101       RETURN
00102       END