Subroutine rhok (nff)
c *** see RH Landau, Program lpott, 1981
c     given a coordinate (r ) space density
c     sub rhok numerically calculates the partial wave amplitudes of the
c     form factor for a k + kp, and stortes them in Common/sec5   for
c     imax partial waves  and nff f.f. s
c     n.b. this Doesnot Do spin distrbs, for he3 know ff*s so ok
c     given a coordinate saspace(r)  dedensity
c     n.b. r is in mev-1
      Implicit Real*8   (a-h,o-z)
      Real*8 rpt(48),wtr(48),r1,r2
      Real*8 kp,k,achp,achn,acmp,acmn,rho,wsp,wsn,r,rhor,acheck
      Dimension rho(4,50), fjkp(60), fjk(60), pl(50), wt(2,50), rhor(4),
     1 rnorm(4)
      Common /sec4/ achp,acmp,wsp,achn,acmn,wsn
      Common /sec5/ drho(4,50),rho,kp,k,imax
ccc   Equivalence (pl(1),fjk(1)), (wtr(1),drho(1,1))
c
      data fjk/60*0. /,ndata/0/,acheck/0./,rnorm/4*0. /
      anew = achp+achn+acmn+acmp+wsp+wsn
      If ((ndata.ne.0).and.(acheck.eq.(anew))) GoTo 30
c     firtst call to rhok,set up roints,determine normalization
      nrs = 24
      r = achp/197.32891
      r2 = 5.*r
      r1 = 1.25*r*r2/(r2-1.25*r)
      Call Gauss (nrs,052,r1,r2,rpt,wtr)
      ndata = 1
      acheck = anew
      Do 20 nr=1,nrs
         r = rpt(nr)
         Call rhoofr (r,rhor,nff)
         Do 10 j=1,nff
            wt(j,nr) = wtr(nr)*r*r*rhor(j)*4.*3.1415927
            rnorm(j) = rnorm(j)+wt(j,nr)
 10      Continue
 20   Continue
      rnorm(3) = rnorm(1)**2
      rnorm(4) = rnorm(2)**2
c     r integration
 30   Do 50 i=1,imax
         Do 40 j=1,nff
            drho(j,i) = 0.
 40      Continue
 50   Continue
c     ---- -major Do loop over r
      Do 90 nr=1,nrs
         r = rpt(nr)
         r1 = kp*r
         Call spbesl (imax,r1,fjkp)
         If (kp.eq.k) GoTo 60
         r1 = k*r
         Call spbesl (imax,r1,fjk)
 60      Do 80 i=1,imax
            rhl1 = fjkp(i)
            rhl2 = fjk(i)
            If (kp.eq.k) rhl2 = rhl1
            Do 70 j=1,nff
               drho(j,i) = drho(j,i)+rhl1*rhl2*wt(j,nr)
 70         Continue
 80      Continue
 90   Continue
      Do 110 i=1,imax
         Do 100 j=1,nff
            drho(j,i) = drho(j,i)*(2.*i-1.)/rnorm(j)
 100     Continue
 110  Continue
      Return
      End