zkscl.f

Go to the documentation of this file.

      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,XZABS,XZLOG
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, xzabs,
     * 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 = xzabs(s1r,s1i)
        acs = -zrr + dlog(as)
        nz = nz + 1
        yr(i) = zeror
        yi(i) = zeroi
        IF (acs.LT.(-elim)) go to 10
        CALL xzlog(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 = xzabs(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 xzlog(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