c/* %W%      latest revision %G% %U% */
      SUBROUTINE RHOK (nfr)
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 k,kp,Mp,MN
      REAL*8 RPT(96),WTR(96)
      dimension drho(4,100), y(100), ff(4), nifty(20)
      DIMENSION FJKP(100), FJK(100), PL(100), WT(4,100), RHOR(4),
     1 RNORM(2)
c     EQUIVALENCE (PL(1),FJK(1)), (WTR(1),DRHO(1,1))
      common /switch/   nifty
      common /sizes/ achp,acmp,wsp,achn,acmn,wsn
      common /woods/ drho,kp,k,imax
      common /params/   hbarc , pi , Mp , MN , nZ , nA , nes , nwaves
      common /inputs/   Tlab , b , y1 , y2 , xnang ,  kode , lmax ,
     $                  ngp , nr
C
      DATA FJK/100*0. /,NDATA/0/,ACHECK/0./,RNORM/2*0. /
      ANEW = ACHP+ACHN+ACMN+ACMP+WSP+WSN
      NFF=nfr
      if(nff.gt.2)nff=2
C     IF ((NDATA.NE.0).AND.(ACHECK.EQ.(ANEW))) GO TO 30
      IF (ACHECK.EQ.ANEW) GO TO 30
C     FIRTST CALL TO RHOK,SET UP ROINTS,DETERMINE NORMALIZATION
      NRS = 96
      R = ACHP/hbarc
      R2 = 5.*R
      R1 = 1.25*R*R2/(R2-1.25*R)
      CALL GAUSS2 (NRS,1,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.*pi
            RNORM(J) = RNORM(J)+WT(J,NR)
 10      CONTINUE
 20   CONTINUE
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
c
c normalize
c
      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
