c234567
      program pwba
      implicit real*8 (a-h,o-z)
c      include 'param.inc'

      COMMON/FFACTORS/F10,F20,GA0,F13,F23,GA3,F18,F28,GA8

      common/num/step,rfin,kmax
      common/uno/hp(0:100)
      common/pfin/akpf,epf
      common/cons/amassp
c
      real*8 fff(3), step,v(-1:1),s0
      character*20 rpiv
      character*60 ipivot
      character*26 infile,outfile,fil1,fil2,fil3,fil4,fil5,fil6
      complex*16 zi,zero,zj1x,zj1y,zj1z,zrho,zeroz,
     , pfmas,pfmin,zpm(-1:1),zpf(-1:1),
     , zalp(-1:4),zbet(-1:4),zgam(-1:4),
     , zcros(-1:1),zsuu,zsdd,zsud,zsdu,zvdd(-1:1),zvuu(-1:1),
     , zvud(-1:1),zvdu(-1:1),qsp(-1:1),zfive(-1:4),zj1yold,
     , zcri(-1:1),zcrf(-1:1),zsfive
c----------------------------------------------------------
c
c      D E C L A R A T I O N S    F O R      T H E
c
c            P I O N   E X C H A N G E
c
c----------------------------------------------------------
      complex*16 z2rho,z2j(-1:1), zk(-1:1),zpa(-1:1),
     , zcrk(-1:1),zcrup(-1:1),zcrdp(-1:1),zvudp(-1:1),zvdup(-1:1),
     , zcruk(-1:1),zcrdk(-1:1),
     , zdotud,zdotdu,zdotup,zdotdp,zdot2,zdot3,zdot4,zdot5,
     , v0t,v3t,v8t,a0t,a3t,a8t ,  onet, twot, vft, onept, twopt,
     , ga0t, ad3t, a3t1t, a3t2t, donet,dtwot,dvft,donept,dtwopt,
     ,dga0t,dad3t,da3t1t,da3t2t,
     , v0v(-1:1),v3v(-1:1),v8v(-1:1),a0v(-1:1),a3v(-1:1),a8v(-1:1),
     , onev(-1:1), twov(-1:1), vfv(-1:1), onepv(-1:1), twopv(-1:1),
     ,donev(-1:1),dtwov(-1:1),dvfv(-1:1),donepv(-1:1),dtwopv(-1:1),
     , ga0v(-1:1), ad3v(-1:1), a3t1v(-1:1), a3t2v(-1:1),
     ,dga0v(-1:1),dad3v(-1:1),da3t1v(-1:1),da3t2v(-1:1),
     , zprop,zprol,ZJ2X,ZJ2Y,ZJ2Z,       af23,
     , a5k,a0k,a0k2,ak



      zi=dcmplx(0.d0,1.d0)
      zero=dcmplx(0.d0,0.d0)
      zeroz=zero
      PI=dacos(-1.d0)
      ALPHA=1.d0/137.036d0
      hc=197.33d0
      amassp=938.3d0/hc
c      amassp=0.1d0
      alph=137.d0
      ua=931.4944/hc
      sqr2=dsqrt(2.d0)
      sqr6=dsqrt(6.d0)
      sqr3=dsqrt(3.d0)
      fourpi=4.d0*pi
c      call cputim(it1)
      cpmfac=1.d0
      call consea

      sin2w=0.233d0

c----------------------------------------------------
c
c     N E W   C O N S T A N T S   F O R   P I O N
C
C            E X C H A N G E
c
C--------------------------------------------------

      GPI2=179.7
      AMPI2=(137.d0/hc)**2





      open(1,file='muon.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='muonx.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


       anom=1.793d0

       if(iabs(icharg).eq.2)then
        anom=-1.913d0
        amassp=939.6d0/hc
       endif

       if(icharg.lt.0)then
       icharg=iabs(icharg)
       amassp=938.3d0/hc
       endif

      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      if(ireac.ne.0)then
c        z0v=-0.5
c
c        z0a=0.5
c        z3v=1.-2*sin2w
c        z3a=-1.
c        z8v=1./sqr3*(1.-2*sin2w)
c        z8a=-1./sqr3
c      else
c        z0v=0.d0
c        z3v=1.d0
c        z8v=1.d0/sqr3
c        z0a=0.d0
c        z3a=0.d0
c        z8a=0.d0
c      endif
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,vcoul
      write(13,9991)tpfin,tpffin,tpfstep,qinp,pmfin,dpm,vcoul

c      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,nloop,ia,ib
c      write(13,1111)step,rfin,nstcos,cosmin,cosmax,imp,nloop,ia,ib
 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     time like part one body current

      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,*)z0vc,z3vc,z8vc,z0ac,z3ac,z8ac
      write(13,*)z0vc,z3vc,z8vc,z0ac,z3ac,z8ac

c     time like part two body current

      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,*)z0v,z3v,z8v,z0a,z3a,z8a
      write(13,*)z0v,z3v,z8v,z0a,z3a,z8a

c     spatial like one body current

      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,*)z0vcs,z3vcs,z8vcs,z0acs,z3acs,z8acs
      write(13,*)z0vcs,z3vcs,z8vcs,z0acs,z3acs,z8acs

c     spatial like two body current

      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,*)z0vs,z3vs,z8vs,z0as,z3as,z8as
      write(13,*)z0vs,z3vs,z8vs,z0as,z3as,z8as



      ammu=dpm/hc
      ammu2=ammu**2

      write(13,*)'final lepton mass=',ammu*hc,'  MeV'
c
c
c     Fin de la entrada de datos
c


c      open(15,file=fil3,status='new')
c      open(24,file='ancrs3.out')
c      open(25,file='ancrs4.out')
c      open(26,file='sigcc1.out')
c      open(11,file=infile,status='old')
c      open(37,file='2int.out')
c      open(38,file='0int.out')
c      open(67,file='fact1.out')
c      open(68,file='nofact1.out')
      open(36,file='0int.out')
       open(39,file='rmuonx.out')
       open(40,file='musigx.out')
      open(57,file='formf.out')
c     open(46,file='response.out')
      open(41,file='integx.out')
	open(71,file='crosstot.out')
	open(81,file='crossdif.out')



      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,*)'-----------------------------------------------------'





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

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



      emiss=emiss/hc
      B=emiss

      eki=eki/hc

c      sekf=dpm/hc
      pmfmin=pmfin/hc
      pfermi=pmfmin

      tpfin=max(tpfin/hc,0.000001d0)
      tpfstep=tpfstep/hc
      tpffin=tpffin/hc
C----------------------------------------------------
C
C      F E R M I    E N E R G Y  A N D   F A C T O R S
C
C      F O R   T H E   L O O P   I N T E G R A L S
C
C-----------------------------------------------------

      if(nloop.gt.0)then

      PPST=PFERMI/nloop

      NLOOP=PFERMI/PPST

      NLOOP=2*(NLOOP/2+1)

      PPST=PFERMI/NLOOP

      PPFERMI=PFERMI+PPST/3.D0

      else

      ppst=000000.

      endif


      write(13,*)'pfermi=',pfermi*hc,PPST*HC,NLOOP

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

      write(13,*)'tfermi=',(enfermi-amassp)*hc

      ei=eki
      ef=300.0d0/hc
      est=step/hc
      fmean=spft
      fsig=60.

      sumflux=0.d0
      i=0
      do 9333 ekk=ei,300.0d0/hc,est
c      do 9333 ekk=ei,ei,est

      flux=dexp(-((ekk*hc-fmean)/fsig)**2)
c	flux=1.0
	 sumflux=sumflux+flux
      if(i.eq.0)fluxi=flux
      i=1
      write(13,*)ekk*hc,flux
c      write(*,*)ekk*hc,flux
9333  continue


      fluxt=(sumflux-flux/2.-fluxi/2)*est
c      print*,' average flux=',fluxt
      write(13,*)'flux total=',fluxt

      resacu=0.d0
      resacu2=0.0d0


c
c     normalization factor for fermi gas
c

      ffg0=3.d0/4.d0/pi/pfermi**3

C-----------------------------------------------
C
C       I S O S P I N   P A R T   O F   T H E
C
C         M E C    C U R R E N T
C
C----------------------------------------------


       ALA2=3.D0

       AL3PLUS=0.5D0

       IF(IABS(ICHARG).EQ.2)AL3PLUS=-0.5D0

       AL3MINS=1.D0

       IF(IABS(ICHARG).EQ.2)AL3MINS=-1.D0

       AL8PLUS=SQR3/2.d0

       iekloop=0

      do 3333 ekiloop=ei,300.d0/hc,est
c      do 3333 ekiloop=ei,58.2d0/hc,est

      eki=ekiloop

      flux=dexp(-((eki*hc-fmean)/fsig)**2)/fluxt
c	flux=1.0

      eki2=eki*eki

      resumt=0.d0
      resumtt=0.d0
      facres=2.d0
      irs=1

c      tpfend=min(tpffin,eki-ammu)

c
c
	ntt=(tpffin-tpfin)/tpfstep
	ntt=2*(ntt/2+1)
c	tpfstep=(tpffin-tpfin)/ntt
c	tpffin=tpffin+tpfstep/3.d0
c
c

c	&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
c	&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
c	&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&


c	vary the kinetic energy of the outgoing nucleon
      do 12345 tpfdo=tpfin,tpffin,tpfstep

c      jnt=jnt+1

c	print*,tpfdo

      Tpf=tpfdo
        f0=1.d0
c       if(tpf.lt.enfermi-amassp) then
	if(tpf.lt.dabs(qinp)/hc) then
	print*,tpf,dabs(qinp)
        f0=0.0d0
	goto 12345
	endif
c        print*,tpfdo*hc,(enfermi-amassp)*hc

      akpf=dsqrt(tpf*tpf+2.d0*amassp*tpf)
      epf=dsqrt(akpf*akpf+amassp*amassp)      
      epf2=epf*epf
      akpf2=akpf*akpf

c	print*,pfermi,epf

c
c     Emissing=Es+Ex=w-Tf-Erec
c





c
c     energia electron incidente
c


c

c
c     calcular q y demas (la cinematica)
c

c
c     step en ekf
c




c      sigaekf=0.d0
c      ichekf=0
c      facekf=2.0d0
c      sigaekftt=0.0d0
c      sigaekfts=0.0d0
c      sigaekfst=0.0d0

c      sigaekfl=0.d0
c      sigaekft=0.d0
c      sigaekftp=0.d0
 

      wmin=-dsqrt(pfermi**2+amassp**2)+epf+B
      wmax=epf+b-amassp

      wakpf=(epf+b)**2-amassp**2-(akpf)**2
      wakpf=wakpf/2.d0/(-akpf+epf+b)
      
      weipf=(epf+b)**2-amassp**2-(eki+eki-akpf)**2
      weipf=weipf/2.d0/(akpf+epf-eki-eki+B)
c      weipf=0.0
c
c    w goes from w min to 
c

	wmini=0.0/hc
      wlow=max(wmin,wmini)
c      wlow=wmin
      weki=2*eki-akpf
c
c      write(68,*)'wmin=',wmin*hc
c      write(68,*)'wmax=',wmax*hc
c      write(68,*)'weipf=',weipf*hc,(eki+eki-akpf)*hc
c      write(68,*)'wakpf=',wakpf*hc,akpf*hc
c      write(68,*)'weki=',weki*hc

      if(akpf.lt.pfermi)      wh1=999.
c      if(akpf.lt.(eki+eki-wmin)) wh1=999.
      if(wmax.lt.weki)wh1=999.
      wh1=999.

      whigh=min(wmax,wh1)

      if(whigh.gt.akpf)      whigh=min(whigh,wakpf)

      if(whigh.gt.(eki+eki-akpf)) whigh=min(whigh,weipf)



c      write(68,*)wlow*hc,whigh*hc
      call wlimit(eki,tpf,ammu,amassp,b,enfermi,wlow,whigh)
c      write(68,*)'nuevos',wlow*hc,whigh*hc

      ekfmax=max(eki-wlow,ammu)

	pmumax=dsqrt(ekfmax*ekfmax-ammu2)

c      write(*,*)' Tpf, akpf',tpf*hc,akpf*hc

c      write(68,*)' ekf max=',tpf*hc,wlow*hc,whigh*hc



c      qmax=eki+pmumax

c      write(*,*)' q max=',qmax*hc

c      write(*,*)' pmmin=',(akpf-qmax)*hc
c      write(*,*)' pmmax=',(akpf+qmax)*hc
  
      ekfmin1=eki-whigh
	ekfmin2=ammu
	ekfmin=max(ekfmin1,ekfmin2)
	pmumin=dsqrt(ekfmin*ekfmin-ammu2)


c
c   number of points should be 2n+1
c
      



      isimps=0

c      write(*,*)ekfmax*hc,ekfmin*hc, 'max,min'

      ekfmin=1.00001*ekfmin
      ekfmax=ekfmax*0.99999d0

      pmin=ekfmin
      pmax=ekfmax
c      print*,' max min ammu ',pmax*hc,pmin*hc,ammu*hc,tpfdo*hc

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


      pmax=pmax+sekf/30.d0


	fac=2.d0
	kstart=1
	ssig=0.d0
	ssigt=0.d0
c	res0=0.d0
c	res0t=0.d0
	isimps=0



	fint=1.d0
        f0=1.d0
	if(pmin.gt.pmax) then
c	write(38,*)'error',tpf*hc
        f0=0.0
	goto 12345
	endif	

      
	ffg=ffg0*z1num

	do 1300 bet=pmin,pmax,sekf

        ekf=bet

	pmu=dsqrt(ekf*ekf-ammu2)
	pmu2=pmu*pmu

        tmuon=ekf-ammu
	
       isimps=isimps+1

	w=eki-ekf
