cbknu.f

Go to the documentation of this file.

      SUBROUTINE cbknu(Z, FNU, KODE, N, Y, NZ, TOL, ELIM, ALIM)
C***BEGIN PROLOGUE  CBKNU
C***REFER TO  CBESI,CBESK,CAIRY,CBESH
C
C     CBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE
C
C***ROUTINES CALLED  CKSCL,CSHCH,GAMLN,I1MACH,R1MACH,CUCHK
C***END PROLOGUE  CBKNU
C
      COMPLEX CCH, CK, COEF, CONE, CRSC, CS, CSCL, CSH, CSR, CSS, CTWO,
     * CZ, CZERO, F, FMU, P, PT, P1, P2, Q, RZ, SMU, ST, S1, S2, Y, Z,
     * ZD, CELM, CY
      REAL AA, AK, ALIM, ASCLE, A1, A2, BB, BK, BRY, CAZ, CC, DNU,
     * DNU2, ELIM, ETEST, FC, FHS, FK, FKS, FNU, FPI, G1, G2, HPI, PI,
     * P2I, P2M, P2R, RK, RTHPI, R1, S, SPI, TM, TOL, TTH, T1, T2, XX,
     * YY, GAMLN, R1MACH, HELIM, ELM, XD, YD, ALAS, AS
      INTEGER I, IDUM, IFLAG, INU, K, KFLAG, KK, KMAX, KODE, KODED, N,
     * NZ, I1MACH, NW, J, IC, INUB
      dimension bry(3), cc(8), css(3), csr(3), y(n), cy(2)
C
      DATA kmax / 30 /
      DATA r1 / 2.0e0 /
      DATA czero,cone,ctwo /(0.0e0,0.0e0),(1.0e0,0.0e0),(2.0e0,0.0e0)/
C
      DATA pi, rthpi, spi ,hpi, fpi, tth /
     1     3.14159265358979324e0,       1.25331413731550025e0,
     2     1.90985931710274403e0,       1.57079632679489662e0,
     3     1.89769999331517738e0,       6.66666666666666666e-01/
C
      DATA cc(1), cc(2), cc(3), cc(4), cc(5), cc(6), cc(7), cc(8)/
     1     5.77215664901532861e-01,    -4.20026350340952355e-02,
     2    -4.21977345555443367e-02,     7.21894324666309954e-03,
     3    -2.15241674114950973e-04,    -2.01348547807882387e-05,
     4     1.13302723198169588e-06,     6.11609510448141582e-09/
C
      xx = real(z)
      yy = aimag(z)
      caz = cabs(z)
      cscl = cmplx(1.0e0/tol,0.0e0)
      crsc = cmplx(tol,0.0e0)
      css(1) = cscl
      css(2) = cone
      css(3) = crsc
      csr(1) = crsc
      csr(2) = cone
      csr(3) = cscl
      bry(1) = 1.0e+3*r1mach(1)/tol
      bry(2) = 1.0e0/bry(1)
      bry(3) = r1mach(2)
      nz = 0
      iflag = 0
      koded = kode
      rz = ctwo/z
      inu = int(fnu+0.5e0)
      dnu = fnu - float(inu)
      IF (abs(dnu).EQ.0.5e0) GO TO 110
      dnu2 = 0.0e0
      IF (abs(dnu).GT.tol) dnu2 = dnu*dnu
      IF (caz.GT.r1) GO TO 110
C-----------------------------------------------------------------------
C     SERIES FOR CABS(Z).LE.R1
C-----------------------------------------------------------------------
      fc = 1.0e0
      smu = clog(rz)
      fmu = smu*cmplx(dnu,0.0e0)
      CALL cshch(fmu, csh, cch)
      IF (dnu.EQ.0.0e0) GO TO 10
      fc = dnu*pi
      fc = fc/sin(fc)
      smu = csh*cmplx(1.0e0/dnu,0.0e0)
   10 CONTINUE
      a2 = 1.0e0 + 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(-gamln(a2,idum))
      t1 = 1.0e0/(t2*fc)
      IF (abs(dnu).GT.0.1e0) GO TO 40
