      Program facform
c***  see RH Landau Program Lpott 1981
c     Computes He3 matter and spin Form Factors for p-matter
c     and n-matter.
c     Calls routine ffhe3
c     produces files ropmat.dat
c 
      Implicit Real*8(a-h,o-z)    
      Common /sec4/achp,acmp,wsp,achn,acmn,wsn          
      achp=0.6540
      acmp=0.4560
      wsp=0.810
      achn=0.6750
      acmn=0.366
      wsn=0.8360
      Open(6,File='ropmat.dat',Status='Unknown')
c     q2 given in MeV**2 ffh3 converts to fm**-2
      q2mev=0.0
      nfcn=0
      Do 10 i=1,00
        Call ffhe3 (q2mev,rp,rn,rpsp,rnsp,nfcn,achp,acmp,achn)
        q2fermis2=q2mev/(197.329**2)
        Write(*,*)q2fermis2, rp
        q2mev=q2mev+5000.0
 10   Continue
      Close(6)
      Stop 'Data are in ropmat.dat'
      End

      Subroutine ffpn (q2,fcp,fcn)
c *** see RH Landau, Program lpott, 1981
c     proton charge ff (also p mag and n mag If nrom(0)=1)
c     neutron charge form factor
c     q in f-1
      Implicit Real*8 (a-h,m,o-z)
c
      data x,un,mn/18.2,-1.91315,4.7548/
      fcp = 1./((1.+q2/x)**2)
      fcn = -(un*q2*0.25/mn/mn)*fcp
      Return
      End


      function ffche3 (q)
      Implicit Real*8(a-h,o-z)
c *** see RH Landau, Program lpott, 1981
c     charge form factor for he3 from mccarthy et al, prl25,884,1970
c     q in f- 1
c           /60.1,70.4,3-e59.5-,24357.,88251.,984664./p,0q,d,2c,2b,2a at
c          /147.,98.3,3-e16.7-,014646.,82611.,988444./p,0q,d,2c,2b,2a at
c
      Common/sec4/a2m,b2m,c2m,a2,b2,c2
      data d,q0,p/ -6.78e-3,3.98,0.90/
cccc  data a2,b2,c2,d,q0,p/.455625,.13396,.698896, -6.78e-3,3.98,0.90/
      q2 = q*q
      rhl1 = 0.
     
      If(q2.gt.24.)GoTo 40
      rhl = a2*q2
      If (rhl.ge.150.) GoTo 10
      rhl1 = exp(-rhl)
 10   rhl2 = 0.
      rhl = c2*q2
      If (rhl.ge.150.) GoTo 20
      rhl2 = -b2*q2*exp(-rhl)
 20   rhl = ((q-q0)/p)**2
      rhl3 = 0.
      If (rhl.ge.150.) GoTo 30
      rhl3 = d*exp(-rhl)
 30   ffche3 = rhl1+rhl2+rhl3
      Return