c
      Epi=-w+epf
c	print*,'w,ekf,epi,epf',w*hc,ekf*hc,epi*hc,epf*hc

c	print*,'CUIDADO%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
c      print*,'epi+b-amassp ',epi+b-amassp
      akr=dsqrt(dabs((epi+b)**2-amassp*amassp))
c	print*,'aqui estoy',akr
      akr=akr+1.d-7
      pm=akr
       pm2=pm*pm

      if(pm-pfermi.gt.1.d-7)then
c      print*,'superrerror',(pm-pfermi)*hc
      write(36,*)'superrerror',(pm-pfermi)*hc
      endif
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

      if(iwoper.lt.0)then
      woper=w
      else
      woper=wbar
      endif



      amaxfq=akr+akpf
      aminfq=dabs(akpf-akr)

	
	aminiq=eki-pmu
	amaxiq=eki+pmu

      qmax=min(amaxfq,amaxiq)

      qmin=max(aminiq,aminfq)
c	print*,'qmax,qmin',qmax*hc,qmin*hc

c      write(36,*)'qmax etc',qmax*hc,qmin*hc
c      write(36,*)'limits etc',amaxfq*hc,aminfq*hc
      if(qmax.lt.qmin) then
c       print*, 'primigenio mistake!!!!'
c       write(38,*)'primi error',tpf*hc,qmin*hc,qmax*hc      
       goto 1300
      endif

      if(qmIN.lt.1.D-5) then
c       print*, 'Q DIVERG.!!!!  ',QMIN*HC
       write(36,*)'Q DIVERG.!!!!  ',QMIN*HC
      endif

      cosnumax=0.999999*min(1.d0,-(qmin**2-eki2-pmu2)/2./eki/pmu)
      cosnumin=0.999999*max(-1.d0,-(-eki2-pmu2+qmax**2)/2./eki/pmu)
c      print*,'cos ',cosnumin,cosnumax
c      cosnumax= 0.999999
c      cosnumin=-0.999999
c      print*,qmax*hc

      alp2=dsqrt(eki*eki+pmu2-2.d0*eki*pmu*cosnumin)
      alp1=dsqrt(eki*eki+pmu2-2.d0*eki*pmu*cosnumax)
      alp1=alp1/eki/pmu
      alp2=alp2/eki/pmu
c	print*,'limites in. y fin.',alp1,alp2

      split=alp2-alp1
C
C     INTEGRATING ON POLAR LEPTON SCATTERING ANGLE COSNU
C
       KSTCOS=1
       FACOS=2.D0
       NCOS=split/TICK
       NCOS=2*(NCOS/2+1)
       STCOS=split/NCOS

       icos=0
       sumsig=0.d0
       sumsigt=0.d0




      sinphi=0.d0
      cosphi=1.d0

       DO 12346 alp=alp1,alp2+STCOS/3.,STCOS
       icos=icos+1
c       PRINT*,COSNU,TICK,NCOS,cosnumax


      q=eki*pmu*alp
      q2=q*q
      cosnu=-(q2-eki*eki-pmu2)/2.d0/eki/pmu


      cosif=cosnu


      cosf=q2+akpf*akpf-pm*pm
c      print*,cosf,q*hc
      cosf=cosf/2.d0/q/akpf
      cosm=akpf/akr*cosf-q/akr

c      print*,cosf,cosm,q2,pm*pm-akpf*akpf






      iysh=1

c

c     un caso de test

c

      cosgam=cosf

c      print*,cosgam,cosm

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


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

      qnu2=w2-q2
      qnu2bar=w2bar-q2



      TET2=DACOS(COSNU)/2.D0
c      print*,tet2,cosnu

      SINIF2=DSIN(TET2)
      COSIF2=DCOS(TET2)

      if(cosif2.ne.0.d0)then
      TGIF2=SINIF2/COSIF2
      endif

      TG2IF2=TGIF2**2
      SIN2IF2=SINIF2**2
      COS2IF2=COSIF2**2


      cos2if=cosif**2
      sin2if=1.d0-cos2if
      sinif=dsqrt(dabs(sin2if))



      pmx=akpf*singam*cosphi
      pmy=akpf*singam*sinphi
      pmz=akpf*cosgam-q



      pfx=pmx
      pfy=pmy
      pfz=akpf*cosgam




c
c     calculamos
c     pf+,pf-,pm+,pm-
c







      qz=q
      qx=0.d0
      qy=0.d0


      pfmas=-(pfx+zi*pfy)/sqr2
      pfmin=(pfx-zi*pfy)/sqr2





      zpf(1)=pfmas
      zpf(-1)=pfmin
      zpf(0)=pfz





      zpm(1)=-(pmx+zi*pmy)/sqr2
      zpm(-1)=(pmx-zi*pmy)/sqr2
      zpm(0)=dcmplx(pmz,0.d0)





      qsp(0)=dcmplx(qz,0.d0)
      qsp(1)=-(qx+zi*qy)/sqr2
      qsp(-1)=(qx-zi*qy)/sqr2

c      write(*,*)' q   cosif   cosgam   cosm'
c      write(*,2493)cosnu,cosif,cosgam,cosm

2493  format(5e13.5)





c
c     los factores longitudinal/espacial etc
c



      q2=q*q

818   format(11f7.1)



C     FACTORES  DE LA SECCION EFICAZ



      sin2phi=sinphi*sinphi
      cos2phi=cosphi*cosphi
      sin2gam=singam*singam
      cos2gam=cosgam*cosgam
      sinphi2=2.d0*sinphi*cosphi
      cosphi2=cos2phi-sin2phi


c
c   time like form factors (one body current)
c
c      call ffactor(idip,icharg,ireac,qnu2,ef1,ef2,eg1,fcoup,
c     , f2s,g1s,    z0ac,z3ac,z8ac,z0vc,z3vc,z8vc)

	inff=9
      call jaffactor(inff,icharg,ireac,q,w,fff(1),eff,ggg)

c	print*,ef1,fff(1)
c	print*,ef2,eff
c	print*,eg1,ggg

      IF(IMP.EQ.1)THEN
       WRITE(57,*)'---------------   ',QNU2,'  ---------'
c       WRITE(57,'(3(3F20.5/))')F10,F20,GA0,F13,F23,GA3,F18,F28,GA8
       write(57,'(10e13.5)')fff(1),eff,ggg,fcoup,f2s,g1s
      ENDIF

      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.0d0
      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)
      g1t0=ggg



c-------------------------------------------
c   the current operator terms:
c   falp is for gamma_mu
c   fbet is for sigma_munu q^nu
c   fgam is for (Pf+Pi)nu
c   g1 is for gamma_nu gamma_5
c------------------------------------------

      if(iabs(icc).eq.2) then

      fgam0=0.d0
      falp0=fff(1)
      fbet0=fff(2)

      endif

      if(iabs(icc).eq.1) then

      fgam0=-fff(2)
      falp0=fff(3)
      fbet0=0.d0

      endif
c
c
      if(imp.eq.1)then
      write(57,*)'time like'
      write(57,'(10e13.5)')fgam0,falp0,fbet0
      endif
c
c   spatial -like form factors
c

      call ffactor(idip,icharg,ireac,qnu2,ef1,ef2,eg1,fcoup,
     , f2s,g1s,  z0acs,z3acs,z8acs,z0vcs,z3vcs,z8vcs)
	inff=9
      call jaffactor(inff,icharg,ireac,q,w,fff(1),eff,ggg)

      IF(IMP.EQ.1)THEN
       WRITE(57,*)'---------------   ',QNU2,'  ---------'
c       WRITE(57,'(3(3F20.5/))')F10,F20,GA0,F13,F23,GA3,F18,F28,GA8
       write(57,'(10e13.5)')fff(1),eff,ggg,fcoup,f2s,g1s
      ENDIF

      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.0d0
      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



c-------------------------------------------
c   the current operator terms:
c   falp is for gamma_mu
c   fbet is for sigma_munu q^nu
c   fgam is for (Pf+Pi)nu
c   g1 is for gamma_nu gamma_5
c------------------------------------------

      if(iabs(icc).eq.2) then
      fgam=0.d0
      falp=fff(1)
      fbet=fff(2)
      endif

      if(iabs(icc).eq.1) then
      fgam=-fff(2)
      falp=fff(3)
      fbet=0.d0
      endif
c

      if(imp.eq.1)then
      write(57,*)'spatial like'
      write(57,'(10e13.5)')fgam,falp,fbet
      endif

C
C     MOTT'S CROSS SECTION in fm^-2 x units of fcoup (see in ffactor)
C
      sigmott=fcoup*pmu/eki/4.d0
      sm=sigmott
      sfac=16*eki*pmu*cos2if2
c	print*,cos2if2
      sigmott=sm*sfac
      sm=sigmott

      Ha=-2*amassp*g1/(qnu2bar-ampi2)
      ha=ha*1.0d0
      Hat0=ha
      goto 1152

C     SECCION EFICAZ DE DEFOREST SIGcc1

      Wl=qnu2*qnu2/q2/q2
      Wper=-qnu2/q2*sin2phi+tg2if2
      Wpar=-qnu2/q2*cos2phi+tg2if2
      Wlt=qnu2/q2*dsqrt(-qnu2/q2+tg2if2)


      fener=0.25d0/Epibar/Epf
      free2=(Epf+Epibar)**2


c      write(57,'(5e13.5)')pm*hc,f1,f2,f3,g1

c
c     f2 es (-kF2(q)/2Mp)**2, f1 es F1(q)**2, f3 es (f1+kF2)**2
c

      rho2=free2*(f1-f2*qnu2bar)-q2*f3
      rho2=rho2*fener
      djper2=-qnu2bar*f3
      djpar2=akpf2*sin2gam*4.d0*(f1-f2*qnu2bar)+djper2
      djper2=djper2*fener
      djpar2=djpar2*fener
      djlt=akpf*singam*4.d0*dsqrt(free2)*(f1-f2*qnu2bar)
      djlt=djlt*fener*cosphi

      djzz=-w2bar*f3+((pmz+pfz)**2)*(f1-f2*qnu2bar)
      djzz=djzz*fener

c
c     cambio de nombres...
c

      rcc1=rho2
      dltcc1=djlt
      percc1=djper2
      parcc1=djpar2
      djzcc1=djzz

c

      waw=( wl*rho2 + wper*djper2 + wpar*djpar2 + wlt*djlt )

      sigcc1=waw*sigmott
      scc1ro=wl*rho2*sigmott
      scc1per=sigmott*wper*djper2
      scc1par=sigmott*wpar*djpar2
      scc1lt=sigmott*wlt*djlt

c      write(26,'(3f9.4,e13.5)')f1,f2,f3,sigmott
c      write(26,3664)pm*hc,gamma,phi,thet,sigcc1,scc1ro,
c     , scc1per+scc1par,scc1lt

3664  format(4f7.1,5e10.3)
1152  CONTINUE

c--------------------------------------
c         factores cinematicos
c--------------------------------------

      Sl=1.d0
      St=tg2if2-qnu2/2.d0/q2
      Slt=-dsqrt(dabs(tg2if2-qnu2/q2))*2.d0

      slt=0.d0

      Sl0=w2/q2
      Sl1=w/q*2.d0
      Att=-asig*2.d0*dsqrt(dabs(tg2if2-qnu2/q2))*tgif2

      Stt1=qnu2/q2
      Stt2=stt1/2.d0

c
c--------------------------------------
c
c
c      write(49,'(7f10.5)')St,Slt,Sl0,Sl1,Att,Stt1,Stt2



      ind=0

      Rho2=0.d0
      adjzr=0.d0
      adjz2=0.d0
      djpar2=0.d0
      djper2=0.d0
      djttp=0.d0
      f1=(epibar+amassp)
      f2=(epf+amassp)

      fnorm=dsqrt(f1*f2)
      srf=fnorm

      dot=(pmx*pfx+pmy*pfy+pmz*pfz)/f1/f2

      zcros(0)=-zi*(zpm(-1)*zpf(1)-zpm(1)*zpf(-1))/f1/f2
      zcros(1)=-zi*(zpm(0)*zpf(1) -zpm(1) *zpf(0))/f1/f2
      zcros(-1)=-zi*(zpm(-1)*zpf(0)-zpm(0) *zpf(-1))/f1/f2


      do 5999 ispi2=-1,1,2
      do 5999 ispf2=-1,1,2

       m2=ispf2-ispi2
       m=m2/2

c
c     las 4 funciones del calculo analitico
c     que dependen de cada valor de spin
c


      v(1)=0.d0
      v(0)=0.d0
      v(-1)=0.d0


      s0=0.d0
      if(ispi2.eq.ispf2)s0=1.d0


      v(m)=sqr6*(-1)**(iabs(1-ispf2)/2)*
     *                  sex3j(2,1,1,m2,ispi2,-ispf2)



      zcri(0)=  zi*(zpm(-1) * v(1)-zpm(1) * v(-1))
      zcri(1)=  zi*(zpm(0)  * v(1)-zpm(1) * v(0))
      zcri(-1)= zi*(zpm(-1) * v(0)-zpm(0) * v(-1))


      zcrf(0)=  zi*(zpf(-1) * v(1)-zpf(1) * v(-1))
      zcrf(1)=  zi*(zpf(0)  * v(1)-zpf(1) * v(0))
      zcrf(-1)= zi*(zpf(-1) * v(0)-zpf(0) * v(-1))





c-------------------------------------------
c     matrix elements (layer one)
c------------------------------------------