C-----------------------------------------------------------------------
C     SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
C-----------------------------------------------------------------------
      ak = 1.0e0
      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 = 0.5e0*(t1+t2)*fc
      g1 = g1*fc
      f = cmplx(g1,0.0e0)*cch + smu*cmplx(g2,0.0e0)
      pt = cexp(fmu)
      p = cmplx(0.5e0/t2,0.0e0)*pt
      q = cmplx(0.5e0/t1,0.0e0)/pt
      s1 = f
      s2 = p
      ak = 1.0e0
      a1 = 1.0e0
      ck = cone
      bk = 1.0e0 - 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
      cz = z*z*cmplx(0.25e0,0.0e0)
      t1 = 0.25e0*caz*caz
   60 CONTINUE
      f = (f*cmplx(ak,0.0e0)+p+q)*cmplx(1.0e0/bk,0.0e0)
      p = p*cmplx(1.0e0/(ak-dnu),0.0e0)
      q = q*cmplx(1.0e0/(ak+dnu),0.0e0)
      rk = 1.0e0/ak
      ck = ck*cz*cmplx(rk,0.0)
      s1 = s1 + ck*f
      a1 = a1*t1*rk
      bk = bk + ak + ak + 1.0e0
      ak = ak + 1.0e0
      IF (a1.GT.tol) GO TO 60
   70 CONTINUE
      y(1) = s1
      IF (koded.EQ.1) RETURN
      y(1) = s1*cexp(z)
      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
      cz = z*z*cmplx(0.25e0,0.0e0)
      t1 = 0.25e0*caz*caz
   90 CONTINUE
      f = (f*cmplx(ak,0.0e0)+p+q)*cmplx(1.0e0/bk,0.0e0)
      p = p*cmplx(1.0e0/(ak-dnu),0.0e0)
      q = q*cmplx(1.0e0/(ak+dnu),0.0e0)
      rk = 1.0e0/ak
      ck = ck*cz*cmplx(rk,0.0e0)
      s1 = s1 + ck*f
      s2 = s2 + ck*(p-f*cmplx(ak,0.0e0))
      a1 = a1*t1*rk
      bk = bk + ak + ak + 1.0e0
      ak = ak + 1.0e0
      IF (a1.GT.tol) GO TO 90
  100 CONTINUE
      kflag = 2
      bk = real(smu)
      a1 = fnu + 1.0e0
      ak = a1*abs(bk)
      IF (ak.GT.alim) kflag = 3
      p2 = s2*css(kflag)
      s2 = p2*rz
      s1 = s1*css(kflag)
      IF (koded.EQ.1) GO TO 210
      f = cexp(z)
      s1 = s1*f
      s2 = s2*f
      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
      coef = cmplx(rthpi,0.0e0)/csqrt(z)
      kflag = 2
      IF (koded.EQ.2) GO TO 120
      IF (xx.GT.alim) GO TO 290
C     BLANK LINE
      a1 = exp(-xx)*real(css(kflag))
      pt = cmplx(a1,0.0e0)*cmplx(cos(yy),-sin(yy))
      coef = coef*pt
  120 CONTINUE
      IF (abs(dnu).EQ.0.5e0) GO TO 300
C-----------------------------------------------------------------------
C     MILLER ALGORITHM FOR CABS(Z).GT.R1
C-----------------------------------------------------------------------
      ak = cos(pi*dnu)
      ak = abs(ak)
      IF (ak.EQ.0.0e0) GO TO 300
      fhs = abs(0.25e0-dnu2)
      IF (fhs.EQ.0.0e0) GO TO 300
C-----------------------------------------------------------------------
C     COMPUTE R2=F(E). IF CABS(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(11))=
C     TOL WHERE B IS THE BASE OF THE ARITHMETIC.
C-----------------------------------------------------------------------
      t1 = float(i1mach(11)-1)*r1mach(5)*3.321928094e0
      t1 = amax1(t1,12.0e0)
      t1 = amin1(t1,60.0e0)
      t2 = tth*t1 - 6.0e0
      IF (xx.NE.0.0e0) GO TO 130
      t1 = hpi
      GO TO 140
  130 CONTINUE
      t1 = atan(yy/xx)
      t1 = abs(t1)
  140 CONTINUE
      IF (t2.GT.caz) GO TO 170
