
c/* %W%      latest revision %G% %U% */
      subroutine fwoods (q2,nff,ff)
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 k,kp,kmax,Mp, MN
      dimension drho(4,100), y(100), ff(4), nifty(20)
      common /switch/   nifty
      common /sizes/ achp,acmp,wsp,achn,acmn,wsn
      common /woods/ drho,kp,k,imax
      common /params/   hbarc , pi , Mp , MN , nZ , nA , nes , nwaves
      common /inputs/   Tlab , b , y1 , y2 , xnang ,  kode , lmax ,
     $                  ngp , nr
c
      data ndata/0/,acheck/0./
      anew = achp+acmp+wsp+achn+acmn+wsn
         h2 = hbarc * hbarc
c     if ((ndata.ne.0).and.(acheck.eq.(anew))) go to 10
      if (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
      imax=lmax
      ndata = 1
      acheck = anew
c
c kmax corresponds to qmax 5.2invfm i.e. 60 degree. p C  547MeV
      kmax=512.652
c kmax corresponds to qmax 10.5invfm i.e. 180 degree. p C  547MeV
c     kmax=1031.0
      kp = kmax
      k = kmax
      q2max=(2*kmax)**2
c     maybe need imax only u/-12 yet jmax=20
      jmax = 20
      if (jmax.lt.imax) jmax = imax+2
      if (jmax.gt.100) jmax=100
      iimax = imax
      imax = jmax
c     call rhok; sets up partial wave cof in woods's rho, for all ff*s if
      call rhok (nff)
      imax = iimax
c     check size of q
c      modifeied so no switch form for one kp,k set
 10   if (q2.ge.q2max) then
c     q too large for ff exsnp to be valid, use a cut off(unphysical reg
      rhl = 0.
       if(nA.ge.40.and.nA.le.100)then
       xx = (q2-q2max)/5.e04
       else if(nA.eq.13)then
       xx = (q2-q2max)/10.6e04
       else if (na.gt.100)then
       xx = (q2-q2max)/2.8e04
       else
       xx = (q2-q2max)/15.0e+4
       endif
      if (xx.lt.150.) rhl = exp(-xx)
      endif
c     calc series for f(q) , first determine costheta, then interpolate
      xx = 1.-2.*q2/q2max
      if(q2.ge.q2max)xx=-1.
      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
c      multiply f(qmax)*(exp(q2-q2max)*R2) and remove proton size
      do 70 i=1,nff
      if(q2.gt.q2max)ff(i) = ff(i)*rhl
      if(nifty(15).eq.0)ff(i) = ff(i)*(1.+q2/18.2/h2)**2
 70   continue
c
c *** Special case for C13 (spin form factor from single neutron)
c
      if(nA.eq.13.and.nff.eq.4)then
      ff(3)=0.0
c for single neutron in 1p1 shell divide by (2j+1) and N (multiplied later by N)
      ff(4)=ff(2)/(nA-nZ)/3.
      endif

      return


      end
