c
c    include some modifications for the (e,e') case...
c


csinpmc234567
      program pwba
      implicit real*8 (a-h,o-z)
c
      real*8 step,r,hst,fac,s2(-1:1),fff(3)
      character*20 rpiv
      character*60 ipivot
      character*26 infile,outfile,fil1,fil2,fil3,fil4,fil5,fil6
      PI=dacos(-1.d0)
      ALPHA=1.d0/137.036d0
      hc=197.33d0
c      print*,'mete escale'
c      read(*,*)esc
c      hc=hc/esc
      amassp=938.3d0/hc
c      amassp=0.1d0
      alph=137.d0
      ua=931.4944/hc
      sqr2=dsqrt(2.d0)
      sqr6=dsqrt(6.d0)
      fourpi=4.d0*pi
c      call cputim(it1)
      cpmfac=1.d0

      open(1,file='rfg.in',status='old')
1234  format(60A)
123   format(26A)
c      call cputim(itbeg)
c
c
c     Entrada de datos
c
c
 9991 format(7e13.5)
      read(1,1234)ipivot      
      read(1,1234)ipivot
      read(1,123)outfile
      open(13,file='rfg.out')
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,*)z1num,z1mass,z2num,z2mass
      write(13,9991)z1num,z1mass,z2num,z2mass
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,*)nb,lb2,jb2,spft,iff,ireac,icharg,asig





      if(ireac.eq.0)then


      if(icharg.eq.1)then
      rpiv='  (e,e''p)  '
      else
      rpiv='  (e,e''n)  '
      endif


      else

      if(icharg.eq.1) then

      if(asig.gt.0.)then
      rpiv='  (nu,nu''p)  '
      else
      rpiv='  (anu,anu''p)  '
      endif

      else


      if(asig.gt.0.)then
      rpiv='  (nu,nu''n)  '
      else
      rpiv='  (anu,anu''n)  '
      endif

      endif

      endif





c
c
c

      write(13,*)nb,lb2,jb2,spft,iff,ireac,icharg,asig
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,*)Tpfin,tpffin,tpfstep,qinp,pmfin,dpm
      write(13,9991)tpfin,tpffin,tpfstep,qinp,pmfin,dpm
      wmini=wmini/hc
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,*)eki,Emiss,tick
      write(13,9991)eki,Emiss,tick
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,123)infile
      write(13,123)infile
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,*)step,rfin,nstcos,cosmin,cosmax,imp
      write(13,1111)step,rfin,nstcos,cosmin,cosmax,imp
 1111 format(2e13.6,i5.3,2e13.5,i5.3)

      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,123)fil1
      read(1,123)fil2
c      read(1,123)fil3
      write(13,*)fil1
      write(13,*)fil2
c      write(13,*)fil3
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,123)fil4
      read(1,123)fil5
      read(1,123)fil6
      write(13,*)fil4
      write(13,*)fil5
      write(13,*)fil6
      read(1,1234)ipivot      
      write(13,*)ipivot
c      read(1,1234)ipivot
c      write(13,*)ipivot
      read(1,123)fil3
      write(13,*)fil3
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,*)ipl,ipt,icc,iwoper,f2s,g1s  
      write(13,*)ipl,ipt,icc,iwoper ,f2s,g1s
c
c
c     Fin de la entrada de datos
c

      open(26,file='sigcc1.out')