C-----------------------------------------------------------------------
C     FORWARD RECURRENCE LOOP WHEN CABS(Z).GE.R2
C-----------------------------------------------------------------------
      etest = ak/(pi*caz*tol)
      fk = 1.0e0
      IF (etest.LT.1.0e0) GO TO 180
      fks = 2.0e0
      rk = caz + caz + 2.0e0
      a1 = 0.0e0
      a2 = 1.0e0
      DO 150 i=1,kmax
        ak = fhs/fks
        bk = rk/(fk+1.0e0)
        tm = a2
        a2 = bk*a2 - ak*a1
        a1 = tm
        rk = rk + 2.0e0
        fks = fks + fk + fk + 2.0e0
        fhs = fhs + fk + fk
        fk = fk + 1.0e0
        tm = abs(a2)*fk
        IF (etest.LT.tm) GO TO 160
  150 CONTINUE
      GO TO 310
  160 CONTINUE
      fk = fk + spi*t1*sqrt(t2/caz)
      fhs = abs(0.25e0-dnu2)
      GO TO 180
  170 CONTINUE
C-----------------------------------------------------------------------
C     COMPUTE BACKWARD INDEX K FOR CABS(Z).LT.R2
C-----------------------------------------------------------------------
      a2 = sqrt(caz)
      ak = fpi*ak/(tol*sqrt(a2))
      aa = 3.0e0*t1/(1.0e0+caz)
      bb = 14.7e0*t1/(28.0e0+caz)
      ak = (alog(ak)+caz*cos(aa)/(1.0e0+0.008e0*caz))/cos(bb)
      fk = 0.12125e0*ak*ak/caz + 1.5e0
  180 CONTINUE
      k = int(fk)
C-----------------------------------------------------------------------
C     BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM
C-----------------------------------------------------------------------
      fk = float(k)
      fks = fk*fk
      p1 = czero
      p2 = cmplx(tol,0.0e0)
      cs = p2
      DO 190 i=1,k
        a1 = fks - fk
        a2 = (fks+fk)/(a1+fhs)
        rk = 2.0e0/(fk+1.0e0)
        t1 = (fk+xx)*rk
        t2 = yy*rk
        pt = p2
        p2 = (p2*cmplx(t1,t2)-p1)*cmplx(a2,0.0e0)
        p1 = pt
        cs = cs + p2
        fks = a1 - fk + 1.0e0
        fk = fk - 1.0e0
  190 CONTINUE
C-----------------------------------------------------------------------
C     COMPUTE (P2/CS)=(P2/CABS(CS))*(CONJG(CS)/CABS(CS)) FOR BETTER
C     SCALING
C-----------------------------------------------------------------------
      tm = cabs(cs)
      pt = cmplx(1.0e0/tm,0.0e0)
      s1 = pt*p2
      cs = conjg(cs)*pt
      s1 = coef*s1*cs
      IF (inu.GT.0 .OR. n.GT.1) GO TO 200
      zd = z
      IF(iflag.EQ.1) GO TO 270
      GO TO 240
  200 CONTINUE