c

      zsuu=s0
      zsdd=s0*dot+zi*(zcros(-m)*v(m)*(-1)**(m))
      zsud=(-1)**m * v(m)* zpm(-m) /f1
      zsdu=(-1)**m * v(m)* zpf(-m) /f2




      do 7557 iv=-1,1,1

      zvdu(iv)=(zpf(iv)*s0-zi*zcrf(iv))/f2

      zvud(iv)=(zpm(iv)*s0+zi*zcri(iv))/f1

      zvuu(iv)=v(iv)

      zvdd(iv)=zpf(iv)*(-1)**m*v(m)*zpm(-m)/f1/f2+
     +         zpm(iv)*(-1)**m*v(m)*zpf(-m)/f1/f2-
     -         v(iv)*dot-
     -         zi*zcros(iv)*s0

7557  continue


C------------------------------------------
C
C        I N C L U D E    M E C     C U R R E N T
C
C-----------------------------------------

       GA=2.D0*GA3
       FWEAK=2.D0*SIN2W

c       INCLUDE 'MECX.INC'

c234567
       wm=wbar

       cf2=wm/2.d0/amassp


       onet=zeroz
       twot=zeroz
       vft=zeroz
       onept=zeroz
       twopt=zeroz
       ga0t=zeroz
       ad3t=zeroz
       a3t1t=zeroz
       a3t2t=zeroz
c
c      new term
c

       a0k=zeroz
       a5k=zeroz
       a0k2=zeroz
       ak=zeroz

       do 211 i=-1,1

       onev(i)=zeroz
       twov(i)=zeroz
       vfv(i)=zeroz
       onepv(i)=zeroz
       twopv(i)=zeroz
       ga0v(i)=zeroz
       ad3v(i)=zeroz
       a3t1v(i)=zeroz
       a3t2v(i)=zeroz

211    continue

       facl=4.d0

       ILOOP=0

       akx=qx/wm
       aky=qy/wm
       akz=qz/wm

       zk(+1)=qsp(+1)/wm
       zk(-1)=qsp(-1)/wm
       zk(0)=qsp(0)/wm

       zcruk( 0)=zi*(zvuu(-1)*zk( 1)-    zvuu( 1)*zk(-1))
       zcruk( 1)=zi*(zvuu( 0)*zk( 1)-    zvuu( 1)*zk( 0))
       zcruk(-1)=zi*(zvuu(-1)*zk( 0)-    zvuu( 0)*zk(-1))


       zcrdk( 0)=zi*(zvdd(-1)*zk( 1)-    zvdd( 1)*zk(-1))
       zcrdk( 1)=zi*(zvdd( 0)*zk( 1)-    zvdd( 1)*zk( 0))
       zcrdk(-1)=zi*(zvdd(-1)*zk( 0)-    zvdd( 0)*zk(-1))

       IF(NLOOP.EQ.0) GOTO 3211

       if(ia.eq.0)then

       z2rho=zeroz
       z2j(1)=zeroz
       z2j(-1)=zeroz
       z2j(0)=zeroz

       goto 4252

       endif

       do 2211 AKP=Ppst,Ppfermi,Ppst

       ILOOP=Iloop+1

       akp2=akp*akp
       EP=DSQRT(AKP2+AMASSP**2)

c
c     diagram a
c

       epq=epf
       akpq2=epq*epq-amassp*amassp
       akpq=dsqrt(akpq2)


       pnum=2*amassp**2-ampi2-2*ep*epq+2*akp*akpq
       pden=2*amassp**2-ampi2-2*ep*epq-2*akp*akpq


c
c     factors: propagator, akp2=p^2 and facsimpson's
c
c       print*,'pnum pden ',pnum,pden
c       print*,ampi2,ep,epq,akp,akpq

       zprop=akp*zi*pi/akpq*dlog(pnum/pden)
       zprol=zprop*facl

       pfc=2.*akp/dlog(pnum/pden)+
     +    (2*ep*epq+ampi2-2*amassp**2)/2.d0/akpq

       pax=pfx/akpf*pfc
       pay=pfy/akpf*pfc
       paz=pfz/akpf*pfc


       dotk=akx*pax+aky*pay+akz*paz

       zpa(1)=zpf(1)/akpf*pfc
       zpa(-1)=zpf(-1)/akpf*pfc
       zpa(0)=zpf(0)/akpf*pfc

       zcrk( 0)=-zi*(zk(-1)*zpa( 1)-zk( 1)*zpa(-1))
       zcrk( 1)=-zi*(zk( 0)*zpa( 1)-zk( 1)*zpa( 0))
       zcrk(-1)=-zi*(zk(-1)*zpa( 0)-zk( 0)*zpa(-1))




       zdot2=zeroz

       do 239 i=-1,1
       zdot2=zdot2+(-1)**i*zvuu(i)*zcrk(-i)
239    continue

       zdot3=zeroz

       do 240 i=-1,1
       zdot3=zdot3+(-1)**i*zvdd(i)*zcrk(-i)
240    continue

       zdot4=zeroz

       do 241 i=-1,1
       zdot4=zdot4+(-1)**i*zvud(i)*zk(-i)
241    continue

       zdot5=zeroz

       do 242 i=-1,1
       zdot5=zdot5+(-1)**i*zvdu(i)*zk(-i)
242    continue




       epa=epibar

       fn2=2.d0*ep*epa+ep*amassp+epa*amassp
       fn1=-2.d0*ep*epa+ep*amassp+epa*amassp

       fn4=-2.d0*ep*epa-ep*amassp+epa*amassp
       fn3=2.d0*ep*epa-ep*amassp+epa*amassp



       fe1=Ep+amassp
       fe2=Ep-amassp

       c0=(ep+epa)/2.d0/ep/ep/epa
       cv=(epa-ep)/2.d0/ep/ep/epa

       d0=1.d0/2.d0/ep/ep/epa/amassp
       dv=1.d0/2.d0/ep/ep/epa/amassp

c
c     f1 term , V0,V3,V8
c

       zdotud=zeroz

       do 237 i=-1,1
       zdotud=zdotud+(-1)**i*(zvud(i))*zpa(-i)
237    continue

       zdotdu=zeroz

       do 236 i=-1,1
       zdotdu=zdotdu+(-1)**i*(zvdu(i))*zpa(-i)
236    continue


       zdotup=zeroz

       do 342 i=-1,1
       zdotup=zdotup+(-1)**i*zvuu(i)*zpa(-i)
342    continue


       zdotdp=zeroz

       do 442 i=-1,1
       zdotdp=zdotdp+(-1)**i*zvdd(i)*zpa(-i)
442    continue



       zcrup( 0)=zi*(zvuu(-1)*zpa( 1)-   zvuu( 1)*zpa(-1))
       zcrup( 1)=zi*(zvuu( 0)*zpa( 1)-   zvuu( 1)*zpa( 0))
       zcrup(-1)=zi*(zvuu(-1)*zpa( 0)-   zvuu( 0)*zpa(-1))


       zcrdp( 0)=zi*(zvdd(-1)*zpa( 1)-   zvdd( 1)*zpa(-1))
       zcrdp( 1)=zi*(zvdd( 0)*zpa( 1)-   zvdd( 1)*zpa( 0))
       zcrdp(-1)=zi*(zvdd(-1)*zpa( 0)-   zvdd( 0)*zpa(-1))


       zvudp( 0)=zi*(zvud(-1)*zpa( 1)-   zvud( 1)*zpa(-1))
       zvudp( 1)=zi*(zvud( 0)*zpa( 1)-   zvud( 1)*zpa( 0))
       zvudp(-1)=zi*(zvud(-1)*zpa( 0)-   zvud( 0)*zpa(-1))


       zvdup( 0)=zi*(zvdu(-1)*zpa( 1)-   zvdu( 1)*zpa(-1))
       zvdup( 1)=zi*(zvdu( 0)*zpa( 1)-   zvdu( 1)*zpa( 0))
       zvdup(-1)=zi*(zvdu(-1)*zpa( 0)-   zvdu( 0)*zpa(-1))


c
c      new terms proportional to q_photon
c

       da0k=c0*(-zdotup+zdotdp+fe1*zsdu-fe2*zsud)
       da0k2=-cv*(zdotup-zdotdp-fe1*zsdu+fe2*zsud)
       da5k=1./2/EP*(zdotud+zdotdu-fe1*zsdd-fe2*zsuu)
       dak=1./2./EP*(-zdotup-zdotdp+fe1*zsud+fe2*zsdu)


       donet=-c0*(fe1*zsdd+fe2*zsuu-zdotud-zdotdu)

       do 238 i=-1,1
       donev(i)=-cv*(fe2*zvud(i)+fe1*zvdu(i)-zpa(i)*(zsdd+zsuu)-
     -           zi*(zcrup(i)+zcrdp(i)))
238     continue

c
c     f2 term V0,V3,V8 remember the w/2/m/m factor
c



       dtwot=     fn1*(zsuu*dotk+zi*zdot2-fe1*zdot5)
       dtwot=dtwot- fn2*(zsdd*dotk+zi*zdot3-fe2*zdot4)

       dtwot=-d0*dtwot


       do 243 i=-1,1

       dtwov(i)=fn3*( fe2*zi*zcruk(i)-zvdu(i)*dotk+
     +            zpa(i)*zdot5-zi*zsdu*zcrk(i))
       dtwov(i)=dtwov(i)+fn4*(-fe1*zi*zcrdk(i)+zvud(i)*dotk-
     -            zpa(i)*zdot4+zi*zsud*zcrk(i))
       dtwov(i)=dtwov(i)+fn1*( zsuu*zpa(i)+zi*zcrup(i)-fe1*zvdu(i))
       dtwov(i)=dtwov(i)+fn2*(-zsdd*zpa(i)-zi*zcrdp(i)+fe2*zvud(i))

       dtwov(i)=-dtwov(i)*dv

243    continue

c
c     vf3 term , (commutator term of V3)
c

       dvft=(-zsuu*fe2+zsdd*fe1-zdotud+zdotdu)
       dvft=dvft/2.d0/ep/amassp

       do 244 i=-1,1
       dvfv(i)=zpa(i)*(zsdd-zsuu)-zi*(zcrup(i)-zcrdp(i))-
     -       fe2*zvud(i)+fe1*zvdu(i)
       dvfv(i)=dvfv(i)/2.d0/ep/amassp
244    continue

c
c     f1 term (V3, commutator term)
c

       donept=-cv*(fe2*zsuu+fe1*zsdd-zdotud-zdotdu)

       do 245 i=-1,1

        donepv(i)=-zpa(i)*(zsuu+zsdd)-zi*(zcrup(i)+zcrdp(i))+
     +            fe2*zvud(i)+fe1*zvdu(i)
        donepv(i)=-c0*donepv(i)

245    continue

c
c      f2 term (remember the w/2/m/m factor)
c

       dtwopt=(zsuu-zsdd)*dotk+zi*(zdot2-zdot3)+fe2*zdot4-fe1*zdot5
       dtwopt=-dtwopt*cv

       fei1=ep+epa
       fei2=epa-ep

       do 246 i=-1,1

       dtwopv(i)=fei1*zi*(fe2*zcruk(i)-fe1*zcrdk(i))+
     +        fei1*(dotk*(zvud(i)-zvdu(i))-(zdot4-zdot5)*zpa(i)+
     +              zi*(zsud-zsdu)*zcrk(i)) +
     +        fei2*(fe2*zvud(i)-fe1*zvdu(i))     +
     +        fei2*((zsuu-zsdd)*zpa(i)+zi*(zcrup(i)-zcrdp(i)))
       dtwopv(i)=-dtwopv(i)/2.d0/ep/ep/epa

246    continue

c******************************************
c
c      the axial contribution
c
c
c     part with Ga0


      dga0t=cv*(zdotup+zdotdp-fe2*zsud-fe1*zsdu)

      do 247 i=-1,1
      dga0v(i)=fe2*zvuu(i)+fe1*zvdd(i)-zpa(i)*(zsud+zsdu)-
     -     zi*(zvdup(i)+zvudp(i))
      dga0v(i)=-c0*dga0v(i)
247   continue


c
c     Ad3 part
c

      dad3t=zdotup-zdotdp+fe2*zsud-fe1*zsdu
      dad3t=dad3t/ep/amassp

      do 248 i=-1,1
      dad3v(i)=fe2*zvuu(i)-fe1*zvdd(i)+(zsud-zsdu)*zpa(i)+
     +        zi*(zvudp(i)-zvdup(i))
      dad3v(i)=dad3v(i)/ep/amassp
248   continue

c
c     ga3 commutator part
c

      da3t1t=c0*(zdotup+zdotdp-fe2*zsud-fe1*zsdu)

      do 249 i=-1,1
      da3t1v(i)=fe2*zvuu(i)+fe1*zvdd(i)-zpa(i)*(zsud+zsdu)-
     -     zi*(zvdup(i)+zvudp(i))
       da3t1v(i)=-cv*da3t1v(i)
249    continue

c
c    f23 part
c

      af23=(zdotup+zdotdp-fe2*zsud-fe1*zsdu)/2.d0/ep/amassp**2

      da3t2t=(ep+epibar)*af23

      do 250 i=-1,1
      da3t2v(i)=(zpa(i)+zpm(i))*af23
250   continue

      onet=onet+donet*zprol
      twot=twot+dtwot*zprol
      vft=vft+dvft*zprol
      onept=onept+donept*zprol
      twopt=twopt+dtwopt*zprol
      ga0t=ga0t+dga0t*zprol
      ad3t=ad3t+dad3t*zprol
      a3t1t=a3t1t+da3t1t*zprol
      a3t2t=a3t2t+da3t2t*zprol