c      open(11,file=infile,status='old')
      open(37,file='proton.out')
      open(38,file='neutron.out')
      open(36,file='0int.out')
      open(39,file='inclus.out')
      open(57,file='formf.out')
      open(40,file='sigma.out')


      write(*,*)' ---------------------------------------'
      write(*,*)' -  Program RFG- INCLUSIV-  J.M. UDIAS -'
      write(*,*)' -AMSTERDAM -1994----   Current v. 2.01-'
      write(*,*)' ---------------------------------------'

      write(39,*)' ---------------------------------------'
      write(39,*)' -  Program RFG- INCLUSIV-  J.M. UDIAS -'
      write(39,*)' -AMSTERDAM -1994----   Current v. 2.01-'
      write(39,*)' ---------------------------------------'


      write(39,*)'-----------------------------------------'
      write(39,*)' Fermi Momentum=',pmfin,'  MeV  '
      write(39,*)' Binding Energy=',b    ,'  MeV  '
      write(39,*)' energy incident=',eki,'  MeV'
      write(39,*)' values of response functions given for  '
      write(39,*)'    w=',tpfin,'   MeV'
      write(39,*)' to w=',tpffin, '   MeV '
      write(39,*)'(or the values constrained by kinematics)'
      write(39,*)' in steps of ',tpfstep,' MeV  '
      write(39,*)' Ef integration done in ',nstcos, '  steps '
      write(39,*)'------------------------------------------'


      write(*,*)'-----------------------------------------'
      write(*,*)' Fermi Momentum=',pmfin,'  MeV  '
      write(*,*)' Binding Energy=',b    ,'  MeV  '
      write(*,*)' energy incident=',  eki,   '  MeV  '
      write(*,*)' values of response functions given for  '
      write(*,*)'    w=',tpfin,'   MeV'
      write(*,*)' to w=',tpffin, '   MeV '
      write(*,*)'(or the values constrained by kinematics)'
      write(*,*)' in steps of ',tpfstep,' MeV  '
      write(*,*)' Ef integration done in ',nstcos, '  steps '
      write(*,*)'------------------------------------------'


      write(57,*) '--------------------------------------'
      write(57,*)' Factores de Forma para el caso',rpiv


      
      idip=1

      if(iff.lt.0)then

      write(57,*)'dipole dependence set to 1'
      idip=0
      iff=iabs(iff)

      endif





      if(iff.eq.2)write(57,*)' F1 set to 0 in the calculation'
      if(iff.eq.3)write(57,*)' F2 set to 0 in the calculation'
      if(iff.eq.4)write(57,*)' F1+kF2 set to 0 in the calc.'
      if(iff.eq.5)write(57,*)' F1=F2=0 only G1=1 enters'



      write(57,*)'---------------------------------------'
      write(57,*)' qnu2 (fm-2)     ff1     ff2        G1      fcoup   ',
     , '  gm/ge     test '    
      write(57,*)'-----------------------------------------------------'





      write(38,*) '--------------------------------------'
      write(38,*)' Funciones de Respuesta para el caso ',rpiv
      write(38,*)'-----------------------------------------------------'
      write(38,*)' Pm (MeV)  ispf2   mb2      RL        RT         RLT',
     ,   '      RTT       RATT   '
      write(38,*)'-----------------------------------------------------'
      write(38,*) '--------------------------------------'

4799  format(' Pm    spf2 mb2 comp         gam                 sig     ',
     ,'              five  ')
4759  format('-----------------------------------------------------')





      write(26,*)' cross section sig_cc1 of de Forest in fm2 ',
     ,                                             rpiv
      write(26,*)'  pm    gamma  phi  theta  cc1',
     ,           '    L       T           LT'
      write(26,*)'-------------------------------------------------'



      if(iff.eq.4)print*,'F_1 + k F_2 = 0'
      if(iff.eq.5)print*,'-k F_2 /2/M = 0'

c
c     masa nucleo residual y blanco
c


c
c     energia total, momento y energia cinetica del proton final
c
 

      emiss=emiss/hc
      B=emiss
      eki=eki/hc
      pmfmin=pmfin/hc


c
c     fermi momentum in fm-1
c
      pfermi=pmfmin
c      print*,'pfermi=',pfermi*hc

c
c     normalization factor for fermi gas
c
      ffg0=3.d0/4.d0/pi/pfermi**3

      eki2=eki*eki
      do 12345 wloop=tpfin,tpffin,tpfstep

      w=wloop/hc
      ekf=eki-w
      ekf2=ekf*ekf


      if(qinp.gt.0)then
        q=qinp/hc
        q2=q*q
        cosnu=(eki2+ekf2-q2)/2/eki/ekf
        ang=dacos(cosnu)/pi*180
      else
        ang=dabs(qinp)
        cosnu=dcos(ang*pi/180.d0)
        q2=eki2+ekf2-2*eki*ekf*cosnu
        q=dsqrt(dabs(q2))
      endif

c       print*,w

       aa=2*(b-w)
       bb=(b-w)**2-q*q
       root=4*aa*aa*bb*bb+4*(4*q*q-aa*aa)*(bb*bb+4*q*q*amassp**2)
       root=dsqrt(root) 

       den=8*q*q-2*aa*aa
       e1=2*aa*bb+root
       e2=2*aa*bb-root
       e1=e1/den
       e2=e2/den
c       root=4*aa*aa*bb*bb+4*(4*q*q-aa*aa)*(bb*bb+4*q*q*amassp**2)