C-----------------------------------------------------------------------
C     COMPUTE P1/P2=(P1/CABS(P2)*CONJG(P2)/CABS(P2) FOR SCALING
C-----------------------------------------------------------------------
      tm = cabs(p2)
      pt = cmplx(1.0e0/tm,0.0e0)
      p1 = pt*p1
      p2 = conjg(p2)*pt
      pt = p1*p2
      s2 = s1*(cone+(cmplx(dnu+0.5e0,0.0e0)-pt)/z)
C-----------------------------------------------------------------------
C     FORWARD RECURSION ON THE THREE TERM RECURSION RELATION WITH
C     SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3
C-----------------------------------------------------------------------
  210 CONTINUE
      ck = cmplx(dnu+1.0e0,0.0e0)*rz
      IF (n.EQ.1) inu = inu - 1
      IF (inu.GT.0) GO TO 220
      IF (n.EQ.1) s1=s2
      zd = z
      IF(iflag.EQ.1) GO TO 270
      GO TO 240
  220 CONTINUE
      inub = 1
      IF (iflag.EQ.1) GO TO 261
  225 CONTINUE
      p1 = csr(kflag)
      ascle = bry(kflag)
      DO 230 i=inub,inu
        st = s2
        s2 = ck*s2 + s1
        s1 = st
        ck = ck + rz
        IF (kflag.GE.3) GO TO 230
        p2 = s2*p1
        p2r = real(p2)
        p2i = aimag(p2)
        p2r = abs(p2r)
        p2i = abs(p2i)
        p2m = amax1(p2r,p2i)
        IF (p2m.LE.ascle) GO TO 230
        kflag = kflag + 1
        ascle = bry(kflag)
        s1 = s1*p1
        s2 = p2
        s1 = s1*css(kflag)
        s2 = s2*css(kflag)
        p1 = csr(kflag)
  230 CONTINUE
      IF (n.EQ.1) s1 = s2
  240 CONTINUE
      y(1) = s1*csr(kflag)
      IF (n.EQ.1) RETURN
      y(2) = s2*csr(kflag)
      IF (n.EQ.2) RETURN
      kk = 2
  250 CONTINUE
      kk = kk + 1
      IF (kk.GT.n) RETURN
      p1 = csr(kflag)
      ascle = bry(kflag)
      DO 260 i=kk,n
        p2 = s2
        s2 = ck*s2 + s1
        s1 = p2
        ck = ck + rz
        p2 = s2*p1
        y(i) = p2
        IF (kflag.GE.3) GO TO 260
        p2r = real(p2)
        p2i = aimag(p2)
        p2r = abs(p2r)
        p2i = abs(p2i)
        p2m = amax1(p2r,p2i)
        IF (p2m.LE.ascle) GO TO 260
        kflag = kflag + 1
        ascle = bry(kflag)
        s1 = s1*p1
        s2 = p2
        s1 = s1*css(kflag)
        s2 = s2*css(kflag)
        p1 = csr(kflag)
  260 CONTINUE
      RETURN
C-----------------------------------------------------------------------
C     IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW
C-----------------------------------------------------------------------
  261 CONTINUE
      helim = 0.5e0*elim
      elm = exp(-elim)
      celm = cmplx(elm,0.0)
      ascle = bry(1)
      zd = z
      xd = xx
      yd = yy
      ic = -1
      j = 2
      DO 262 i=1,inu
        st = s2
        s2 = ck*s2+s1
        s1 = st
        ck = ck+rz
        as = cabs(s2)
        alas = alog(as)
        p2r = -xd+alas
        IF(p2r.LT.(-elim)) GO TO 263
        p2 = -zd+clog(s2)
        p2r = real(p2)
        p2i = aimag(p2)
        p2m = exp(p2r)/tol
        p1 = cmplx(p2m,0.0e0)*cmplx(cos(p2i),sin(p2i))
        CALL cuchk(p1,nw,ascle,tol)
        IF(nw.NE.0) GO TO 263
        j=3-j
        cy(j) = p1
        IF(ic.EQ.(i-1)) GO TO 264
        ic = i
        GO TO 262
  263   CONTINUE
        IF(alas.LT.helim) GO TO 262
        xd = xd-elim
        s1 = s1*celm
        s2 = s2*celm
        zd = cmplx(xd,yd)
  262 CONTINUE
      IF(n.EQ.1) s1 = s2
      GO TO 270
  264 CONTINUE
      kflag = 1
      inub = i+1
      s2 = cy(j)
      j = 3 - j
      s1 = cy(j)
      IF(inub.LE.inu) GO TO 225
      IF(n.EQ.1) s1 = s2
      GO TO 240
  270 CONTINUE
      y(1) = s1
      IF (n.EQ.1) GO TO 280
      y(2) = s2
  280 CONTINUE
      ascle = bry(1)
      CALL ckscl(zd, fnu, n, y, nz, rz, ascle, tol, elim)
      inu = n - nz
      IF (inu.LE.0) RETURN
      kk = nz + 1
      s1 = y(kk)
      y(kk) = s1*csr(1)
      IF (inu.EQ.1) RETURN
      kk = nz + 2
      s2 = y(kk)
      y(kk) = s2*csr(1)
      IF (inu.EQ.2) RETURN
      t2 = fnu + float(kk-1)
      ck = cmplx(t2,0.0e0)*rz
      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
      s1 = coef
      s2 = coef
      GO TO 210
  310 CONTINUE
      nz=-2
      RETURN
      END