ccc
c     the new q_photon terms
ccc

      a0k=a0k+da0k*zprol
      a0k2=da0k2*zprol+a0k2
      a5k=a5k+da5k*zprol
      ak=ak+dak*zprol

      do 1211 i=-1,1

      onev(i)=onev(i)+donev(i)*zprol
      twov(i)=twov(i)+dtwov(i)*zprol
      vfv(i)=vfv(i)+dvfv(i)*zprol
      onepv(i)=onepv(i)+donepv(i)*zprol
      twopv(i)=twopv(i)+dtwopv(i)*zprol
      ga0v(i)=ga0v(i)+dga0v(i)*zprol
      ad3v(i)=ad3v(i)+dad3v(i)*zprol
      a3t1v(i)=a3t1v(i)+da3t1v(i)*zprol
      a3t2v(i)=a3t2v(i)+da3t2v(i)*zprol

1211   continue


      facl=6.d0-facl

2211  continue

c   3211  CONTINUE

      faci=ppst/3.d0


      ak=faci*(ak-dak*zprop)
      a5k=faci*(a5k-da5k*zprop)
      a0k=faci*(a0k-da0k*zprop)
      a0k2=faci*(a0k2-da0k2*zprop)

      onet=faci*(onet-donet*zprop)
      twot=faci*(twot-dtwot*zprop)
      vft=faci*(vft-dvft*zprop)
      onept=faci*(onept-donept*zprop)
      twopt=faci*(twopt-dtwopt*zprop)
      ga0t=faci*(ga0t-dga0t*zprop)
      ad3t=faci*(ad3t-dad3t*zprop)
      a3t1t=faci*(a3t1t-da3t1t*zprop)
      a3t2t=faci*(a3t2t-da3t2t*zprop)

      do 1212 i=-1,1

      onev(i)=faci*(onev(i)-donev(i)*zprop)
      twov(i)=faci*(twov(i)-dtwov(i)*zprop)
      vfv(i)=faci*(vfv(i)-dvfv(i)*zprop)
      onepv(i)=faci*(onepv(i)-donepv(i)*zprop)
      twopv(i)=faci*(twopv(i)-dtwopv(i)*zprop)
      ga0v(i)=faci*(ga0v(i)-dga0v(i)*zprop)
      ad3v(i)=faci*(ad3v(i)-dad3v(i)*zprop)
      a3t1v(i)=faci*(a3t1v(i)-da3t1v(i)*zprop)
      a3t2v(i)=faci*(a3t2v(i)-da3t2v(i)*zprop)

1212   continue




      twot=cf2*twot
      twopt=cf2*twopt

      do 252 i=-1,1

      twov(i)=cf2*twov(i)
      twopv(i)=cf2*twopv(i)

252   continue


      v0t=-zi*dsqrt(2.d0/3.d0)/2.d0*ala2*(f10*onet+f20*twot)

      v3t=zi*(f13*onet+f23*twot)*al3plus
      v3t=v3t+zi*(2*f13+2*ga3/ga)*vft*al3mins
      v3t=v3t+zi*(f13*onept+f23*twopt)*al3mins


c       print*,a5k,da5k,zprol
      v8t=zi*(f18*onet+f28*twot)*al8plus

      a0t=zi*dsqrt(2.d0/3.d0)/2.d0*ga0*ga0t*ala2

      a3t=zi*ga3*ga0t*al3plus-zi/ga*(f13+f23+ga*ga3)*ad3t*al3mins
      a3t=a3t+zi*ga3*a3t1t*al3mins+zi*f23/ga*a3t2t*al3mins

      v3t=v3t+zi*4*(ga3/ga+f13)*a5k*al3mins*wbar/qnu2bar
      a3t=a3t+zi*ha*(2*ak/amassp-a0k)*al3plus*wbar
      a3t=a3t-zi*ha*a0k2*al3mins*wbar

      a8t=zi*ga8*ga0t*al8plus



      do 251 i=-1,1

      v0v(i)=-zi*dsqrt(2.d0/3.d0)/2.d0*ala2*(f10*onev(i)+f20*twov(i))

      v3v(i)=zi*(f13*onev(i)+f23*twov(i))*al3plus
      v3v(i)=v3v(i)+zi*(2*f13+2*ga3/ga)*vfv(i)*al3mins
      v3v(i)=v3v(i)+zi*(f13*onepv(i)+f23*twopv(i))*al3mins

      v3v(i)=v3v(i)+zi*4*ga3/ga*a5k*al3mins*q/qnu2bar

      v8v(i)=zi*(f18*onev(i)+f28*twov(i))*al8plus

      a0v(i)=zi*dsqrt(2.d0/3.d0)/2.d0*ga0*ga0v(i)*ala2

      a3v(i)=zi*ga3*ga0v(i)*al3plus-
     -       zi/ga*(f13+f23+ga*ga3)*ad3v(i)*al3mins
      a3v(i)=a3v(i)+zi*ga3*a3t1v(i)*al3mins+
     +       zi*f23/ga*a3t2v(i)*al3mins

      v3v(i)=v3v(i)+zi*4*(ga3/ga+f13)*a5k*al3mins*q/qnu2bar
      a3v(i)=a3v(i)+zi*ha*(2*ak/amassp-a0k)*al3plus*q
      a3v(i)=a3v(i)-zi*ha*a0k2*al3mins*q

      a8v(i)=zi*ga8*ga0v(i)*al8plus

251   continue


c
c     now include the g**2 factor and add to the
c     one body current...
c

c      z2rho=v3t-a3t-fweak*(v3t+v8t/sqr3)-
c     -      0.5d0*(v0t-v8t/sqr3)+0.5*(a0t-2.d0/sqr3*a8t)


       z2rho=z0a*a0t+z3a*a3t+z8a*a8t+z0v*v0t+z3v*v3t+z8v*v8t



c       write(47,'(2(7e11.3/))')v0t,v3t,v8t,a0t,a3t,a8t

      do 4251 i=-1,1

c      z2j(i)=v3v(i)-a3v(i)-fweak*(v3v(i)+v8v(i)/sqr3)-
c     -      0.5d0*(v0v(i)-v8v(i)/sqr3)+0.5*(a0v(i)-2.d0/sqr3*a8v(i))

      z2j(i)=z3vs*v3v(i)+z3as*a3v(i)+z8vs*v8v(i)+
     +       z0vs*v0v(i)+z0as*a0v(i)+z8as*a8v(i)
c       write(47,'(i4,2(7e11.3/))')i,v0v(i),v3v(i),v8v(i),
c     ,         a0v(i),a3v(i),a8v(i)

4251  continue

4252  continue

c---------------------------------------------------------------
c       N O W   A D D     D I A G R A M   B
c---------------------------------------------------------------

       onet=zeroz
       twot=zeroz
       vft=zeroz
       onept=zeroz
       twopt=zeroz
       ga0t=zeroz
       ad3t=zeroz
       a3t1t=zeroz
       a3t2t=zeroz

c      new term
c

       a0k=zeroz
       a5k=zeroz
       a0k2=zeroz
       ak=zeroz

       do 9211 i=-1,1

       onev(i)=zeroz
       twov(i)=zeroz
       vfv(i)=zeroz
       onepv(i)=zeroz
       twopv(i)=zeroz
       ga0v(i)=zeroz
       ad3v(i)=zeroz
       a3t1v(i)=zeroz
       a3t2v(i)=zeroz

9211    continue

       facl=4.d0

       ILOOP=0


       if(ib.eq.0)then

       goto 34252

       endif




       do 32211 AKP=Ppst,Ppfermi,Ppst

       ILOOP=Iloop+1

       akp2=akp*akp
       EP=DSQRT(AKP2+AMASSP**2)

c
c     diagram B
c

       epq=epIBAR
       akpq2=epq*epq-amassp*amassp
       akpq=dsqrt(akpq2)


       pnum=2*amassp**2-ampi2-2*ep*epq+2*akp*akpq
       pden=2*amassp**2-ampi2-2*ep*epq-2*akp*akpq


c
c     factors: propagator, akp2=p^2 and facsimpson's
c
c       print*,'pnum pden ',pnum,pden
c       print*,ampi2,ep,epq,akp,akpq

       zprop=akp*zi*pi/akpq*dlog(pnum/pden)




       pfc=2.*akp/dlog(pnum/pden)+
     +    (2*ep*epq+ampi2-2*amassp**2)/2.d0/akpq

c       if(akpq.lt.0.001)then
c        pfc=0.d0
c        zprop=akpf2*zi*4.d0*pi/(2*amassp**2-2*ep*epq-ampi2)
c       endif


       zprol=zprop*facl
       pax=pMx/pM*pfc
       pay=pMy/pM*pfc
       paz=pMz/pM*pfc


       dotk=akx*pax+aky*pay+akz*paz

       zpa(1)=zpM(1)/PM*pfc
       zpa(-1)=zpM(-1)/PM*pfc
       zpa(0)=zpM(0)/PM*pfc

       zcrk( 0)=-zi*(zk(-1)*zpa( 1)-zk( 1)*zpa(-1))
       zcrk( 1)=-zi*(zk( 0)*zpa( 1)-zk( 1)*zpa( 0))
       zcrk(-1)=-zi*(zk(-1)*zpa( 0)-zk( 0)*zpa(-1))


c       write(47,'(7e11.4)')pfc,pnum,pden,ep,epq,akpq,akp

C *********************   HERE WE FOLLOW ****************


       zdot2=zeroz

       do 3239 i=-1,1
       zdot2=zdot2+(-1)**i*zvuu(i)*zcrk(-i)
3239    continue

       zdot3=zeroz

       do 3240 i=-1,1
       zdot3=zdot3+(-1)**i*zvdd(i)*zcrk(-i)
3240    continue

       zdot4=zeroz

       do 3241 i=-1,1
       zdot4=zdot4+(-1)**i*zvud(i)*zk(-i)
3241    continue

       zdot5=zeroz

       do 3242 i=-1,1
       zdot5=zdot5+(-1)**i*zvdu(i)*zk(-i)
3242    continue




       epa=epf

       fn2=2.d0*ep*epa+ep*amassp+epa*amassp
       fn1=-2.d0*ep*epa+ep*amassp+epa*amassp

       fn4=-2.d0*ep*epa-ep*amassp+epa*amassp
       fn3=2.d0*ep*epa-ep*amassp+epa*amassp



       fe1=Ep+amassp
       fe2=Ep-amassp

       c0=(ep+epa)/2.d0/ep/ep/epa
       cv=(epa-ep)/2.d0/ep/ep/epa

       d0=1.d0/2.d0/ep/ep/epa/amassp
       dv=1.d0/2.d0/ep/ep/epa/amassp

c
c     f1 term , V0,V3,V8
c

       zdotud=zeroz

       do 3237 i=-1,1
       zdotud=zdotud+(-1)**i*(zvud(i))*zpa(-i)
3237    continue

       zdotdu=zeroz

       do 3236 i=-1,1
       zdotdu=zdotdu+(-1)**i*(zvdu(i))*zpa(-i)
3236    continue


       zdotup=zeroz

       do 3342 i=-1,1
       zdotup=zdotup+(-1)**i*zvuu(i)*zpa(-i)
3342    continue


       zdotdp=zeroz

       do 3442 i=-1,1
       zdotdp=zdotdp+(-1)**i*zvdd(i)*zpa(-i)
3442    continue



       zcrup( 0)=zi*(zvuu(-1)*zpa( 1)-   zvuu( 1)*zpa(-1))
       zcrup( 1)=zi*(zvuu( 0)*zpa( 1)-   zvuu( 1)*zpa( 0))
       zcrup(-1)=zi*(zvuu(-1)*zpa( 0)-   zvuu( 0)*zpa(-1))


       zcrdp( 0)=zi*(zvdd(-1)*zpa( 1)-   zvdd( 1)*zpa(-1))
       zcrdp( 1)=zi*(zvdd( 0)*zpa( 1)-   zvdd( 1)*zpa( 0))
       zcrdp(-1)=zi*(zvdd(-1)*zpa( 0)-   zvdd( 0)*zpa(-1))


       zvudp( 0)=zi*(zvud(-1)*zpa( 1)-   zvud( 1)*zpa(-1))
       zvudp( 1)=zi*(zvud( 0)*zpa( 1)-   zvud( 1)*zpa( 0))
       zvudp(-1)=zi*(zvud(-1)*zpa( 0)-   zvud( 0)*zpa(-1))


       zvdup( 0)=zi*(zvdu(-1)*zpa( 1)-   zvdu( 1)*zpa(-1))
       zvdup( 1)=zi*(zvdu( 0)*zpa( 1)-   zvdu( 1)*zpa( 0))
       zvdup(-1)=zi*(zvdu(-1)*zpa( 0)-   zvdu( 0)*zpa(-1))



c
c      new terms proportional to q_photon
c

       da0k=c0*(-zdotup+zdotdp+fe1*zsdu-fe2*zsud)
       da0k2=-cv*(zdotup-zdotdp-fe1*zsdu+fe2*zsud)
       da5k=1./2/EP*(zdotud+zdotdu-fe1*zsdd-fe2*zsuu)
       dak=1./2./EP*(-zdotup-zdotdp+fe1*zsud+fe2*zsdu)





       donet=-c0*(fe1*zsdd+fe2*zsuu-zdotud-zdotdu)

       do 3238 i=-1,1
       donev(i)=-cv*(fe2*zvdu(i)+fe1*zvud(i)-zpa(i)*(zsdd+zsuu)+
     +           zi*(zcrup(i)+zcrdp(i)))
3238     continue

