c/* %W%      latest revision %G% %U% */
      subroutine rcwfn (rho,eta,minl,maxl,fc,fcp,gc,gcp,accur,step)
c *** see r h landau, program lpott, 1981
      implicit  real*8   (a-h,o-z)
      real*8 k,k1,k2,k3,k4,m1,m2,m3,m4
      real*8 rrpace
      dimension fc(*), fcp(*), gc(*), gcp(*)
c *** coulomb wavefunctions calculated at r = rho by the
c *** continued-fraction method of steed   minl,maxl are actual l-values
c *** see barnett feng steed and goldfarb computer physics commun 1974
      pace = step
      acc = accur
      if (pace.lt.100.0) pace = 100.0
      if (acc.lt.1.0e-15.or.acc.gt.1.0e-6) acc = 1.0e-6
      r = rho
      ktr = 1
      lmax = maxl
      lmin1 = minl+1
      xll1 = minl*lmin1
      eta2 = eta*eta
      turn = eta+sqrt(eta2+xll1)
      if (r.lt.turn.and.abs(eta).ge.1.0e-6) ktr = -1
      ktrp = ktr
      go to 20
 10   r = turn
      tf = f
      tfp = fp
      lmax = minl
      ktrp = 1
 20   etar = eta*r
      rho2 = r*r
      pl = lmax+1
      pmx = pl+0.5
c *** continued fraction for fp(maxl)/f(maxl)  xl is f  xlprime is fp **
      fp = eta/pl+pl/r
      dk = etar*2.0
      del = 0.0
      d = 0.0
      f = 1.0
      k = (pl*pl-pl+etar)*(2.0*pl-1.0)
      if (pl*pl+pl+etar.ne.0.0) go to 30
      r = r+1.0e-6
      go to 20
 30   h = (pl*pl+eta2)*(1.0-pl*pl)*rho2
      k = k+dk+pl*pl*6.0
      d = 1.0/(d*h+k)
      del = del*(d*k-1.0)
      if (pl.lt.pmx) del = -r*(pl*pl+eta2)*(pl+1.0)*d/pl
      pl = pl+1.0
      fp = fp+del
      if (d.lt.0.0) f = -f
      if (pl.gt.20000.) go to 110
      if (abs(del/fp).ge.acc) go to 30
      fp = f*fp
      if (lmax.eq.minl) go to 50
      fc(lmax+1) = f
      fcp(lmax+1) = fp
c *** downward recursion to minl for f and fp, arrays gc,gcp are storage
      l = lmax
      do 40 lp=lmin1,lmax
         pl = float(l)
         gc(l+1) = eta/pl+pl/r
         gcp(l+1) = sqrt(eta2+pl*pl)/pl
         fc(l) = (gc(l+1)*fc(l+1)+fcp(l+1))/gcp(l+1)
         fcp(l) = gc(l+1)*fc(l)-gcp(l+1)*fc(l+1)
         l = l-1
 40   continue
      f = fc(lmin1)
      fp = fcp(lmin1)
 50   if (ktrp.eq.-1) go to 10
c *** repeat for r = turn if rho lt turn
c *** now obtain p + i.q for minl from continued fraction (32)
c *** real arithmetic to facilitate conversion to ibm using r*8
      p = 0.0
      q = r-eta
      pl = 0.0
      ar = -(eta2+xll1)
      ai = eta
      br = 2.0*q
      bi = 2.0
      wi = 2.0*eta
      dr = br/(br*br+bi*bi)
      di = -bi/(br*br+bi*bi)
      dp = -(ar*di+ai*dr)
      dq = (ar*dr-ai*di)
 60   p = p+dp
      q = q+dq
      pl = pl+2.0
      ar = ar+pl
      ai = ai+wi
      bi = bi+2.0
      d = ar*dr-ai*di+br
      di = ai*dr+ar*di+bi
      t = 1.0/(d*d+di*di)
      dr = t*d
      di = -t*di
      h = br*dr-bi*di-1.0
      k = bi*dr+br*di
      t = dp*h-dq*k
      dq = dp*k+dq*h
      dp = t
      if (pl.gt.46000.) go to 110
      if (abs(dp)+abs(dq).ge.(abs(p)+abs(q))*acc) go to 60
      p = p/r
      q = q/r
c *** solve for fp,g,gp and normalise f  at l=minl
      g = (fp-p*f)/q
      gp = p*g-q*f
      w = 1.0/sqrt(fp*g-f*gp)
      g = w*g
      gp = w*gp
      if (ktr.eq.1) go to 80
      f = tf
      fp = tfp
      lmax = maxl
c *** runge-kutta integration of g(minl) and gp(minl) inwards from turn
c ***             see fox and mayers 1968 pg 202
      if (rho.lt.0.2*turn) pace = 999.0
      r3 = 1.0/3.0d0
      h = (rho-turn)/(pace+1.0)
      h2 = 0.5*h
      rrpace = pace+0.001
      i2 = rrpace
      etah = eta*h
      h2ll = h2*xll1
      s = (etah+h2ll/r)/r-h2
 70   rh2 = r+h2
      t = (etah+h2ll/rh2)/rh2-h2
      k1 = h2*gp
      m1 = s*g
      k2 = h2*(gp+m1)
      m2 = t*(g+m1)
      k3 = h*(gp+m2)
      m3 = t*(g+k2)
      m3 = m3+m3
      k4 = h2*(gp+m3)
      rh = r+h
      s = (etah+h2ll/rh)/rh-h2
      m4 = s*(g+k3)
      g = g+(k1+k2+k2+k3+k4)*r3
      gp = gp+(m1+m2+m2+m3+m4)*r3
      r = rh
      i2 = i2-1
      if (abs(gp).gt.1.0e+36) go to 110
c     nb changed exponent to 36 for vax
      if (i2.ge.0) go to 70
      w = 1.0/(fp*g-f*gp)
c *** upward recursion from gc(minl) and gcp(minl),stored values are r,s
c *** renormalise fc,fcp for each l-value
 80   gc(lmin1) = g
      gcp(lmin1) = gp
      if (lmax.eq.minl) go to 100
      do 90 l=lmin1,lmax
         t = gc(l+1)
         gc(l+1) = (gc(l)*gc(l+1)-gcp(l))/gcp(l+1)
         gcp(l+1) = gc(l)*gcp(l+1)-gc(l+1)*t
         fc(l+1) = w*fc(l+1)
         fcp(l+1) = w*fcp(l+1)
 90   continue
      fc(lmin1) = fc(lmin1)*w
      fcp(lmin1) = fcp(lmin1)*w
      return
 100  fc(lmin1) = w*f
      fcp(lmin1) = w*fp
      return
 110  w = 0.0
      g = 0.0
      gp = 0.0
      go to 80
      end