c/* %W%      latest revision %G% %U% */
      Subroutine FFC13( q2, nff, ff )
c***********************************************************************

c *** Calculates the form factors the Carbon-13
c ***
c *** Based on formula given by Landau in LPOTT paper (CPC)

c *** Written by Milt Sagen February 2, 1988

c***********************************************************************

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

      dimension ff(4)
      dimension nifty(20)

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

c***********************************************************************

      nn = na - nz

c     *** convert q2 from MeVs to inverse fermis

      qfm2 = q2/(hbarc * hbarc)

c     *** use the harmonic oscillator model for the proton ***

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

c     *** set spin form factor of proton to zero

      ff(3) = 0.

c     *** now do special form for neutron ***
c
c     *** matter is Gaussian for N-1 neutrons ***

      factor = 1.- (nn-3.)/(6*(nn-1.)) * (qfm2*acmn**2)
      factor = (nn-1.)/nn * factor
      arg=-qfm2*achn**2/4.0
      if(arg.gt.-150.0)then
c     ffgn = factor * exp(-qfm2*achn**2/4.)
      ffgn = factor * exp(arg)
      else
      ffgn=0.0
      endif

c     *** and now the extra term ***

      factor = (1. - (qfm2 * wsn**2)/6.)/nn
      arg=-qfm2*wsn**2/4.0
      if(arg.gt.-150.0)then
c     ffadd = factor * exp(-qfm2 * wsn**2/4.)
      ffadd = factor * exp(arg)
      else
      ffadd=0.0
      endif
      ff(2) = ffgn + ffadd

c     *** now do neutron spin form factor ***

      ff(4) = ffadd/3.

      Return
      End