      SUBROUTINE RHOK (NFF)
C *** SEE R H 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
      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))) GO TO 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) GO TO 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
