
c/* %W%      latest revision %G% %U% */
      subroutine FFMMM( qf2 , f1c , f2c , psp , pd )
************************************************************************

c *** See R.H. Landau's program LPOTT, 1981
c     fmmm interpolates Malcolm McMillan form factors and forms f1c and
c     f2f from them (based on approx triton wave functns)
c     psp = s prime prob     pd= dstate prob     qf2= q2 in ferm-2

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

      dimension q2(18), f3456(4,18), y(4), q20(27), f0(1,27), y0(1)
      dimension a1(18),a2(18),a3(18),a4(18)

      data q2  / 0.0d+0,  0.1d+0,  1.1d+0,  2.1d+0,  3.1d+0,  4.1d+0,
     $           5.1d+0,  6.1d+0,  7.1d+0,  8.1d+0,  9.1d+0, 10.1d+0,
     $          11.1d+0, 12.1d+0, 13.1d+0, 14.1d+0, 15.1d+0, 16.1d+0  /

      data q20 / 0.d+0,   0.3d+0,  0.9d+0,  1.5d+0,  2.1d+0,  2.7d+0,
     $           3.3d+0,  3.9d+0,  4.5d+0,  5.1d+0,  5.7d+0,  6.3d+0,
     $           6.9d+0,  7.5d+0,  8.0d+0,  9.0d+0, 10.0d+0, 11.0d+0,
     $          12.0d+0, 13.0d+0, 14.0d+0, 15.0d+0, 16.0d+0, 17.0d+0,
     $          18.0d+0, 19.0d+0, 20.0d+0                             /

      data f0  / 1.0d+0,       0.881240d+0,  0.693630d+0,  0.553330d+0,
     $           0.445690d+0,  0.361530d+0,  0.294750d+0,  0.241140d+0,
     $           0.197700d+0,  0.162210d+0,  0.133050d+0,  0.108960d+0,
     $           0.088970d+0,  0.072320d+0,  0.605673d-1,  0.417046d-1,
     $           0.276478d-1,  0.171613d-1,  0.935277d-2,  0.356770d-2,
     $          -0.679445d-3, -0.375343d-2, -0.592979d-2, -0.741903d-2,
     $          -0.838265d-2, -0.894523d-2, -0.920333d-2              /

c *** The following data statements were added by M.S. to
c *** eliminate need for block data.

      data a1  / 0.0d+0,        0.294322d-01,  0.208823d+00,
     $           0.274661d+00,  0.291498d+00,  0.285677d+00,
     $           0.269357d+00,  0.248563d+00,  0.226397d+00,
     $           0.204474d+00,  0.183614d+00,  0.164207d+00,
     $           0.146402d+00,  0.130215d+00,  0.115587d+00,
     $           0.102426d+00,  0.906173d-01,  0.800462d-01           /

      data a2  / 2.0d+00,       0.191074d+01,  0.127537d+01,
     $           0.894039d+00,  0.647425d+00,  0.480273d+00,
     $           0.363230d+00,  0.279278d+00,  0.217938d+00,
     $           0.172452d+00,  0.138314d+00,  0.112430d+00,
     $           0.926325d-01,  0.773725d-01,  0.655261d-01,
     $           0.562674d-01,  0.489832d-01,  0.432135d-01           /

      data a3  / 0.526184d-01,  0.505840d-01,  0.341021d-01,
     $           0.236391d-01,  0.166030d-01,  0.116959d-01,
     $           0.818813d-02,  0.563686d-02,  0.375893d-02,
     $           0.236579d-02,  0.132789d-02,  0.553974d-03,
     $          -0.215559d-04, -0.446686d-03, -0.757024d-03,
     $          -0.979299d-03, -0.113380d-02, -0.123603d-02           /

      data a4  / 0.0d+00,      -0.256572d-03, -0.182162d-02,
     $          -0.232969d-02, -0.237360d-02, -0.221288d-02,
     $          -0.196895d-02, -0.170076d-02, -0.143704d-02,
     $          -0.119130d-02, -0.969153d-03, -0.772119d-03,
     $          -0.599600d-03, -0.449968d-03, -0.321150d-03,
     $          -0.210971d-03, -0.117277d-03, -0.380834d-04           /

************************************************************************

c *** The following statements were added by M.S. to make
c *** the above data equivalent to the previous block data

      do 900 j=1,18
         f3456(1,j) = a1(j)
         f3456(2,j) = a2(j)
         f3456(3,j) = a3(j)
         f3456(4,j) = a4(j)
 900  continue

      if (qf2.le.16.1) go to 20
      q = sqrt(qf2)
      if ((0.4111*qf2.lt.150.).and.(1.91*(q-4.4721)**2.lt.150.))
     1 go to 10
c     q2 too large for exp, set=0
      f1c = 0.
      f2c = 0.
      return
c     q2 outside of table, use annaltic forms fitted to last oints
 10   y0(1) = exp(-0.4111*qf2)-0.94722e-2*exp(-1.906*(q-4.4721)**2)
      y0(1) = y0(1)*1.23469
      y(1) = 0.589748*exp(-0.124043*qf2)
      y(2) = 0.325014*exp(-0.125324*qf2)
      y(3) = (1./19.0048)*(exp(-0.4667*qf2)-0.25369e-01*exp(-1.94305*(q-
     14.1833)**2))
      y(4) = -0.831500*exp(-0.587182*qf2)
      go to 30
 20   continue
      if(qf2.lt.0.01)go to 40
c     interpolate on form factors,,n.b. may have to cahge this for more
c     accuracy and interpolate on logs of ff and the expeonentiate
      call lagrng (qf2,q2,y,f3456,18,4,2,18,4)
      call lagrng (qf2,q20,y0,f0,27,1,2,27,1)
 30   ps = 1.-psp-pd
      f1c = ps*y0(1)+0.5*psp*y(2)+19.0048*pd*y(3)
      f2c = sqrt(psp*ps)*y(1)+19.0048*pd*y(4)
c     write(6,900) qf2,y0(1),y(1),y(2),y(3),y(4),f1c,f2c
c     900 format(" ", f10.3,7e15.5)
      return
c     q2 lt 1st point, set f= 1,0
 40   f1c=1.
      f2c=0.
      return
c**  this program valid on ftn4 and ftn5 **
      end
