Subroutine rcwfn (rho,eta,minl,maxl,fc,fcp,gc,gcp,accur,step) c *** see RH 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 GoTo 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) GoTo 30 r = r+1.0e-6 GoTo 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.) GoTo 110 If (abs(del/fp).ge.acc) GoTo 30 fp = f*fp If (lmax.eq.minl) GoTo 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) GoTo 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.) GoTo 110 If (abs(dp)+abs(dq).ge.(abs(p)+abs(q))*acc) GoTo 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) GoTo 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) GoTo 110 c nb changed exponent to 36 for vax If (i2.ge.0) GoTo 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) GoTo 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 GoTo 80 End