c       e3=
c       print*,'e1,e1 ',e1*hc,e2*hc

       enfermi=dsqrt(amassp**2+pfermi**2)

       efmin=max(e1,w+amassp-b,amassp,enfermi)
       efmax=w+dsqrt(amassp**2+pfermi**2)-b 

      

      pmin=efmin
      pmax=efmax

      sekf=(pmax-pmin)/nstcos
      nsimps=(pmax-pmin)/sekf
      nsimps=2*(nsimps/2+1)
      sekf=(pmax-pmin)/nsimps

      pmax=pmax+sekf/30.d0

c      print*,'amin, amax, sa',pmin*hc,pmax*hc,sekf*hc
      if(pmin.gt.pmax) goto 12345
      print*,'-------------------------------------------'
      print*,' at w=',w*hc,' MeV'
      write(*,*)' the 3-momentum transferred is q=',q*hc, '  MeV'
      print*,' the electron scattering angle is tet=',ang, '  degrees'
      print*,' the min. and max. value of Tf are:'
      print*,(pmin-amassp)*hc,(pmax-amassp)*hc
      suml0=0.
      sumt0=0.
      sumttp0=0.
      suml0t=0.
      sumt0t=0.
      sumttpt0=0.

      do 1400 icharg=1,2

      ffg=ffg0*z1num
      if(iabs(icharg).eq.2)ffg=ffg0*(z1mass-z1num)

      fac=2.d0
      suml=0.d0
      sumt=0.d0
      sumttp=0.d0

      sumlt=0.d0
      sumtt=0.d0
      sumttpt=0.d0

      kstart=1
      do 1300 bet=pmin,pmax,sekf

      self=0.8*450-amassp/bet*550
      self=self/hc

      epf=bet
      write(42,*)bet*hc,epf*hc,w*hc,self*hc
      epf2=epf**2
c      print*,epf2-amassp**2
      akpf=dsqrt(epf2-amassp**2)
      tpf=epf-amassp
      akpf2=akpf*akpf




       isimps=isimps+1

c       self=0.
c
      ekf=eki-w
      Epi=-w+epf+self
c      printi,'epi+b-amassp ',epi+b-amassp,epi+b,epf,w   
      akr=dsqrt(dabs((epi+b)**2-amassp*amassp))

      akr=akr+1.d-11
      pm=akr
      pm2=pm*pm


c      print*,'ekf ',ekf*hc

c
c     calculo de w y Epi
c

      w2=w*w
      Epi2=Epi*Epi

c
c     calculo de wbar y ebar (energias iniciales si proton libre)
c
c
c     empiezan a calcularse las variables en el sistema en que pm=0
c     ojo, todo calculado SIN RECOIL y en el sistema LABORATORIO     
c
c
      Epibar=dsqrt(pm2+amassp*amassp)
      wbar=(Epf-Epibar)
      w2bar=wbar*wbar




c
c     energia electron final
c

      ekf2=ekf*ekf




      facos=2.0d0
      sigacos=0.0d0
      ichcos=0
      sigacost=0.0d0


2393  format(/' e lept i=',f9.3,' e lept f=',f9.3,' akpf   =',f9.3,
     , '   MeV'/
     ,       '  Epf     =',f9.3,'   w     =',f9.3,'  pm    =',f9.3,
     , '   MeV'/
     ,       '  Epibar  =',f9.3,'   Tb    =',f9.3,'  Epi   =',f9.3,
     , '   MeV'/
     ,       '  Wbar    =',f9.3,'  eki   =',f9.3,' ekf  =',f9.3,
     , '   MeV'/
     ,      ' Emiss     =',f9.3,'  q     =',f9.1,'  cosgam =',f9.6,
     ,      '  cosm     =',f9.6 /)

c
c



      cosif=cosnu
      cosf=q2+akpf*akpf-pm*pm
      cosf=cosf/2.d0/q/akpf      
      cosm=akpf/akr*cosf-q/akr


      if(imp.eq.1)then
      write(*,2393)eki*hc,ekf*hc,akpf*hc,epf*hc,
     ,w*hc,pm*hc,epibar*hc,tpf*hc,epi*hc,wbar*hc,eki*hc,
     , ekf*hc,emiss*hc,q*hc,cosf,cosm
      endif

      cosgam=cosf


      cos2gam=cosgam*cosgam
      singam=dsqrt(dabs(1.d0-cos2gam))
      sin2gam=singam*singam