c     arnold et al fch(he30 for largest q
 40   ffche3=0.
      rhl=0.1059*q2
      If(rhl.ge.150.)Return
      tau=q2/4./202.5
      rhl1=Sqrt((1.+tau)/(1.+4.527*tau))
      ffche3=rhl1*0.034*exp(-rhl)
      Return
c**  this Program valid on ftn4 and ftn5 **
      End

      function ffmhe3 (q)
      Implicit Real*8(a-h,o-z)
c *** see RH Landau, Program lpott, 1981; mods by m sagen
c     magnetuc ff of he3 from mcarcthy et al
c     q in f-1
c     cccc data a2,b2,c2/0.4277 ,0.207936,0.674041/
      Common /sec4/ a2,b2,c2,achn,acmn,wsn
c                                 /16727.0,522532.0,486954.0/2c,2b,2a at
c                                   /125226.0,33281.0,9693.0/2c,2b,2a at
      q2 = q*q
      rhl1 = 0.
      rhl = a2*q2
      If (rhl.ge.150.) GoTo 10
      rhl1 = exp(-rhl)
 10   rhl = c2*q2
       rhl2 = 0.
      If (rhl.ge.150.) GoTo 20
      rhl2 = -b2*q2*exp(-rhl)
 20   ffmhe3 = rhl1+rhl2
      Return
c**  this Program valid on ftn4 and ftn5 **
      End


      Subroutine ffhe3 (q2mev,rp,rn,rpsp,rnsp,nfcn,ache,amhe,ach)
c *** see RH Landau, Program lpott, 1981, mods by m sagen
c     m paez eliminated use of fmmm- McMillan 3-nucleon ff
c     ffhe3 calculates the four form factors needed to describe he3with
c     (fsd is only calc as a cpmparison, it should be small and is
c     the apptox, or effectove, sd contib to the magnetic ff
c     the ff*s are calc by using
c       1 fche3- expt of 7ccarthy et al for all q2,have fitted
c       2 ff*s calc by m McMillan,for fcn=0 these fit h3 data
c         of collard and for larger q2( gt 6 f-2) they rep a
c         reasonable extrapltn which has a correct zero as he3
c     3 fmhe3 use the form of 7ccarthy et al but vary this since the exh
c     chage contib is included in e scatter, but presmbly not in our pi
c     scatt (nb we use mu he3 = un)
c     these ff include s,s* + d states exactly in f1c,f2c,
c     the d state is only put in approximatly in f1m,f2m as Gibson has
c     shown this to be small (ignore too the dd terms too)
c     nfcn=1 , include neutron form factor in escat data, =0 , Don*t
c     data seems best interpreted with no neutron ?charge? form factor
c     possibly this is due to polarization or meson exch efects &&
c
c     q2mev2 is in mev**2 not Fermi-2
      Implicit Real*8(a-h,o-z)
      data up,un,uhe3/2.79278,-1.91315,-1.91315/
      q2f2 = q2mev/(197.329**2)
      qf = Sqrt(q2f2)
c     nucleon ff
      fn=0.
      fp=1.
      If(nfcn.eq.0) Call ffpn (q2f2,fp,fn)
      If ((nfcn.ne.1).and.(nfcn.ne.0)) Write (6,40) nfcn
c     he3 and h3 charge ff
      fche3 = ffche3(qf)
c     use mmm wf*s for h3,this fits collard (if fn=0) and is reasnb,e
c     sp= s* prob, pd= d state prob, these a can be varied
      psp = 0.02
      pd = 0.09
c      Call fmmm (q2f2,f1c,f2c,psp,pd)
c      fch3 = (2.*fn+fp)*f1c+(fp-fn)*f2c
c     he3 mag ff
      fmhe3 = ffmhe3(qf)
      GoTo 30
      entry fhe3gs(q2mev,rp,rn,rpsp,rnsp,nfcn,ache,amhe,ach)
      q2f2 = q2mev/(197.32**2)
      fn = 0.
      fp = 1.
c     Gaussian form factors
      fche3 = 0.
      rhl = ache*ache*q2f2
      If (rhl.gt.150.) GoTo 10
      fche3 = exp(-rhl)
 10   fch3 = 0.
      rhl = ach*ach*q2f2
      If (rhl.ge.150.) GoTo 20
      fch3 = exp(-rhl)
 20   fmhe3 = 0.
      rhl = amhe*amhe*q2f2
      If (rhl.ge.150.) GoTo 30
      fmhe3 = exp(-rhl)
 30   Continue

c
c     direct calc of ff for fcn=0
      rp=fche3/fp
     
      rn=fch3/fp
      rhl=un/(up+2.*un)/fp
      rhl1=(fch3-4.*fche3)/3.
      rpsp=rhl*(fmhe3+rhl1)/2.
      rnsp=rhl*(2.*fmhe3-rhl1*(up/un))
      Return
c
 40   format (55h ********** neutron form factor switch is strange,nfcn=
     1,i3,11h **********,/)
      End





