c
c     f2 term V0,V3,V8 remember the w/2/m/m factor
c



       dtwot=       fn1*(-zsuu*dotk+zi*zdot2+fe1*zdot4)
       dtwot=dtwot+ fn2*(zsdd*dotk-zi*zdot3-fe2*zdot5)

       dtwot=-d0*dtwot


       do 3243 i=-1,1

       dtwov(i)=fn3*( fe2*zi*zcruk(i)+zvud(i)*dotk-
     -            zpa(i)*zdot4-zi*zsud*zcrk(i))
       dtwov(i)=dtwov(i)+fn4*(-fe1*zi*zcrdk(i)-zvdu(i)*dotk+
     +            zpa(i)*zdot5+zi*zsdu*zcrk(i))
       dtwov(i)=dtwov(i)+fn1*(-zsuu*zpa(i)+zi*zcrup(i)+fe1*zvud(i))
       dtwov(i)=dtwov(i)+fn2*( zsdd*zpa(i)-zi*zcrdp(i)-fe2*zvdu(i))

       dtwov(i)=-dtwov(i)*dv

3243    continue

c
c     vf3 term , (commutator term of V3)
c

       dvft=(zsuu*fe2-zsdd*fe1-zdotud+zdotdu)
       dvft=dvft/2.d0/ep/amassp

       do 3244 i=-1,1
       dvfv(i)=-zpa(i)*(zsdd-zsuu)-zi*(zcrup(i)-zcrdp(i))+
     +           fe2*zvdu(i)-fe1*zvud(i)
       dvfv(i)=dvfv(i)/2.d0/ep/amassp
3244    continue

c
c     f1 term (V3, commutator term)
c

       donept=-cv*(-fe2*zsuu-fe1*zsdd+zdotud+zdotdu)

       do 3245 i=-1,1

        donepv(i)=zpa(i)*(zsuu+zsdd)-zi*(zcrup(i)+zcrdp(i))-
     -            fe2*zvdu(i)-fe1*zvud(i)
        donepv(i)=-c0*donepv(i)

3245    continue

c
c      f2 term (remember the w/2/m/m factor)
c

       dtwopt=(zsuu-zsdd)*dotk-zi*(zdot2-zdot3)+fe2*zdot5-fe1*zdot4
       dtwopt=-dtwopt*cv

       fei1=ep+epa
       fei2=epa-ep

       do 3246 i=-1,1

       dtwopv(i)=fei1*zi*(-fe2*zcruk(i)+fe1*zcrdk(i))+
     +        fei1*(-dotk*(zvud(i)-zvdu(i))+(zdot4-zdot5)*zpa(i)+
     +              zi*(zsud-zsdu)*zcrk(i)) +
     +        fei2*(-fe2*zvdu(i)+fe1*zvud(i))     +
     +        fei2*((zsuu-zsdd)*zpa(i)-zi*(zcrup(i)-zcrdp(i)))
       dtwopv(i)=-dtwopv(i)/2.d0/ep/ep/epa

3246    continue

c******************************************
c
c      the axial contribution
c
c
c     part with Ga0


      dga0t=cv*(zdotup+zdotdp-fe1*zsud-fe2*zsdu)

      do 3247 i=-1,1
      dga0v(i)=fe2*zvuu(i)+fe1*zvdd(i)-zpa(i)*(zsud+zsdu)+
     +     zi*(zvdup(i)+zvudp(i))
      dga0v(i)=-c0*dga0v(i)
3247   continue


c
c     Ad3 part
c

      dad3t=-zdotup+zdotdp+fe1*zsud-fe2*zsdu
      dad3t=dad3t/ep/amassp

      do 3248 i=-1,1
      dad3v(i)=-fe2*zvuu(i)+fe1*zvdd(i)+(zsud-zsdu)*zpa(i)-
     -        zi*(zvudp(i)-zvdup(i))
      dad3v(i)=dad3v(i)/ep/amassp
3248   continue

c
c     ga3 commutator part
c

       da3t1t=-c0*(zdotup+zdotdp-fe1*zsud-fe2*zsdu)

       do 3249 i=-1,1
       da3t1v(i)=fe2*zvuu(i)+fe1*zvdd(i)-zpa(i)*(zsud+zsdu)-
     -     zi*(zvdup(i)+zvudp(i))
       da3t1v(i)=cv*da3t1v(i)
3249    continue

c
c    f23 part
c

      af23=(-zdotup-zdotdp+fe1*zsud+fe2*zsdu)/2.d0/ep/amassp**2

      da3t2t=(ep+epf)*af23

      do 3250 i=-1,1
      da3t2v(i)=(zpa(i)+zpf(i))*af23
3250   continue

      onet=onet+donet*zprol
      twot=twot+dtwot*zprol
      vft=vft+dvft*zprol
      onept=onept+donept*zprol
      twopt=twopt+dtwopt*zprol
      ga0t=ga0t+dga0t*zprol
      ad3t=ad3t+dad3t*zprol
      a3t1t=a3t1t+da3t1t*zprol
      a3t2t=a3t2t+da3t2t*zprol

ccc
c     the new q_photon terms
ccc

      a0k=a0k+da0k*zprol
      a0k2=da0k2*zprol+a0k2
      a5k=a5k+da5k*zprol
      ak=ak+dak*zprol

      do 31211 i=-1,1

      onev(i)=onev(i)+donev(i)*zprol
      twov(i)=twov(i)+dtwov(i)*zprol
      vfv(i)=vfv(i)+dvfv(i)*zprol
      onepv(i)=onepv(i)+donepv(i)*zprol
      twopv(i)=twopv(i)+dtwopv(i)*zprol
      ga0v(i)=ga0v(i)+dga0v(i)*zprol
      ad3v(i)=ad3v(i)+dad3v(i)*zprol
      a3t1v(i)=a3t1v(i)+da3t1v(i)*zprol
      a3t2v(i)=a3t2v(i)+da3t2v(i)*zprol

31211   continue

c       write(47,*)ispi2,ispf2
c       write(47,'(2(7e11.3/))')donet,donept,dtwot,dtwopt,dvft
c     ,                         dga0t,dad3t



      facl=6.d0-facl

32211  continue

3211  CONTINUE

      faci=ppst/3.d0


      ak=faci*(ak-dak*zprop)
      a5k=faci*(a5k-da5k*zprop)
      a0k=faci*(a0k-da0k*zprop)
      a0k2=faci*(a0k2-da0k2*zprop)

      onet=faci*(onet-donet*zprop)
      twot=faci*(twot-dtwot*zprop)
      vft=faci*(vft-dvft*zprop)
      onept=faci*(onept-donept*zprop)
      twopt=faci*(twopt-dtwopt*zprop)
      ga0t=faci*(ga0t-dga0t*zprop)
      ad3t=faci*(ad3t-dad3t*zprop)
      a3t1t=faci*(a3t1t-da3t1t*zprop)
      a3t2t=faci*(a3t2t-da3t2t*zprop)

      do 31212 i=-1,1

      onev(i)=faci*(onev(i)-donev(i)*zprop)
      twov(i)=faci*(twov(i)-dtwov(i)*zprop)
      vfv(i)=faci*(vfv(i)-dvfv(i)*zprop)
      onepv(i)=faci*(onepv(i)-donepv(i)*zprop)
      twopv(i)=faci*(twopv(i)-dtwopv(i)*zprop)
      ga0v(i)=faci*(ga0v(i)-dga0v(i)*zprop)
      ad3v(i)=faci*(ad3v(i)-dad3v(i)*zprop)
      a3t1v(i)=faci*(a3t1v(i)-da3t1v(i)*zprop)
      a3t2v(i)=faci*(a3t2v(i)-da3t2v(i)*zprop)

31212   continue




      twot=cf2*twot
      twopt=cf2*twopt

      do 3252 i=-1,1

      twov(i)=cf2*twov(i)
      twopv(i)=cf2*twopv(i)

3252   continue

c       write(47,'(2(7e11.3/))')cf2,faci,onet,onept,twot,twopt,
c     , vft,ga0t

      v0t=-zi*dsqrt(2.d0/3.d0)/2.d0*ala2*(f10*onet+f20*twot)

      v3t=zi*(f13*onet+f23*twot)*al3plus
c	write(6,*)ga
      v3t=v3t-zi*(2*f13+2*ga3/ga)*vft*al3mins
      v3t=v3t-zi*(f13*onept+f23*twopt)*al3mins

      v8t=zi*(f18*onet+f28*twot)*al8plus

      a0t=zi*dsqrt(2.d0/3.d0)/2.d0*ga0*ga0t*ala2
c	write(6,*)ga,a3t2t,al3mins

      a3t=zi*ga3*ga0t*al3plus+zi/ga*(f13+f23+ga*ga3)*ad3t*al3mins
c	write(6,*)ga,a3t2t,al3mins
      a3t=a3t-zi*ga3*a3t1t*al3mins-zi*f23/ga*a3t2t*al3mins

      a8t=zi*ga8*ga0t*al8plus

      v3t=v3t+zi*4*(ga3/ga+f13)*a5k*al3mins*wbar/qnu2bar
      a3t=a3t+zi*ha*(2*ak/amassp-a0k)*al3plus*wbar
      a3t=a3t-zi*ha*a0k2*al3mins*wbar

      do 3251 i=-1,1

      v0v(i)=-zi*dsqrt(2.d0/3.d0)/2.d0*ala2*(f10*onev(i)+f20*twov(i))

      v3v(i)=zi*(f13*onev(i)+f23*twov(i))*al3plus
      v3v(i)=v3v(i)-zi*(2*f13+2*ga3/ga)*vfv(i)*al3mins
      v3v(i)=v3v(i)-zi*(f13*onepv(i)+f23*twopv(i))*al3mins

      v8v(i)=zi*(f18*onev(i)+f28*twov(i))*al8plus

      a0v(i)=zi*dsqrt(2.d0/3.d0)/2.d0*ga0*ga0v(i)*ala2

      a3v(i)=zi*ga3*ga0v(i)*al3plus+
     +       zi/ga*(f13+f23+ga*ga3)*ad3v(i)*al3mins
      a3v(i)=a3v(i)-zi*ga3*a3t1v(i)*al3mins-
     -       zi*f23/ga*a3t2v(i)*al3mins


      v3v(i)=v3v(i)+zi*4*(ga3/ga+f13)*a5k*al3mins*q/qnu2bar
      a3v(i)=a3v(i)+zi*ha*(2*ak/amassp-a0k)*al3plus*q
      a3v(i)=a3v(i)-zi*ha*a0k2*al3mins*q

      a8v(i)=zi*ga8*ga0v(i)*al8plus

3251   continue


c
c     now include the g**2 factor and add to the
c     one body current...
c

c      z2rho=z2rho+v3t-a3t-fweak*(v3t+v8t/sqr3)-
c     -      0.5d0*(v0t-v8t/sqr3)+0.5*(a0t-2.d0/sqr3*a8t)

       z2rho=z2rho+z0a*a0t+z3a*a3t+z8a*a8t+z0v*v0t+z3v*v3t+z8v*v8t
c       write(47,'(2(7e11.3/))')v0t,v3t,v8t,a0t,a3t,a8t

      do 34251 i=-1,1

c      z2j(i)=z2j(i)+v3v(i)-a3v(i)-fweak*(v3v(i)+v8v(i)/sqr3)-
c    -      0.5d0*(v0v(i)-v8v(i)/sqr3)+0.5*(a0v(i)-2.d0/sqr3*a8v(i))


      z2j(i)=z2j(i)+z3vs*v3v(i)+z3as*a3v(i)+z8vs*v8v(i)+
     +       z0vs*v0v(i)+z0as*a0v(i)+z8as*a8v(i)

c       write(47,'(i4,2(7e11.3/))')i,v0v(i),v3v(i),v8v(i),
c     ,         a0v(i),a3v(i),a8v(i)
34251  continue

34252  continue

       ftwobody=srf*gpi2/(2*pi)**3



C---------------------------------------------------------------
C      P R E P A R E   T W O - B O D Y   C O N T R I B U T I O N
C---------------------------------------------------------------

       Z2RHO=Z2RHO*FTWOBODY

       ZJ2X=-1/SQR2*(Z2J(1)-Z2J(-1))*ftwobody

       ZJ2Y=ZI/SQR2*(Z2J(1)+Z2J(-1))*ftwobody

       ZJ2Z=Z2J(0)*ftwobody

c       WRITE(47,*)' TWO BODY A+B',pm*hc,ftwobody
c       write(47,'(8f9.3)')z2rho,zj2x,zj2y,zj2z
c       WRITE(47,*)' I. A. ',q*hc







c-----------------------------------------
c
c
c
c
c      parte con gamma_mu
c
c
c      temporal  (gamma_0)
c



       zalp(4)=srf*(zsuu+zsdd)

c
c      espacial  gamma_i
c



       do 558 im2=-2,2,2

       zalp(im2/2)=srf*(zvud(im2/2)+zvdu(im2/2))

  558  continue

c____________________________________________________
c
c      parte con sigma_munu q_nu
c
c      temporal
c

      zbet(4)=-srf*
     *        ((zvdu(0)-zvud(0))*qsp(0)-
     -         (zvdu(1)-zvud(1))*qsp(-1)-
     -         (zvdu(-1)-zvud(-1))*qsp(1))

