c/* %W%      latest revision %G% %U% */
      subroutine FFACT (q2, ff, nff)
************************************************************************

c *** See R.H. Landau's program LPOTT, 1981
c *** Calculates the form factor
c *** via special direct formular

c *** Can now handle rhop .ne. rhon and spin,
c *** nff = f.f.*s needed, ff(i)= rho p,n,psp,nsp for i=1,2,3,4

c *** REVISION HISTORY
c
c        2/2/87      add form factor for Carbon-13    (ms)
c       30 Nov 88     add Ray's form factor for Carbon-13 (NIFTY(16)=8)
c        6 Dec 88    add constant form factor (=1)

      implicit real*8 (a-h, o-z)
      real*8 Mp, MN

      dimension nifty(20), ff(4)

c
      common /params/   hbarc, pi, Mp, MN, nz, na, nes, nwaves
      common /switch/   nifty
      common /sizes/    achp, acmp, wsp, achn, acmn, wsn

c ***

      data ndata/0/

c >>> FIRST EXECUTABLE STATEMENT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

c
c FORM FACTORS CURRENTLY AVAILABLE FOR A=2, 3, 4, 12, 13, 16 + WS for A>20
c
c for A=13 Nifty(16)=8 for Ray's form factor PRC37 (1988) p1169
c          Nifty(16)=1,2 Woods-Saxon form factor
c          Nifty(16)=0 etc. Harmonic Oscillator form factor
c
c for all  Nifty(16)=9 gives constant form factor (=1)

      if(Nifty(16).eq.9)then
      ff(1)=0.
      ff(2)=0.
      ff(3)=0.
      ff(4)=0.
      if(nff.ge.1)ff(1)=1.
      if(nff.ge.2)ff(2)=1.
      if(nff.ge.3)ff(3)=1.
      if(nff.ge.4)ff(4)=1.
      return
      endif

      if (na .eq. 4) then
c     *** Special form factor for He4
c     ***  ach = beta*2       acm = achp

         h2 = hbarc * hbarc

c        *** Proton distribution function
c        *** Set spin ff = 0 for He4

         ff(3) = 0.
         ff(4) = 0.
         betap = achp/2.
         barnie = 4.0e-06
         aarnie = -0.14
         ff(1) = 0.
         if (q2 .lt. 8.0e05) then
            rhl = q2 * betap * betap/h2
            ff(1) = ( 1. - (acmp * acmp * q2/h2)**6 ) * exp(-rhl)
         else
            rhl = q2 * barnie
            if (rhl .gt. 150.) then
               ff(1) = 0.
            else
               ff(1) = aarnie * exp(-rhl)
            end if
         end if

c        *** Divide out dipole fit to proton form factor

         if (nifty(15) .eq. 0) then
            ff(1) = ff(1) * ( 1. + q2/18.2/h2 )**2
         endif
         if (nff .eq. 1) return

c        *** Neutron form factor.

         betan = achn/2.
         ff(2) = 0.
         if (q2 .lt. 8.0e05) then
            rhl = q2 * betan * betan/h2
            ff(2) = ( 1. - (acmn * acmn * q2/h2)**6 ) * exp(-rhl)
         else
            rhl = q2 * barnie
            if (rhl .gt. 150.) then
               ff(2) = 0.
            else
               ff(2) = aarnie * exp(-rhl)
            end if
         end if

c        *** Divide out proton form factor

         if (nifty(15) .eq. 0) then
            ff(2) = ff(2) * ( 1. + q2/18.2/h2 )**2
         endif
c        if (nff .gt. 2) write(6,130) nff

      else if (na .eq. 3) then

c        *** He3 nucleus ***

         nfcn = nifty(15)
         if (nifty(16) .eq. 0) go to 90
         if (nifty(19) .eq. 1) go to 85
         if (nifty(19) .eq. 2) go to 85

c        *** Special case for He3 (McCarthy et al PRL 25, 884(1970))
c        *** ach = 2a   acm = 2c, other sizes built in, but may want to
c        *** change

         call FFHE3( q2,ff(1),ff(2),ff(3),ff(4),nfcn,ache,amhe,ach)

         if (nifty(16) .eq. 4 .or. nifty(16) .eq. 7) go to 81
         if (nifty(16) .eq. 5 .or. nifty(16) .eq. 6) go to 82

c        *** Special H3 case

 83      if (nz .eq. 2) RETURN
         rhl = ff(1)
         ff(1) = ff(2)
         ff(2) = rhl
         rhl = ff(3)
         ff(3) = ff(4)
         ff(4) = rhl
         RETURN

c        *** Nifty(16) = 4  ..fnsp = fnmat, fpsp = 0

c81      continue
 81      ff(3) = 0.
         ff(4) = ff(2)
         go to 83

c        *** Nifty(16) =5     fnsp = fnmat = fpmat, fpsp = 0

  82     ff(2) = ff(1)
         go to 81

  85     n19 = nifty(19)

         call FFHADJ (q2, ff1, ff2, ff3, ff4, n19)
         ff(1) = ff1
         ff(2) = ff2
         ff(3) = ff3
         ff(4) = ff4
         go to 83

c        *** Use Gaussian form factors for He3 and H3 charge and He3 mag
c        *** achp---He3 charge      achn---H3 ch          acmp---He3 mag

 90      ache = achp
         amhe = acmp
         ach = achn
         nfcn = 0

         call FHE3GS( q2,ff(1),ff(2),ff(3),ff(4),nfcn,ache,amhe,ach)

      else if (na .eq. 2) then

c        *** Deuteron

         a = achp * hbarc
         b = acmp * hbarc
         anorm = 1./(0.5/a + 0.5/b - 2./(a+b))
         if (q2 .eq. 0.) then
            ff(1) = 1.
            ff(2) = 1.
            ff(3) = 1.
            ff(4) = 1.
         else
            q = sqrt(q2)
            ff(1) = anorm/q * ( atan(q/(2.*a)) + atan(q/(2.*b))
     $                          - 2. * atan(q/(a+b)) )
            ff(2) = ff(1)
            ff(3) = ff(1)
            ff(4) = ff(1)
         endif
         RETURN

      else if (na .eq. 13 ) then

c        *** Special form for Carbon-13 ***

c for Ray's form factor
         if(NIFTY(16).eq.8)then
         call fray13( q2, nff, ff )
c for W-S form factor
         else if(NIFTY(16).eq.1.or.NIFTY(16).eq.2)then
         call fwoods( q2, nff, ff )
         else
c for standard HO form factor
         call  FFC13( q2, nff, ff )
         endif

      else if (na .lt. 20) then

c        *** Special form for harmonic oscillator nuclei

         call FFHMO (q2, nff, ff(1), ff(2))

      else

c        *** Special form for Woods-Saxon density (na>20)

         call fwoods (q2, nff, ff)


      endif

c *** FORMATS ***

 130  format (29h xxxxbuggy in ffact nff gt2,=,i5)

c ***
      end