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

c  this subroutine calculates the matter and spin form factors for proton*
c  and neutrons in 13C using Ray's tables of core form factors and rhl's  *
c  fit to Ray's valence n's psi(r), and it's analytic Fourier transform   *
c                                                                         *
c  f1,2,3,4  : pmatter, nmatter, pspin, nspin
c                                                                         *
c *************************************************************************

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

      dimension formfs(2,21),y(2),q2ray(21)
      dimension fcorep(21),fcoren(21)
      dimension a2(4),ff(4)

      common /params/    hbarc, pi, Mp, MN, nz, na, nes, nwaves

      data q2ray  / 0.0000d+0,  0.0625d+0,  0.2500d+0,   0.5630d+0,
     $              1.0000d+0,  1.5630d+0,  2.2500d+0,   3.0630d+0,
     $              3.4230d+0,  4.0000d+0,  5.0630d+0,   6.2500d+0,
     $              7.5600d+0,  9.0000d+0,  1.0563d+1,   1.2250d+1,
     $              1.4060d+1,  1.6000d+1,   1.6500d+1,  1.7000d+1,
     $              1.7600d+1  /

      data a2     / 0.9820d+0,  0.9820d+0,  0.7592d+0,   0.7592d+0/

      data fcorep  / 1.000d+00,   9.450d-01,  7.962d-01,   5.950d-01,
     $               3.893d-01,   2.157d-01,  9.175d-02,   1.778d-02,
     $               6.160d-05,  -1.656d-02, -2.538d-02,  -2.165d-02,
     $              -1.454d-02,  -8.613d-03, -4.867d-03,  -2.637d-03,
     $              -1.197d-03,  -2.700d-04, -1.833d-04,  -9.667d-05,
     $              -5.330d-05/

      data fcoren  / 1.000d+00,   9.486d-01,  8.091d-01,   6.166d-01,
     $               4.147d-01,   2.391d-01,  1.089d-01,   2.724d-02,
     $               6.592d-03,  -1.387d-02, -2.712d-02,  -2.530d-02,
     $              -1.816d-02,  -1.112d-02, -5.980d-03,  -2.659d-03,
     $              -5.377d-04,   5.831d-04,  6.842d-04,   5.841d-04,
     $               4.005d-04/


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

       q2f2 = q2/(hbarc * hbarc)

c *** if q2f2 outside of table then extrapolate with ***
c ***
c ***          fcore(q2) = fcore(q02)exp -(q2-q02)/a2        ***
c ***
c *** where f(q02) = last point in table
c *** (ie. q02=17.6 fm-2)

C     store form factor as 2D array for interpolation
C
      do 70 j=1,21
         formfs(1,j) = fcorep(j)
         formfs(2,j) = fcoren(j)
  70    continue
      if (q2f2 .lt. 17.6) go to 40
         lpt = 17.6
c        beyond table, extrapolate, set large q2a2 to 150
         do 20 i=1,2
      q2a2 = (q2f2-lpt)/a2(i)
      if (q2a2  .gt. 150.) q2a2 = 150.
      ff(i) = formfs(i,21)*exp(-q2a2)
  20       continue
         ff(3) = 0.
         ff(4) = 0.
         go to 30
c     within tabulated region, interpolate with 2 or 4 points
  40  continue
      if ((q2f2 .ge. 17.0) .or. (q2f2 .lt. 0.625)) go to 50
c
c
c*************************************************************************
c                                                                        *
c  q2f2 lt 17.0, or gt 0.625 use 4 pt interpolation                      *
c                                                                        *
c*************************************************************************
c
c
      call lagrng(q2f2,q2ray,y,formfs,21,2,4,21,2)
      ff(1) = y(1)
      ff(2) = y(2)
      ff(3) = 0.
      ff(4) = 0.
c    add on n valence piece "fval", norm 1, then rescale
 30   continue
      aa = a2(4)
      rhl = 4. * aa + q2f2
      rhl2 = rhl*rhl
      rhl3 = rhl*rhl2
      a4= aa**2
      a6= a4*aa
      fval = (12.8*a6/rhl3)*(1.-48.*aa/rhl + 256.*a4/rhl2)
c   for c13 take neutron # = N = 7
c   now add valence piece
      n=7
      ff(2) = ((n-1)*ff(2) + fval)/n
      if(nff.gt.2)ff(4) = fval/n/3.
      return

  50  continue
c
c
c**************************************************************************
c                                                                         *
c  q2f2 lt 17.6.25 but gt than 17.0, use 2-point                          *
c  interpolation since in between the last two data points  of table      *
c                                                                         *
c**************************************************************************
c
c
      call lagrng(q2f2,q2ray,y,formfs,21,2,2,21,2)
      ff(1) = y(1)
      ff(2) = y(2)
      ff(3) = 0
      ff(4) = 0.
      go to 30
      end