ckscl.f

Go to the documentation of this file.

      SUBROUTINE ckscl(ZR, FNU, N, Y, NZ, RZ, ASCLE, TOL, ELIM)
C***BEGIN PROLOGUE  CKSCL
C***REFER TO  CBKNU,CUNK1,CUNK2
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  CUCHK
C***END PROLOGUE  CKSCL
      COMPLEX ck, cs, cy, czero, rz, s1, s2, y, zr, zd, celm
      REAL aa, ascle, acs, as, csi, csr, elim, fn, fnu, tol, xx, zri,
     * elm, alas, helim
      INTEGER i, ic, k, kk, n, nn, nw, nz
      dimension y(n), cy(2)
      DATA czero / (0.0e0,0.0e0) /
C
      nz = 0
      ic = 0
      xx = REAL(zr)
      nn = min0(2,n)
      DO 10 i=1,nn
        s1 = y(i)
        cy(i) = s1
        as = cabs(s1)
        acs = -xx + alog(as)
        nz = nz + 1
        y(i) = czero
        IF (acs.LT.(-elim)) go to 10
        cs = -zr + clog(s1)
        csr = REAL(cs)
        csi = aimag(cs)
        aa = exp(csr)/tol
        cs = cmplx(aa,0.0e0)*cmplx(cos(csi),sin(csi))
        CALL cuchk(cs, nw, ascle, tol)
        IF (nw.NE.0) go to 10
        y(i) = cs
        nz = nz - 1
        ic = i
   10 CONTINUE
      IF (n.EQ.1) RETURN
      IF (ic.GT.1) go to 20
      y(1) = czero
      nz = 2
   20 CONTINUE
      IF (n.EQ.2) RETURN
      IF (nz.EQ.0) RETURN
      fn = fnu + 1.0e0
      ck = cmplx(fn,0.0e0)*rz
      s1 = cy(1)
      s2 = cy(2)
      helim = 0.5e0*elim
      elm = exp(-elim)
      celm = cmplx(elm,0.0e0)
      zri =aimag(zr)
      zd = zr
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
        cs = s2
        s2 = ck*s2 + s1
        s1 = cs
        ck = ck + rz
        as = cabs(s2)
        alas = alog(as)
        acs = -xx + alas
        nz = nz + 1
        y(i) = czero
        IF (acs.LT.(-elim)) go to 25
        cs = -zd + clog(s2)
        csr = REAL(cs)
        csi = aimag(cs)
        aa = exp(csr)/tol
        cs = cmplx(aa,0.0e0)*cmplx(cos(csi),sin(csi))
        CALL cuchk(cs, nw, ascle, tol)
        IF (nw.NE.0) go to 25
        y(i) = cs
        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
        xx = xx-elim
        s1 = s1*celm
        s2 = s2*celm
        zd = cmplx(xx,zri)
   30 CONTINUE
      nz = n
      IF(ic.EQ.n) nz=n-1
      go to 45
   40 CONTINUE
      nz = kk - 2
   45 CONTINUE
      DO 50 k=1,nz
        y(k) = czero
   50 CONTINUE
      RETURN
      END