c-----------------------------------------------
c     espacial
c-----------------------------------------------



      zbet(0)=-woper*(zvdu(0)-zvud(0))
      zbet(1)=-woper*(zvdu(1)-zvud(1))
      zbet(-1)=-woper*(zvdu(-1)-zvud(-1))



      zbet(0)=zbet(0)+  (-zvuu(-1)+zvdd(-1))*qsp(+1)-
     +                  (-zvuu(1)+zvdd(1))*qsp(-1)


      zbet(-1)=zbet(-1)+  (-zvuu(-1)+zvdd(-1))*qsp(0)-
     +                  (-zvuu(0)+zvdd(0))*qsp(-1)


      zbet(1)=zbet(1)+  (-zvuu(0)+zvdd(0))*qsp(1)-
     +                  (-zvuu(1)+zvdd(1))*qsp(0)



      zbet(0)=srf*zbet(0)
      zbet(1)=srf*zbet(1)
      zbet(-1)=srf*zbet(-1)


c-----------------------------------------------------
c  Parte con P_f + P_i
c
c   temporal
c
c----------------------------------


      zgam(4)=srf*(zsuu-zsdd)*(Epibar+Epf)


c---------------------------------------
c     espacial
c---------------------------------------



      do 6666 i=-1,1

      zgam(i)=srf*(zsuu-zsdd)*(zpm(i)+zpf(i))



6666  continue


c------------------------------------------------
c     parte con gamma_mu gamma_5
c
c     temporal
c
c-----------------------------------------------



      zfive(4)=srf*(zsud+zsdu)



      do 9985 iv=-1,1,1

      zfive(iv)=srf*(zvuu(iv)+zvdd(iv))

9985  continue

c -----------------------------------------
c  induced pseudoscalar gamma_5 term
c -----------------------------------------

      zsfive=srf*(zsud-zsdu)

      zrho=falp0*zalp(4)+fbet0*zbet(4)+fgam0*zgam(4)+
     +     g1t0*zfive(4)
c
c     H_A term
c
      zrho=zrho+HAt0*woper*zsfive

      zj1x=-1.d0/dsqrt(2.d0)*( falp*(zalp(1) - zalp(-1))+
     +                         fbet*(zbet(1) - zbet(-1))+
     +                         fgam*(zgam(1) - zgam(-1))+
     +                         g1    *(zfive(1)-zfive(-1)))

      zj1x=zj1x+ha*qx*zsfive

      zj1y=zi/dsqrt(2.d0)  *(  falp*(zalp(1) + zalp(-1))+
     +                         fbet*(zbet(1) + zbet(-1))+
     +                         fgam*(zgam(1) + zgam(-1))+
     +                         g1*    (zfive(1)+zfive(-1)))

      zj1y=zj1y+ha*qy*zsfive

      zj1z=falp*zalp(0)+fbet*zbet(0)+fgam*zgam(0)+g1*zfive(0)

      zj1z=zj1z+ha*qz*zsfive
C-------------------------------
C      ADD THE  TWO BODY  PART
C-------------------------------

      ZRHO=ZRHO+Z2RHO

      ZJ1X=ZJ1X+ZJ2X
      ZJ1Y=ZJ1Y+ZJ2Y
      ZJ1Z=ZJ1Z+ZJ2Z

c-------------------------------------------------------
c  Impresion de la corriente nuclear
c
c

c      write(47,'(i4,4e13.5)')m,s0,v(-1),v(0),v(1)
c      write(47,'(6e13.5)')pmx,pmy,pmz,pfx,pfy,pfz
c      write(47,'(f6.1,3i4.2,6f10.4)')pm*hc,
c     ,          ispf2,ispi2,4,zalp(4),zgam(4),zfive(4)

c      do 1047 imq=-1,1,1
c      write(47,'(14x,i4.2,6f10.4)')imq,zalp(imq),zgam(imq),zfive(imq)
c 1047  continue



c      write(47,9947)zrho,q/w*zj1z,q/wbar*zj1z
c      write(47,*)'          zrho              zj1x       ',
c     , '        zj1y              zj1z    '
c       write(47,'(8f9.3)')zrho,zj1x,zj1y,zj1z
c      write(47,*)'-------------------------------------'

9947  format('--- test ---',6e10.3)



c
c  now rotate zj1x and zj1y to the nucleon scattering plane, after
c  that:
c                zj1x=zjpar,
c                zj1y=zjper
c





      zj1yold=zj1y
      zj1y=-sinphi*zj1x+cosphi*zj1yold
      zj1x=cosphi*zj1x+sinphi*zj1yold


c------------------------------------------
c     how to deal with current (no) conservation
c-----------------------------------------

      if(icc.lt.0)then

      zj1z=zrho*w/q

      endif


c-------------------------------------------------------
c   construccion de las funciones de respuesta
c------------------------------------------------------

      rho2=rho2+dble(zrho*dconjg(zrho))
      adjz2=adjz2+dble(zj1z*dconjg(zj1z))
      adjzr=adjzr+ 2.d0*dble(dconjg(zrho)*zj1z)

      djpar2=djpar2+dble(zj1x*dconjg(zj1x))
      djper2=djper2+dble(zj1y*dconjg(zj1y))

      djttp=djttp+Dimag(zj1x*dconjg(zj1y))

c---------------------------------------------------------
c    impresion de estas...
c
c      write(46,'(f6.1,2i8.2,5f10.4)')pm*hc,
c     , ispf2,ispi2,rl,rt,rlt,rtt,ratt



5999  continue



6999  format(/,' Acumuladas   ',6e11.3/
     ,   '------------------------')
      if(imp.eq.1)then
      write(46,6999) pm*hc,rlac*sl,rtac*st,rltac*slt,0.0,att*rattac
      endif




c-------------------------------------------------------
c   construccion de la seccion eficaz
c

c       Waverage=Rlac*Sl+Rtac*St+Rltac*Slt+Rttac+ Att*Rattac


      wawl=((eki+ekf)**2-q2-ammu2)*rho2+
     +     ((eki2-pmu2)**2/q2-w**2+ammu2)*adjz2-
     -     ((eki+ekf)*(eki2-pmu2)/q-w*q)*adjzr
      wawl=4.*wawl

      wawt=16*eki*pmu*
     *   (eki*pmu*sin2if/2/q2-0.5*(-ekf/pmu+cosif))*
     *     (djper2+djpar2)

      wawttp=-asig*16*eki*pmu/q*
     *     (eki*ekf/pmu+pmu-(eki+ekf)*cosif)*djttp


      waw=wawl+wawt+wawttp



      factor=ffg*(2.d0*pi)**2/4.d0
c
c
c     factor for comparison with Horowitz
c
c
c      factor=epibar/epi*factor

c     ffg is the corresponding normalization for RFG
c     factor is in fm^4 so sig is in the units of fcoup
c     I.E. IN 10^(-42) CM^2 FM /MEV
C

c	print*,sm,factor,waw

      SIG=   sm*factor*waw/sfac/2.d0

      SUMSIG=SUMSIG+SIG*FACOS
      SUMSIGT=SUMSIGT+SIG
c	print*,sfac,sumsig,sumsigt,facos

      if(imp.eq.1)then
      write(38,'(10e13.6)')smott,waverage,st
      write(38,'(6e13.5/6e13.5/)')eki*hc,ekf*hc,q*hc,
     , w*hc,pmu*hc,sig,
     , cosnu,cosm,cosgam,pm*hc,akpf*hc,tpf*hc
      endif

       if (kstCOS.eq.1)then
       SIGSIG0=SIG
       kstCOS=0
       endif

       FACOS=6.D0-FACOS

c      write(36,'(6e12.4,i5)')w*hc,SIG,rho2/2.,adjz2/2.,adjzr/2.,
c     ,      djttp/2.,icos
12346 CONTINUE

      SIGS=(SUMSIG-SIG-SIGSIG0)*dabs(STCOS)/3.d0
      SIGT=(sumSIGt-0.5*SIG-0.5*SIGSIG0)*dabs(STCOS)
	
c	print*,sigs,sigt

      beta=pmu/ekf
      eta=-asig*z1num/137./beta
      fcoul=2*pi*eta/(dexp(2*pi*eta)-1.d0)
c      fcoul=1.d0
       if (kstart.eq.1)then
       sigs0=sigs
       sigs0c=sigt
       kstart=0
       endif
      write(89,'(6e13.5)')tpf*hc,tmuon*hc,sigs,sigt
     

      ssig=ssig+sigs*fac
      ssigt=ssigt+sigt
      fac=6.d0-fac
1300   continue

c	write(88,*)tpf*hc,fcoul
      ssig=(ssig-sigs0-sigs)*sekf/3.d0
      tsig=(ssigt-sigs0c/2.-sigs/2.)*sekf
      
c	print*,ssig,tsig
C
C    SSIG,TSIG IN 10^-42 CM^2/MEV
C
      WRITE(40,'(2F9.2,4E13.5,i5)')tpf*HC,eki*hc,
     ,sSIG,tSIG,tpfend,wloop,isimps

c	write(81,'(3e13.5,2i5)')tpf*hc,flux*ssig,flux*tsig,icos,isimps

      resumt=resumt+ssig*facres
      resumtt=resumtt+tsig
      if(irs.eq.1)then
       irs=0
       res0=ssig
       res0t=tsig
       print*,'irs=',irs,res0,res0t,sekf,sigs,sigt
      endif
c      write(37,*)
       facres=6.d0-facres
       
12345   continue
	f0=1.d0
      print*,ssig,res0
       res=f0*(resumt-(ssig+res0))*tpfstep*hc/3.d0
       rest=f0*(resumtt-tsig/2.-res0t/2.)*tpfstep*hc
c      res=f0*(resumt-(ssig+res0)/2.d0)*tpfstep*hc
c      rest=f0*(resumtt-(tsig+res0t)/2.d0)*tpfstep*hc

C     REST TOTAL CROSS-SECTION IN 10-42 CM^2
C
c      print*,'res=',eki*hc,res,rest
c      write(40,*)
c      write(40,*)'integrated=',res,(enfermi-amassp)*hc
      write(41,'(6e12.4)')eki*hc,res,rest,flux,
     ,res*flux/1000,rest*flux/1000
	write(71,*)eki*hc,res/1000.
	
      resacu=resacu+res*flux
      resacu2=resacu2+rest*flux


      if(iekloop.eq.0) then
      iekloop=1
      resacui=res*flux
      resacu2i=rest*flux
      endif
3333  continue

      write(41,*)
      write(41,*)'est=',est,est*hc
      write(41,*)'folded total cross-section',
     ,(resacu-resacui/2.-res*flux/2.)*est/1000.,' 10-39 cm^2'
      write(41,*)'folded total cross-section',
     ,(resacu2-resacu2i/2.-rest*flux/2.)*est/1000.,' 10-39 cm^2'




      PRINT*,' ILOOP=',ILOOP


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

      subroutine ffactor(idip,icharg,ireac,qnu2,ef1,ef2,eg1,fcoup,
     ,     f2s,g1s,z0a,z3a,z8a,z0v,z3v,z8v)
      implicit real*8 (a-h,o-z)

      COMMON/FFACTORS/F10,F20,G10,F13,F23,G13,F18,F28,G18

      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

      amasse2=(0.5d0/197.3)**2
      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=1032.d0**2

      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.645

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

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

      ef1=ff1(icharg)
      ef2=ff2(icharg)
      eg1=gg1(icharg)

      if(ireac.ne.0)then
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
c
c    em case:  fcoup is in fm^4
c
      fcoup=alph*alph/(qnu2)**2
      endif
c      write(57,'(4f10.5,f10.6,2f10.5)')qnu2,ef1,ef2,eg1,fcoup,
c     ,                                 dip,gmt/get

      end

      DOUBLE PRECISION FUNCTION sex3j(LA,LB,LC,NLF,NBT,NGA)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IF((NLF+NBT+NGA).NE.0.OR.IABS(LA-LB).GT.LC.OR.(LA+LB).LT.LC.
     .OR.IABS(NLF).GT.LA.OR.IABS(NBT).GT.LB.OR.IABS(NGA).GT.LC)THEN
      sex3j=0.0D00
      ELSE
      A1=ea((LA+NLF)/2)+ea((LB+NBT)/2)+ea((LB-NBT)/2)+ea((LC+NGA)/2)+
     +ea((LC-NGA)/2)+ea((LA-NLF)/2)
      D1=eelta(LA,LB,LC)/2.d0
      CNO=((-1)**(iabs(LA-LB-NGA)/2))
C
C    THIS IS DONE TO PREVENT OVERFLOW
C
      SUM=0.0D00
      N1=MAX0(LB-NLF-LC,LA+NBT-LC,0)/2
      N2=MIN0(LA+LB-LC,LA-NLF,LB+NBT)/2
      DO 100 N=N2,N1,-1
      B1=ea(N)+ea((LC-LB+NLF)/2+N)+ea((LC-LA-NBT)/2+N)+
     + ea((LA+LB-LC)/2-N)+ea((LA-NLF)/2-N)+ea((LB+NBT)/2-N)
      SUM=SUM+((-1)**n)*dexp(a1/2-b1)
100   CONTINUE
      sex3j=CNO*SUM*dexp(d1)
      ENDIF
      RETURN
      END
C
*************************** D E L T A ********************************
*--------------------------------------------------------------------*
*         AUXILIARY FUNCTION IN sex3j                                *
*--------------------------------------------------------------------*
C
      DOUBLE PRECISION FUNCTION eelta(LA,LB,LC)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IF((LA+LB+LC).LT.0)THEN
      eelta=0.0D00
      RETURN
      ELSE
      eelta=ea((LA+LB-LC)/2)+ea((LB+LC-LA)/2)+ea((LC+LA-LB)/2)-
     -ea((LA+LB+LC)/2+1)
      RETURN
      ENDIF
      END
C
C
C
************************** F A *****************************
C
      SUBROUTINE CONSea
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/eaCT/Aea(350)
      common/dbfc/dbf(350)
