c @(#)ffact	1.1	2/4/86 12:22:07

      Subroutine ffact (q2,ff,nff)
c *** see RH Landau, Program lpott, 1981; m sagen mods
c     ffact calc form factor via profjection of partial waves from 1 exs
c      or via special direct formular If ff is known analytically
c     can now handle rhop .ne. rhon and spin,
c     nff= = f.f.*s needeed, ff(i)= rho p,n,psp,nsp for i=1,2,3,4
      Implicit Real*8(a-h,o-z)
      Real*8 imtij,mpi,mn,k,kp
      Dimension rho(4,50), y(50), ff(4)
      Dimension nifty(20), bnuc(20), bnucf(20), bn(16), bnf(16), retij
     1(14),imtij(14)
      Common /sec2/ bnuc,bn,bnucf,bnf,retij,imtij,hbarc,pi,mpi,mn,nz,nes
     1,nwaves,nifty,na
      Common /sec4/ achp,acmp,wsp,achn,acmn,wsn
      Common /sec5/ drho(4,50),rho,kp,k,imax
c
      data ndata/0/,acheck/0./
      anew = achp+acmp+wsp+achn+acmn+wsn
      If (na.eq.4) GoTo 70
      If (na.eq.3) GoTo 80
      If (na.eq.2) GoTo 100
      If (na.lt.20) GoTo 120
      If ((ndata.ne.0).and.(acheck.eq.(anew))) GoTo 10
c     1st call to ffact,or nuclear size is changed
c      use numerical expns, determine expns coefs for a very large q
      ndata = 1
      acheck = anew
      y(1) = kp
      y(2) = k
      kp = 345.31
      k = 345.31
c     maybe need imax only u/-12 yet jmax=20
      jmax = 20
      If (jmax.lt.imax) jmax = imax+2
      iimax = imax
      imax = jmax
c     Call rhok; sets up partial wave cof in sec5*s rho, for all ff*s If
      Call rhok (nff)
      kp = y(1)
      k = y(2)
      imax = iimax
c     check size of q
c      modifeied so no switch form for one kp,k set
 10   If ((kp+k).lt.700.) GoTo 30
c     q too large for ff exsnp to be valid, use a cut off(unphysical reg
      rhl = 0.
      xx = q2/5.e04
      If (na.gt.100) xx = q2/2.8e04
      If (xx.lt.150.) rhl = exp(-xx)
      Do 20 j=1,nff
         ff(j) = rhl
c     scale for rho**2
         If (j.ge.3) ff(j) = rhl*3.*((197.32/1.2)**3)/4./pi/na
 20   Continue
      Return
c     calc series for f(q) , first deteermine costheta, then use legrend
 30   xx = 1.-q2/23.8478e04
      y(1) = 1.
      y(2) = xx
      Do 40 j=1,nff
         ff(j) = drho(j,1)+xx*drho(j,2)
 40   Continue
      ii = jmax-1
      Do 60 i=2,ii
         g = xx*y(i)
         y(i+1) = 2.*g-y(i-1)-(g-y(i-1))/i
         Do 50 j=1,nff
            ff(j) = ff(j)+drho(j,i+1)*y(i+1)
 50      Continue
 60   Continue
      Return
c     special form factor for he 4
c       ach=beta*2       acm= achp
 70   h2 = hbarc*hbarc
c     proton distrb functin
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 toproton form factor
      If(nifty(15).eq.0)ff(1) = ff(1)*(1.+q2/18.2/h2)**2
      If (nff.eq.1) Return
c     neutron f.f.
      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 ff
      If(nifty(15).eq.0)ff(2) = ff(2)*(1.+q2/18.2/h2)**2
      If (nff.gt.2) Write (6,130) nff
      Return
 80   nfcn = nifty(15)
      If (nifty(16).eq.0) GoTo 90
      If (nifty(19) .eq. 1) GoTo 85
      If (nifty(19) .eq. 2) GoTo 85
c     special case for he 3 (mc?arthyet al prl 25,884(1970))
c      ach=2a   acm=2c,  other sizes built in, but may want to change
c     Call ffhe3 (q2,ff(1),ff(2),ff(3),ff(4),nfcn)
c%%%%%%%%%%note provisional last comment%%%%%%%%%%%
      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) GoTo 81
      If(nifty(16).eq.5.or.nifty(16).eq.6) GoTo 82
 83   If (nz.eq.2) Return
c     special h3 case
      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=snmat     ,fpsp=0
 81   ff(3)=0.
      ff(4)=ff(2)
      GoTo 83
c     nifty(16)=5..fnsp=fnmat=fpm   ,fpsp=0
  82  ff(2)=ff(1)
      GoTo 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
       GoTo 83
c     use Gaussian form factors for he 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)
      Return
c     deuteron , hulthen Double exponential
 100  a = achp*hbarc
      b = acmp*hbarc
      anorm = 1./(0.5/a+0.5/b-2./(a+b))
      If (q2.ne.0.) GoTo 110
      ff(1) = 1.
      ff(2) = 1.
      ff(3) = 1.
      ff(4) = 1.
      Return
 110  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)
      Return
c     special form for harmonic oscillator nuclei
 120  Call harmnuc(q2,nff,ff(1),ff(2))
      Return
c
 130  format (29h xxxxbuggy in ffact nff gt2,=,i5)
c**  this Program valid on ftn4 and ftn5 **
      End