c
c     los factores longitudinal/espacial etc
c



      Epi2=Epi*Epi
      Epf2=Epf*Epf
      cosif=cosnu

      qnu2=w2-q2
      qnu2bar=w2bar-q2
      sin2if2=-qnu2/4.d0/eki/ekf
      cos2if2=1.d0-sin2if2
      tg2if2=sin2if2/cos2if2
      tgif2=dsqrt(dabs(tg2if2))

818   format(11f7.1)

c
c     rewrite with my notation...
c     some kinematical factors (only the terms that
c     survived to the azimuthal integration)
c
c
      Wl=qnu2*qnu2/q2/q2
      Wpar=qnu2/q2/2.d0 * (-1.d0) +tg2if2
      Wper=qnu2/q2/2.d0 * (-1.d0) +tg2if2



      Wtt=0.0d0
      Wlt=0.0d0
      Wttp=-asig* 2.d0*tgif2 * dsqrt(dabs(tg2if2-qnu2/q2))
      Wltp=0.0d0
c
c     some factors
c
      free2=(Epf+Epibar)**2
c
c     se traen los factores de forma y la constante
c     de acoplo+propagador. Si ireac=0 es el
c     caso EM si es 1 es para neutrinos
c     icharg es el indice del nucleon (1-proton, 2-neutron)
c
       call ffactor(idip,icharg,ireac,qnu2,fff(1),eff,ggg,fcoup,f2s,g1s)

C
C     MOTT'S CROSS SECTION in fm^-2 x units of fcoup (see in ffactor)
C
      sigmott=4.d0*fcoup*ekf2*cos2if2
      sm=sigmott

      if(iff.eq.2) fff(1)=0.d0
      if(iff.eq.3) eff=0.d0
      if(iff.eq.4) fff(1)=-eff
      if(iff.eq.5) then
        fff(1)=1.d-20
        eff=1.d-20
        ggg=1.d0
      endif

      if(iff.eq.6)then
      ggg=0.d0
      endif

       
      fff(2)=eff/2.d0/amassp

      f2=fff(2)**2
      f1=fff(1)**2
      f3=(eff+fff(1))**2
      fff(3)=eff+fff(1)

      g1=ggg

      g12=g1*g1
      gf=g1*(eff+fff(1))

      pfz=akpf*cosgam
      pmz=pfz-q


C      write(57,'(7e11.3)')pm*hc,f1,f2,f3,g1,g12,gf
c
c     f2 es (-kF2(q)/2Mp)**2, f1 es F1(q)**2, f3 es (f1+kF2)**2
c     The free response functions (cc1 prescription)
c

      rho2=free2*(f1-f2*qnu2bar)-q2*(f3 ) 
      rho2=rho2+g12*(-4.d0*amassp**2 -q2+free2)


      djper2=-qnu2bar*f3     
      djper2=djper2+g12*(4.*amassp**2-qnu2bar)

      djpar2=akpf2*sin2gam*4.d0*(f1-f2*qnu2bar)+djper2
      djpar2=djpar2+akpf2*sin2gam*4.d0*g12


c      djper2=0.d0

      djlt=0.d0
      djtt=0.d0
      djltp=0.0d0

      djttp=-4.d0* gf * ((epf-epibar)* akpf * cosgam - epf* q)

      adjz2=(f1-f2*qnu2bar)*(pmz+pfz)**2-
     -     f3*wbar**2+g12*(4*amassp**2-wbar**2+(pmz+pfz)**2)

      adjzr=g12*4*(epibar*pfz+epf*pmz)+
     +       f3*2.d0*(epibar-epf)*q+
     +       (f1-f2*qnu2bar)*2.d0*(epibar+epf)*(pmz+pfz)

c
c     cambio de nombres...
c      
      rcc1=rho2
      dltcc1=djlt
      percc1=djper2
      parcc1=djpar2
      djzcc1=djzz
c
c     waw and the rest of structure functions have 1/fm^2 dimension
c
c   
c      if icc<0 impose current conservation by removing J_z
c
      if(icc.lt.0)then
      adjz2=w2/q2*rho2
      adjzr=2.d0*w/q*rho2
      endif


      rho2=rho2+w2/q2*adjz2-w/q*adjzr


   
      waw=( rho2+  wper*djper2 + wpar*djpar2 + wlt*djlt )
      waw=waw+(wtt*djtt+wltp*djltp+wttp*djttp)


      sigcc1=waw*sigmott

      scc1ro=wl*rho2*sigmott
      scc1per=sigmott*wper*djper2
      scc1par=sigmott*wpar*djpar2
      scc1lt=sigmott*wlt*djlt
      scc1ltp=sigmott*wltp*djltp
      scc1ttp=sigmott*wttp*djttp
      scc1tt=0.d0

