      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*4 RRPACE
      DIMENSION FC(1), FCP(1), GC(1), GCP(1)
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 = FLOAT(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 = FLOAT(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 = IFIX(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
