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

      subroutine ffact (q2,ff,nff)
c *** see r h 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) go to 70
      if (na.eq.3) go to 80
      if (na.eq.2) go to 100
      if (na.lt.20) go to 120
      if ((ndata.ne.0).and.(acheck.eq.(anew))) go to 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.) go to 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) go to 90
      if (nifty(19) .eq. 1) go to 85
      if (nifty(19) .eq. 2) go to 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) go to 81
      if(nifty(16).eq.5.or.nifty(16).eq.6) go to 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)
      go to 83
c     nifty(16)=5..fnsp=fnmat=fpm   ,fpsp=0
  82  ff(2)=ff(1)
      go to 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
       go to 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.) go to 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