c
c     so sigcc1 has the units of fcoup x fm^(-4)
c

c      acumulate response functions
   
       suml=suml+rho2*fac
       sumt=sumt+(djper2*wper+djpar2*wpar)*fac
       sumttp=sumttp+wttp*djttp*fac

c       print*,'wper',wper,tg2if2

       sumlt=sumlt+rho2
       sumtt=sumtt+djper2*wper+djpar2*wpar
       sumttpt=sumttpt+wttp*djttp

       fac=6.d0-fac

       write(36,'(2f9.3,3e13.5,f11.5)')w*hc,pm*hc,
     ,                 rho2,djper2+djpar2,djttp,cosgam

c      
c      keep values at starting point
c
       if (kstart.eq.1)then
       a0l=rho2
       a0t=djper2*wper+djpar2*wpar
       a0tp=wttp*djttp
       kstart=0
       endif

1300   continue


c   
c     seccion eficaz en 10 -42 cm2/MeV (only for neutral current)
c

       factor=ffg*(2.d0*pi)/q/4.d0

      

c
c     ffg is the corresponding normalization for RFG
c     factor is in fm^4 so esig is in the units of fcoup
c
c      factor=factor*fcoup
c
c     convert to MeV-1
c
      factor=factor/hc
c
c     Simpsons 3/4 integrals
c


      suml=suml0+factor*(suml-rho2-a0l)*sekf/3.d0/wl
      sumt=sumt0+factor*(sumt-wper*djper2-wpar*djpar2-a0t)*
     *        sekf/3.d0/wper
      sumttp=sumttp0+factor*(sumttp-wttp*djttp-a0tp)*sekf/3.d0/wttp

c
c     same as before but now with trapezoidal's rule
c
      sumlt=suml0t+factor*(sumlt-0.5*rho2-0.5*a0l)*sekf/wl
      sumtt=sumt0t+factor*(sumtt-0.5*(wper*djper2+wpar*djpar2)-0.5*a0t)*
     *              sekf/wper
      sumttpt=sumttpt0+
     + factor*(sumttpt-0.5*wttp*djttp-0.5*a0tp)*sekf/wttp
     

c
c     include a factor to compare with Horowitz's
c

      if(icharg.eq.1)then
        write(37,'(f8.3,6e12.4)')w*hc,suml,sumlt,sumt,sumtt,
     ,                                sumttp,sumttpt
      endif
      if(icharg.eq.2)then
        write(39,'(f8.3,6e12.4)')w*hc,suml,sumlt,sumt,sumtt,
     ,                                sumttp,sumttpt
        write(38,'(f8.3,6e12.4)')w*hc,suml-suml0,sumlt-suml0t,
     ,   sumt-sumt0,sumtt-sumt0t,
     ,                                sumttp-sumttp0,
     ,  sumttpt-sumttpt0

      print*,'and the total repsonse function is',suml+wper*
     *sumt+wttp*sumttp,'   1./MeV'
c
c    now printing the cross-section in fm^2/MeV
c
      write(40,'(f12.3,4e13.5)')w*hc,
     ,               sm*(wl*suml+wper*sumt+wttp*sumttp),sm,wl,suml

      endif
       suml0=suml
       suml0t=sumlt
       sumt0=sumt
       sumt0t=sumtt
       sumttpt0=sumttpt
       sumttp0=sumttp

1400   continue
      



    1 continue


12345   continue

      close(39)
      close(37)
      close(38)
      close(46)
      close(47)
      close(13)
      close(57)
      close(67)
      close(68)
      end



      subroutine ffactor(idip,icharg,ireac,qnu2,ef1,ef2,eg1,fcoup,
     ,                                                     f2s,g1s)
      implicit real*8 (a-h,o-z)
      real*8  eff1(2),eff2(2),
     ,         ff1(2),ff2(2),f1(0:8),f2(0:8),g1(0:8),gg1(2)

 

      akp=1.79285d0
      akn=-1.91304d0


      alph=1.d0/137.04
      pi=3.1415926    

     

      sin2w=0.2330