*------------------------------------------*
*     CALCULATES eaCTORIAL(N)              *
*------------------------------------------*
      FB=0.0D00
      DO 10 I=1,350
      FB=FB+dlog(Dfloat(I))
      Aea(I)=FB
 10   CONTINUE
      fd=1.d0
      do 20 j=1,250
      fd=1.d0
      do 30 k=1,j,2
      fd=fd*dfloat(k)
 30   continue
      dbf(j)=fd
 20   continue
      RETURN
      END
      double precision function dbeac(n)
      implicit double precision(a-h,o-z)
      common/dbfc/dbf(350)
      dbeac=dbf(n)
      return
      end
c
c
c
      DOUBLE PRECISION FUNCTION ea(N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/eaCT/Aea(350)
C     IF(N.EQ.-1)THEN
C     ea=1.E100
C     ELSE
      IF(N.EQ.0)THEN
      ea=0.0D0
      ELSE
      ea=Aea(N)
      ENDIF
C     ENDIF
      RETURN
      END


      subroutine jaffactor(inff,icharg,ireac,q,w,ef1,ef2,eg1)
      implicit real*8 (a-h,o-z)
      common/cons/amassp
	pi=dacos(-1.d0)
	hb=197.3d0
c	amass=amassp
	amass=939.d0/hb
	TAU=(Q*Q-W*W)/4.D0/AMASS/AMASS

	IF(IREAC.EQ.1.AND.ICHARG.EQ.1) THEN
	CALL NPIFF(inff,1,Q,W,GEPROT,GMPROT,GENEUT,GMNEUT)
	EF1p=(GEPROT+TAU*GMPROT)/(1.D0+TAU)
	EF2p=(GMPROT-GEPROT)/(1.D0+TAU)
	EF1n=(GENEUT+TAU*GMNEUT)/(1.D0+TAU)
	EF2n=(GMNEUT-GENEUT)/(1.D0+TAU)
	CALL WEAKFF(inff,1,Q,W,GEVPROT,GMVPROT,GEVNEUT,
     1	GMVNEUT,GAPROT,GANEUT)
	ef1=(ef1p-ef1n)*0.974
	ef2=(ef2p-ef2n)*0.974
	EG1=gaprot*0.974

	ENDIF
C
	IF(IREAC.EQ.1.AND.ICHARG.EQ.2) THEN
	CALL NPIFF(inff,1,Q,W,GEPROT,GMPROT,GENEUT,GMNEUT)
	EF1p=(GEPROT+TAU*GMPROT)/(1.D0+TAU)
	EF2p=(GMPROT-GEPROT)/(1.D0+TAU)
	EF1n=(GENEUT+TAU*GMNEUT)/(1.D0+TAU)
	EF2n=(GMNEUT-GENEUT)/(1.D0+TAU)
	CALL WEAKFF(inff,1,Q,W,GEVPROT,GMVPROT,GEVNEUT,
     1	GMVNEUT,GAPROT,GANEUT)
	ef1=-(ef1p-ef1n)*0.974
	ef2=-(ef2p-ef2n)*0.974
	EG1=-gaprot*0.974

	ENDIF

C
c{        write(57,'(8e10.3)')q,w,gevprot,gmvprot,gevneut,
c     2  gmvneut,gaprot,ganeut
	RETURN
	END
C+++++++++++++++++++++++++++++++++++++++++
      SUBROUTINE NPIFF(inff,icharg,Q,ENERGY,GEPROT,GMPROT,GENEUT,GMNEUT)
C+++++++++++++++++++++++++++++++++++++++++
C  E.M. FORM FACTORS OF THE PION AND THE NUCLEON
C
       IMPLICIT REAL*8 (A-H,O-Z)
      common/cons/amassp
	hb=197.3d0

	pi=dacos(-1.d0)

c	amass=amassp
	amass=939.d0/hb

C
	XROMEM=776.D0
	DFN=18.23D0
C
      OMEGA=ENERGY
      Q2=Q*Q-OMEGA*OMEGA
C
C------- NUCLEON ---------
C
      F1=1.0/(1.0 + Q2/DFN )**2
C
      IF(INFF.EQ.0) THEN
        GEPROT= F1
        GMPROT= F1*2.793
        GENEUT= 0.0
        GMNEUT= - F1*1.913
      ENDIF
C
C ---- JANSSENS ET AL. PHYS. REV. 142 (1966) 922 ----
C
      IF(INFF.EQ.1) THEN
        GES=0.5*( 2.5/(1.0+Q2/15.7) - 1.60/(1.0+Q2/26.7) + 0.10 )
        GEV=0.5*( 1.16/(1.0+Q2/8.19)  - 0.16)
        GMS=0.44*( 3.33/(1.0+Q2/15.7) - 2.77/(1.0+Q2/26.7) + 0.44 )
        GMV=2.353*( 1.11/(1.0+Q2/8.19) - 0.11)
C
        GEPROT= GES + GEV
        GMPROT= GMS + GMV
        GENEUT= GES - GEV
        GMNEUT= GMS - GMV
      ENDIF
C
C ---- BERTOZZI ET AL. PHYS. LETT. 41B (1972) 408 ---
C
      IF(INFF.EQ.2) THEN
        AA=0.058333
        BB=0.048333
        F1=1.0 + AA*Q2
        F2=1.0 + BB*Q2
C
        GEPROT= 1/F1/F1
        GENEUT= 1/F2/F2 - 1/F1/F1
        GMPROT= 2.793*GEPROT
        GMNEUT= -1.913*GEPROT
      ENDIF
C
C ---- IACHELLO ET AL. PHYS. LETT. 43B (1973) 191 ---
C
      IF(INFF.EQ.3) THEN
        T=-Q2*(HB*0.001)**2
        SMURHO=0.765
        SMUOME=0.784
        SMUPHI=1.019
        SMUNUC=0.939
        SMUPI=0.139
        GAMMA=0.48
        BRHO=0.746
        BOME=1.289
        BPHI=0.0
        APHI=0.0
        GT=0.5/(1.0-GAMMA*T)
        ALPHA=2.0*SQRT((-T+4.0*SMUPI**2)/(-T))/PI*
     ^  LOG( (SQRT(-T+4.0*SMUPI**2) + SQRT(-T))/(2.0*SMUPI) )
        GAMRHO=0.112
        H=(SMURHO**2+8.0*GAMRHO*SMUPI/PI)/
     ^   ( (SMURHO**2-T)+(4.0*SMUPI**2-T)*GAMRHO*ALPHA/SMUPI ) 
C
        F1S=GT*( (1.0 - BOME - BPHI) + BOME*(SMUOME**2/(SMUOME**2-T))
     ^            + BPHI*(SMUPHI**2/(SMUPHI**2-T)) )
        F1V=GT*( (1.0 - BRHO) + BRHO*H )
        F2S=GT*( (-0.12 - APHI)*( SMUOME**2/ (SMUOME**2-T) ) 
     ^            + APHI*( SMUPHI**2 / (SMUPHI**2-T)  ) )
        F2V=GT*( 3.706*H )
C
        GEPROT= F1V + F1S + (F2V + F2S)*T/(4.0*SMUNUC**2)
        GMPROT= F1V + F1S + F2V + F2S
        GMNEUT= -1.0*( F1V - F1S + F2V - F2S )
        GENEUT= F1S - F1V + (F2S - F2V)*T/(4.0*SMUNUC**2) 
      ENDIF
C
      IF(INFF.EQ.4) THEN
        T=-Q2*(HB*0.001)**2
        SMURHO=0.765
        SMUOME=0.784
        SMUPHI=1.019
        SMUNUC=0.939
        SMUPI=0.139
        GAMMA=0.65
        BRHO=0.588
        BOME=0.397
        BPHI=1.065
        APHI=0.166
        GT=0.5/(1.0-GAMMA*T)
        ALPHA=2.0*SQRT((-T+4.0*SMUPI**2)/(-T))/PI*
     ^  LOG( (SQRT(-T+4.0*SMUPI**2) + SQRT(-T))/(2.0*SMUPI) )
        GAMRHO=0.112
        H=(SMURHO**2+8.0*GAMRHO*SMUPI/PI)/
     ^   ( (SMURHO**2-T)+(4.0*SMUPI**2-T)*GAMRHO*ALPHA/SMUPI ) 
C
        F1S=GT*( (1.0 - BOME - BPHI) + BOME*(SMUOME**2/(SMUOME**2-T))
     ^            + BPHI*(SMUPHI**2/(SMUPHI**2-T)) )
        F1V=GT*( (1.0 - BRHO) + BRHO*H )
        F2S=GT*( (-0.12 - APHI)*( SMUOME**2/ (SMUOME**2-T) ) 
     ^            + APHI*( SMUPHI**2 / (SMUPHI**2-T)  ) )
        F2V=GT*( 3.706*H )
C
        GEPROT= F1V + F1S + (F2V + F2S)*T/(4.0*SMUNUC**2)
        GMPROT= F1V + F1S + F2V + F2S
        GMNEUT= -1.0*( F1V - F1S + F2V - F2S )
        GENEUT= F1S - F1V + (F2S - F2V)*T/(4.0*SMUNUC**2) 
      ENDIF
C
      IF(INFF.EQ.5) THEN
        T=-Q2*(HB*0.001)**2
        SMURHO=0.765
        SMUOME=0.784
        SMUPHI=1.019
        SMUNUC=0.939
        SMUPI=0.139
        GAMMA=0.25
        BRHO=0.687
        BOME=1.197
        BPHI=0.0
        APHI=0.0
        GT=0.5/(1.0-GAMMA*T)**2
        ALPHA=2.0*SQRT((-T+4.0*SMUPI**2)/(-T))/PI*
     ^  LOG( (SQRT(-T+4.0*SMUPI**2) + SQRT(-T))/(2.0*SMUPI) )
        GAMRHO=0.112
        H=(SMURHO**2+8.0*GAMRHO*SMUPI/PI)/
     ^   ( (SMURHO**2-T)+(4.0*SMUPI**2-T)*GAMRHO*ALPHA/SMUPI ) 
C
        F1S=GT*( (1.0 - BOME - BPHI) + BOME*(SMUOME**2/(SMUOME**2-T))
     ^            + BPHI*(SMUPHI**2/(SMUPHI**2-T)) )
        F1V=GT*( (1.0 - BRHO) + BRHO*H )
        F2S=GT*( (-0.12 - APHI)*( SMUOME**2/ (SMUOME**2-T) ) 
     ^            + APHI*( SMUPHI**2 / (SMUPHI**2-T)  ) )
        F2V=GT*( 3.706*H )
C
        GEPROT= F1V + F1S + (F2V + F2S)*T/(4.0*SMUNUC**2)
        GMPROT= F1V + F1S + F2V + F2S
        GMNEUT= -1.0*( F1V - F1S + F2V - F2S )
        GENEUT= F1S - F1V + (F2S - F2V)*T/(4.0*SMUNUC**2) 
      ENDIF
C
      IF(INFF.EQ.6) THEN
        T=-Q2*(HB*0.001)**2
        SMURHO=0.765
        SMUOME=0.784
        SMUPHI=1.019
        SMUNUC=0.939
        SMUPI=0.139
        GAMMA=0.25
        BRHO=0.672
        BOME=1.102
        BPHI=0.112
        APHI=-0.052
        GT=0.5/(1.0-GAMMA*T)**2
        ALPHA=2.0*SQRT((-T+4.0*SMUPI**2)/(-T))/PI*
     ^  LOG( (SQRT(-T+4.0*SMUPI**2) + SQRT(-T))/(2.0*SMUPI) )
        GAMRHO=0.112
        H=(SMURHO**2+8.0*GAMRHO*SMUPI/PI)/
     ^   ( (SMURHO**2-T)+(4.0*SMUPI**2-T)*GAMRHO*ALPHA/SMUPI ) 
C
        F1S=GT*( (1.0 - BOME - BPHI) + BOME*(SMUOME**2/(SMUOME**2-T))
     ^            + BPHI*(SMUPHI**2/(SMUPHI**2-T)) )
        F1V=GT*( (1.0 - BRHO) + BRHO*H )
        F2S=GT*( (-0.12 - APHI)*( SMUOME**2/ (SMUOME**2-T) ) 
     ^            + APHI*( SMUPHI**2 / (SMUPHI**2-T)  ) )
        F2V=GT*( 3.706*H )
C
        GEPROT= F1V + F1S + (F2V + F2S)*T/(4.0*SMUNUC**2)
        GMPROT= F1V + F1S + F2V + F2S
        GMNEUT= -1.0*( F1V - F1S + F2V - F2S )
        GENEUT= F1S - F1V + (F2S - F2V)*T/(4.0*SMUNUC**2) 
      ENDIF
C
C---- SIMON ET AL. NUCL. PHYS. A333 (1980) 381 ----
C
      IF(INFF.EQ.7) THEN
        GEPROT = 0.312/(1.0+Q2/6.0) + 1.312/(1.0+Q2/15.02)
     ^         - 0.709/(1.0+Q2/44.08) + 0.085/(1.0+Q2/154.2)
        GMPROT = 2.793*( 0.694/(1.0+Q2/8.5) + 0.719/(1.0+Q2/15.02) 
     ^          -0.418/(1.0+Q2/44.08) + 0.005/(1.0+Q2/355.4) )     
        GENEUT= 0.0
        GMNEUT= - 1.913/(1.0 + Q2/18.23 )**2
      ENDIF
C
C---- HOEL. ET AL. NUCL. PHYS. B114 (1976) 505 ----
C
      IF(INFF.EQ.8) THEN
       QT=Q2*(HB*0.001)**2
       XMNL=4.7580196
       XMS1= .7824
       XMS2=1.02
       XMS3=1.80
       XMV1=1.21
       XMV2=2.45
       XMV3=2.95
       A11S= .71
       A21S=-.64
       A31S=-.13
       A11V= .05
       A21V=-.52
       A31V= .28
       A12S=-.11
       A22S= .13
       A32S=-.02
       A12V=-1.99
       A22V=  .20
       A32V=  .19
       F1RHO=.5*(  .955+.090/(1.+QT/.355)**2 )/( 1.+QT/.536 )
       F2RHO=.5*( 5.335+.962/(1.+QT/.268)    )/( 1.+QT/.603 )
       F1S=A11S/(XMS1**2+QT)+
     ^     A21S/(XMS2**2+QT)+
     ^     A31S/(XMS3**2+QT)
       F2S=A12S/(XMS1**2+QT)+
     ^     A22S/(XMS2**2+QT)+
     ^     A32S/(XMS3**2+QT)
       F1V=F1RHO+A11V/(XMV1**2+QT)+
     ^           A21V/(XMV2**2+QT)+
     ^           A31V/(XMV3**2+QT)
       F2V=F2RHO+A12V/(XMV1**2+QT)+
     ^           A22V/(XMV2**2+QT)+
     ^           A32V/(XMV3**2+QT)
       F1S=2.*F1S
       F2S=2.*F2S
       F1V=2.*F1V
       F2V=2.*F2V
       GES=F1S-Q2*F2S/(4.*XMNL**2)
       GEV=F1V-Q2*F2V/(4.*XMNL**2)
       GMS=F1S+F2S
       GMV=F1V+F2V
       GEPROT=0.5*(GES+GEV)
       GMPROT=0.5*(GMS+GMV)
       GENEUT=0.5*(GES-GEV)
       GMNEUT=0.5*(GMS-GMV)
      ENDIF
C
C------- GALSTER ----------------
C
      IF(INFF.EQ.9) THEN
C
	AMED=AMASSP
	TAU=Q2/4.D0/AMED/AMED
	RMUP=2.79285D0
	RMUN=-1.91304D0
	RGVD=1.D0/((1.D0+4.97D0*TAU)**2.D0)
	RXIN=1.D0/(1.D0+5.6D0*TAU)
       GEPROT=RGVD
       GMPROT=RMUP*RGVD
       GENEUT=-RMUN*TAU*RGVD*RXIN
       GMNEUT=RMUN*RGVD
      ENDIF
c
c-------GARI AND KRUMPLEMANN [Z. PHYS. A322 (1985) 689
C
      IF(INFF.EQ.10) THEN
C
	AMED=AMASSP
	TAU=Q2/4.D0/AMED/AMED
	Q2=DABS(Q2)*(HB*0.001)**2
C
	RLAMBDA1=0.795D0
	RLAMBDA2=2.27D0
	RLQCD=0.29D0
C
	QH2=Q2*DLOG((RLAMBDA2**2+Q2)/RLQCD**2)/
     1	DLOG(RLAMBDA2**2/RLQCD**2)
C
C
	F1=RLAMBDA1**2.D0/(RLAMBDA1**2.D0+QH2)*RLAMBDA2**2.D0/
     1	(RLAMBDA2**2.D0+QH2)
	F2=RLAMBDA1**2.D0/(RLAMBDA1**2.D0+QH2)*(RLAMBDA2**2.D0/
     1	(RLAMBDA2**2.D0+QH2))**2.D0
C	------------------------------------------------------------
C	ASYMPTOTIC BEHAVIOUR
C	------------------------------------------------------------
c	FOR Q2 SMALL (Q2.LT.1.0 GEV)
c	F1=RLAMBDA1*RLAMBDA1/(RLAMBDA1**2.D0+Q2)
c	F2=F1
C	FOR Q2 LARGE
c	F1=1.D0/(Q2*DLOG(Q2/RLQCD**2.D0))**2.D0
c	F2=F1/Q2
C 	------------------------------------------------------------
C
	GWFW=0.411D0
	GRFR=0.377
	RKR=6.62
	RKW=0.163
C
	RKV=3.706
	RKS=-0.12
	RMR=0.776
	RMW=0.784
C
	COCR=RMR*RMR/(RMR*RMR+Q2)
	COCW=RMW*RMW/(RMW*RMW+Q2)
C
	F1IV=(COCR*GRFR+1.D0-GRFR)*F1
	F2IV=(COCR*RKR*GRFR+RKV-RKR*GRFR)*F2/RKV
	F1IS=(COCW*GWFW+1.D0-GWFW)*F1
	F2IS=(COCW*RKW*GWFW+RKS-RKW*GWFW)*F2/RKS
C
	F1P=0.5D0*(F1IS+F1IV)
	F1N=0.5D0*(F1IS-F1IV)
	F2P=0.5D0*(RKS*F1IS+RKV*F2IV)
	F2N=0.5D0*(RKS*F2IS-RKV*F2IV)
C
	GEPROT=F1P-TAU*F2P
	GMPROT=F1P+F2P
	GENEUT=F1N-TAU*F2N
	GMNEUT=F1N+F2N
	ENDIF
      RETURN
C
 1000 FORMAT(1X,' +++ DIPOLE FORM FACTOR +++')
 1001 FORMAT(1X,' +++ HOFSTADTER FORM FACTOR +++')
 1002 FORMAT(1X,' +++ IACHELLO FORM FACTOR #1 +++')
 1003 FORMAT(1X,' +++ IACHELLO FORM FACTOR #2 +++')
 1004 FORMAT(1X,' +++ IACHELLO FORM FACTOR #3 +++')
 1005 FORMAT(1X,' +++ IACHELLO FORM FACTOR #4 +++')
      END 
C


C+++++++++++++++++++++++++++++++++++++++++
      SUBROUTINE WEAKFF(inff,icharg,Q,ENERGY,GEVPROT,GMVPROT,
     1	GEVNEUT,GMVNEUT,GAPROT,GANEUT)
C+++++++++++++++++++++++++++++++++++++++++
C
C  WEAK FORM FACTORS OF THE NUCLEON
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      common/cons/amassp
	hb=197.3d0

	pi=dacos(-1.d0)

c	amass=amassp
	amass=939.d0/hb
C
      OMEGA=ENERGY
      Q2=Q*Q-OMEGA*OMEGA
	AM=AMASS
C
      QNU2=-Q2
      TAU = q2/(4.D0*AM*AM)


      sinw=0.233d0
      xiv0=-1.d0
      xivt0=-4.d0*sinw
      xivt1=2.d0-4.d0*sinw
      xia0=1.d0
      xiat0=0.d0
      xiat1=-2.d0

	CALL NPIFF(inff,icharg,Q,omega,GEPROT,GMPROT,GENEUT,GMNEUT)


      get1=0.5*(geprot-geneut)
      get0=0.5*(geprot+geneut)
      gmt1=0.5*(gmprot-gmneut)
      gmt0=0.5*(gmprot+gmneut)
      
c factor de forma de extragneza vectorial

      rhos=0.d0
      smu=-0.d0
      gvd=1.d0/(1+tau*4.97)**2
      slandae=5.6d0
      slandam=0.d0
      xies=1.d0/(1.0+slandae*tau)
      xims=1.d0/(1.0+slandam*tau)
      gest=rhos*tau*gvd*xies
      gmst=smu*gvd*xims
c	gest=0.d0
c	gmst=0.d0


      gevprot=xivt1*get1+xivt0*get0+xiv0*gest
      gmvprot=xivt1*gmt1+xivt0*gmt0+xiv0*gmst
      gevneut=-xivt1*get1+xivt0*get0+xiv0*gest
      gmvneut=-xivt1*gmt1+xivt0*gmt0+xiv0*gmst

c factores de forma axial isovector e isoescalar

      dd=1.262/1.64
      ff=0.64*dd
      gad=1.0/(1.0+tau*3.32)**2
      ga3=0.5*gad*(dd+ff)
      ga8=(3.d0*ff-dd)*gad*0.5/dsqrt(3.d0)
      
c factor de forma de extragneza axial

	slandas=0.d0
	xias=1.d0/(1.d0+slandas*tau)
      gas0=-0.d0
      gas=gas0*gad*xias

      gaprot=xiat1*ga3+xiat0*ga8+xia0*gas
      ganeut=-xiat1*ga3+xiat0*ga8+xia0*gas

      RETURN
      END


      subroutine wlimit(e,Tf,mi,M,B,EImax,w1,w2)
      implicit real*8 (a-h,o-z)

      real*8 M,mi,mi2,m2
c
c     computing limits for w transfer. 
c     
c     mi=mass of final lepton
c     M= mass of final nucleon
c     e= energy of incident lepton
c     B=binding energy (or any other energy loss)
c     EImax= maximum energy (M+T) posible for initial nucleon
c            (for instance at fermi level or adequate limit)
c     Tf=required kinetic energy of outgoing nucleon 
c
      mi2=mi*mi
      M2=M*M
      hc=197.33d0
c
c      absolute 'Q value' limits
c
       wabsL=B+TF+M-EIMAX
       wabsH=B+TF

       wi=max(wabsl,0.d0)
      
 
c       print*,'limites absolutos'
c       print*,wabsl*hc,wabsh*hc
c
c      kinematic 'Pi' and 'pf'>0 values
c

       wkinH1=e-mi
       wkinH2=TF+B     
 
c       print*,'limites cionematicos'
c       print*,wkinh1*hc,wkinh2*hc


       wkinh=min(wkinh1,wkinh2)
       wf=wkinh
       Ef=M+TF

       pf=dsqrt(ef*ef-M*M)




c      computing Pf




      sss1=ss1(e,ef,pf,mi,m,b,wkinh)
      sss2=ss2(e,ef,pf,mi,m,b,wkinh)
c      print*,'for wkinh=',wkinh*hc,sss1,sss2

      weipf=(ef+b)**2-M**2-(2*e-pf)**2
      weipf=weipf/2.d0/(pf+ef-e-e+B)

c      call sss(e,ef,pf,mi,m,b,weipf,ss1,ss2)

c      print*,'for weipf=',weipf*hc

      if(sss1.lt.0.and.sss2.gt.0)then
            write(68,*)'limit given by wkinh=',wkinh*hc
            wsol=wkinh
      endif
      if(sss1.gt.0)then
        write(68,*)'limit given by eq. --1--'
        call solve1(e,ef,pf,mi,m,b,wi,wf,wsol)
       write(68,*)'up limit  by eq. --1-- is ',wsol*hc
      endif
      if(sss2.lt.0)then
        write(68,*)'limit given by eq. --2--'     
        call solve2(e,ef,pf,mi,m,b,wi,wf,wsol)
        write(68,*)'up limit given by eq. --2-- is ',wsol*hc
      endif

      w1=max(wabsL,0.d0)
      w2=min(wabsH,wkinH1,wkinH2,wsol)

c      print*,'max =',w2*hc,tf*hc
c      print*,'min =',w1*hc,tf*hc
      return
      end

      subroutine solve1(e,ef,pf,mi,m,b,w1,w2,wsol)
      implicit real*8 (a-h,o-z)
      real*8 mi,m

      kstep=10
      hc=197.33d0

      wa=w1
      wb=w2
90    continue
      step=(wb-wa)/kstep     
      

      do 100 w=wa,wb-step*0.7,step

      wl=w
      wh=w+step

      a1=ss1(e,ef,pf,mi,m,b,wl)
      a2=ss1(e,ef,pf,mi,m,b,wh)

      if(a1*a2.lt.0.d0)then
           if(step*hc.lt.1.d-10)goto 200
           wa=wl
           wb=wh
           goto 90
      endif
100   continue
      print*,'cinematica no posile para ',(ef-m)*hc,' no hay solucion'
      wsol=-1e6
      return
200   continue
      wsol=(wa+wb)/2.d0
      return
      end

      subroutine solve2(e,ef,pf,mi,m,b,w1,w2,wsol)
      implicit real*8 (a-h,o-z)
      real*8 mi,m
      hc=197.33d0
      kstep=10

      wa=w1
      wb=w2
c      print*,w1*hc,w2*hc
90    continue
      step=(wb-wa)/kstep     
      
c      print*,step
      do 100 w=wa,wb-step+step/3.d0,step

      wl=w
      wh=w+step

      a1=ss2(e,ef,pf,mi,m,b,wl)
      a2=ss2(e,ef,pf,mi,m,b,wh)
 
c      write(*,'(5e15.8)')step*hc,wl*hc,wh*hc,a1,a2

      if(a1*a2.lt.0.d0)then
           if(step*hc.lt.1.d-10)goto 200
           wa=wl
           wb=wh
           goto 90
      endif
100   continue
      print*,'cinematica no posible para ',(ef-m)*hc,' no hay solucion'
      wsol=-1e6
      return
200   continue
      wsol=(wa+wb)/2.d0
      return
      end

      real*8 function ss1(e,ef,pf,mi,M,B,w)
      implicit real*8 (a-h,o-z)
      real*8 mi,m
      ss1=e-dsqrt(dabs((e-w)**2-mi**2))-pf-dsqrt(dabs((ef-w+b)**2-m**2))
c      ss2=e+dsqrt((e-w)**2-mi**2)-pf+dsqrt((ef-w+b)**2-m**2)
c      print*,ss1,ss2

      return
      end
      real*8 function ss2(e,ef,pf,mi,M,B,w)
      implicit real*8 (a-h,o-z)
      real*8 mi,m
      ss2=e+dsqrt(dabs((e-w)**2-mi**2))-pf+dsqrt(dabs((ef-w+b)**2-m**2))
c      print*,ss1,ss2

      return
      end