c
c     Gf  in MeV^(-2)
c

      Gf=1.1664e-11
      sqr3=sqrt(3.0)
      amass=939.d0
      hc=197.3d0

      am2v=840.d0**2
      am2a=1030.d0**2

c      f2s=-0.21
c      g1s=-0.19
c      f2s=0.
c      g1s=0.
      f1s=0.
      f1(0)=1.+f1s
      f2(0)=f2s+(akp+akn)
      g1(0)=g1s+0.55/3.



      f1(3)=0.5
      f2(3)=0.5*(akp-akn)
      g1(3)=0.63

      f1(8)=0.5*sqr3
      f2(8)=0.5*sqr3*(akp+akn)
      g1(8)=0.275/sqr3

      z0v=-0.5
      z0a=0.5

      z3v=1.-2*sin2w
      z3a=-1.

      z8v=1./sqr3*(1.-2*sin2w)
      z8a=-1./sqr3

c      print*,z0v,z3v,z8v,f2(0),f2(3),f2(8)
c      write(10,*)'proton-nu   proton anu  neutron nu  neutron-anu'

      Q2=-qnu2*hc*hc

      dip=1.+q2/am2v
      dip=1./dip/dip

      dips=1.+Q2/am2a
      dips=1./dips/dips

      tau=Q2/4./amass/amass
c------------------------------------
c
c     factores de forma EM del neutron y proton
c
      if (icharg.eq.1)then
        anom=akp
        ge0=1.0
      else
        anom=akn
        ge0=0.0
      endif

      Gm0=ge0+anom

      if(idip.eq.0)then
        dip=1.d0
        dips=1.d0
      endif

      gm=gm0*dip
      ge=ge0*dip

      if(icharg.eq.2)then
      ge=-anom*tau/(1.+5.6*tau)*dip
      endif

      eff2(icharg)=(gm-ge)/(1.+tau)/anom
      eff1(icharg)=(ge+tau*gm)/(1.+tau)

      gmt=eff1(icharg)+anom*eff2(icharg)
      get=eff1(icharg)-tau*anom*eff2(icharg)

      f20=f2(0)/(1.+tau)*dip
      f23=f2(3)/(1.+tau)*dip
      f28=f2(8)/(1.+tau)*dip

      f10=f1(0)*dip+tau*f20
      f13=f1(3)*dip+tau*f23 
      f18=f1(8)*dip+tau*f28
     
      g10=g1(0)*dips
      g13=g1(3)*dips
      g18=g1(8)*dips

c
c    fac is for weak case (Gf/sqrt(2))^2 so for the moment is
c    in MeV^(-4)
c



      fac=(gf/sqrt(2.d0)/4.d0/pi)**2

      ff1(1)=z0v*f10+z3v*f13+z8v*f18
      ff1(2)=z0v*f10-z3v*f13+z8v*f18

      ff2(1)=z0v*f20+z3v*f23+z8v*f28
      ff2(2)=z0v*f20-z3v*f23+z8v*f28

      gg1(1)=z0a*g10+z3a*g13+z8a*g18
      gg1(2)=z0a*g10-z3a*g13+z8a*g18

c      print*,'ff1(1),ff1(2),gg1(1),gg1(2)',ff1(1),ff1(2),gg1(1),gg1(2)
c      print*,'ff2 ',ff2(1),ff2(2)

      iii=icharg

c      write(13,'(i4,5e13.5)')iii,q2,ff1(iii),ff2(iii),gg1(iii)
c     i=1, proton
c     i=2, neutron

      if(ireac.ne.0)then
      ef1=ff1(icharg)
      ef2=ff2(icharg)
      eg1=gg1(icharg)
c
c     neutral case, fac is in MeV-4 convert to 
c     cm^2/MeV fm with a factor
c     from MeV-4 to fm ^3/MeV   a factor (hc)^3
c     from fm^3/MeV to fm cm^2/MeV a factor 1.d-26
c
      fcoup=4.d0*fac*(hc)**3*1.d-26
c
c     scale to units 1.d-42 cm^2 fm / MeV
c
      fcoup=fcoup*1.d42

      else

      eg1=0.d0
      ef1=eff1(icharg)
      ef2=eff2(icharg)*anom
c
c    em case:  fcoup is in fm^4
c
      fcoup=alph*alph/qnu2/qnu2
      endif      

c      write(57,'(4f10.5,f10.6,2f10.5)')qnu2,ef1,ef2,eg1,fcoup,
c     ,                                 dip,gmt/get

      end