csinpmc234567

      program pwba
      implicit real*8 (a-h,o-z)
      include 'param.inc'
      common/num/step,rfin,kmax
      common/uno/hp(0:100)
      common/pfin/akpf,epf
      common/cons/amassp
      common/harm/ysh(0:npw+1,-npw:npw)
c
      real*8 yspm(0:npw+1,-npw:npw),
     ,xx(npunt1),fff(3),ux(npunt1),wx(npunt1),uy(npunt1),
     ,wy(npunt1),ub(npunt1),wb(npunt1)
      real*8 step,r,hst,fac,s2(-1:1)
      character*20 rpiv
      character*60 ipivot
      character*26 infile,outfile,fil1,fil2,fil3,fil4,fil5,fil6
      complex*16 zh,zi,zero,zj1x,zj1y,zj1z,zj2x,zj2y,zj2z,zrho,
     , zjplus,zjmins,zj0,zeroz,zjpar,zjper,zlt,zjz,
     , pfmas,pfmin,zpm(-1:1),zpf(-1:1),zwx,zwy,zwl,
     , zs3(-1:1),zs2,zalp(-1:4),zbet(-1:4),zgam(-1:4),zs1,zesp1,zesp2,
     , ztran1,ztran2,ztot1,ztot2,zw0,zw1,zw2,zws1,zwt1,zwt2,zwt4,
     , s0,s0b,v(-1:1),vb(-1:1),sv,svb,vv(-1:1),vvb(-1:1),zcros(-1:1),
     , zcrosb(-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

      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
c      print*,'mete escale'
c      read(*,*)esc
c      hc=hc/esc
      amassp=938.d0/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)


      open(1,file='andwia.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=outfile)
      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.913d0

       if(icharg.eq.1)      anom=1.793d0


      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





      if(lb2.gt.jb2)then

       kb2=lb2

       lb2bar=lb2-2

       

      else

       kb2=-lb2-2

       lb2bar=lb2+2

      endif

      kb=kb2/2

      skb=kb/iabs(kb)

      lbar=lb2bar/2

      l=lb2/2

c

c

c

c

      write(13,*)nb,lb2,jb2,spft,iff,ireac,icharg,asig
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,*)Tpf,pmin,pmfin,dpm
      write(13,9991)tpf,pmin,pmfin,dpm
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,*)eki,Emiss,qcons
      write(13,9991)eki,Emiss,qcons
      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,n2,npmx,nemx,lmx,imp,ida
      write(13,1111)step,rfin,n2,npmx,nemx,lmx,imp,ida
 1111 format(2e13.6,7i5.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 
      write(13,*)ipl,ipt,icc
c
c

c     Fin de la entrada de datos

c

      kred=3

      lmx2=lmx*2

      kemx2=nemx*2

      kpmx=npmx

      kpmx2=kpmx*2

c

c     vamos a leer la funcion ligada

c

      fcs=step/3.d0

      kmax=int(rfin/step+0.5d0)

      kmax=2*((kmax+1)/2)

      rfin=kmax*step

      print*,' new value for Rfin=',rfin

c

c      open(15,file=fil3,status='new')

      open(24,file='ancrs3.out')
      open(25,file='ancrs4.out')
      open(26,file='sigcc1.out')
      open(11,file=infile,status='old')
      open(37,file='fact.out')
      open(38,file='nofact.out')
      open(67,file='fact1.out')
      open(68,file='nofact1.out')
      open(39,file='ancomp.out')
      open(40,file='anmpl4.out')
      open(57,file='formf.out')
      open(46,file='response.out')
      open(47,file='ncurrent.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,*)'-----------------------------------------------------'





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

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

      read(11,'(60a)')ipivot
      read(11,*)nlines,it2,lb22,jb22

      r=0.0d0
      fac=4.d0
      sum=0.0d0
      sum2=0.0d0
      nlin=1

      do 99 i=1,nlines

      if(xx(nlin).gt.rfin)goto 99

      nlin=i

      read(11,*)xx(i),ux(i),wx(i)
      fx1=ux(i)*ux(i)+wx(i)*wx(i)
      fx1=fx1*xx(i)*xx(i)
      sum2=sum2+fx1*xx(i)*xx(i)*fac
      sum=sum+fx1*fac

      fac=6.d0-fac

  99  continue

      nlines=nlin

      print*,'nlines=',nlines

      anorm=(sum-fx1)*(xx(2)-xx(1))/3.d0

c      if(imp.eq.1)then

        print*,'norma=',anorm

        print*,'rbound=',xx(nlines)

        write(13,*)'norma=',anorm

        write(13,*)'rbound=',xx(nlines)

c      endif





 451  continue

 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'

c
c     masa nucleo residual y blanco
c

      amassb=z2mass*ua
      amassa=z1mass*ua

c
c     energia total, momento y energia cinetica del proton final
c
 
      Tpf=tpf/hc

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

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

      Emiss=Emiss/hc



      amassb=amassb+emiss

c
c     energia electron incidente
c

      eki=eki/hc
      eki2=eki*eki

c

c
c     calcular q y demas (la cinematica)
c

c
c     step en ekf
c

      sekf=dpm/hc


      sigaekf=0.d0
      ichekf=0
      facekf=4.0d0
      if(amassa.lt.(akpf+2*eki+epf))then

        print*,' AMASS is TOO LOW!!!!'
        write(37,*),' AmASS is TOO LOW!!!'
        write(38,*),' AMASS is TOO LOW!!!'
    
      endif

      wmin=amassb-amassa+epf+emiss
      wpm=dsqrt(akpf**2+amassb**2)-amassa+epf
      weki=2*eki-akpf
      wplus=(amassa-epf)**2-amassb**2-(weki)**2
      wplus=wplus/2/(akpf-2*eki-amassa+epf)
      wmins=(akpf+2*eki)**2+amassb**2-(amassa-epf)**2
      wmins=wmins/2/(amassa+2*eki+akpf-epf)
      wpf=akpf
      wpfpm=amassb**2-(amassa-epf)**2+akpf**2
      wpfpm=wpfpm/2/(amassa-epf+akpf)


      print*,'wpm=',wpm*hc
      print*,'weki=',weki*hc
      print*,'wplus=',wplus*hc
      print*,'wmins=',wmins*hc
      print*,'wpfpm=',wpfpm*hc
      print*,'wmin=',wmin*hc

      iw=0

      wlow=max(wplus,weki,wmin)

      if(wpm.gt.wlow) then
         whigh=wpm
        iw=1
      print*,' case 1 ',wlow*hc,whigh*hc
      endif


      wl0=min(wpm,weki,wmins)

      if (wl0.gt.wmin) then


        if (wl0.lt.akpf) then

          wlow=wmin
          whigh=wl0
          iw=2
          print*,' case 2' 
        else

         print*, ' We didn''t think that this case was needed !!!  '

        endif

      endif
 

      if(iw.ne.0)then
      wlim=whigh
      wlowlow=wlow
      else
      wlim=0
      endif 

      if(max(wpm,wmin).lt.min(wmins,akpf))then

      wlow=max(wpm,wmin)
      whigh=min(wmins,akpf)
      iw=3
      print*,' case 3 ',wlow*hc,whigh*hc
      endif

      if(wlim.ge.wlow)then
      wlow=wlowlow
      iw=6
      endif
      print*,' case ',iw,wlow*hc,whigh*hc

      ekfmax=eki-wlow

      write(*,*)' eki, akpf',eki*hc,akpf*hc

      write(*,*)' ekf max=',ekfmax*hc



      qmax=eki+ekfmax

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

      write(*,*),' pmmin=',(akpf-qmaxf)*hc


      ekfmin=eki-whigh


      sekf=(ekfmax-ekfmin)/100.d0
c      sekf=0.001/hc
      ekfmin=ekfmin-sekf/3.

c      ekfmax=ekfmax*0.999999d0

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



      do 1 ekf=ekfmax,ekfmin,-sekf



c      print*,ekf,sekf



      tick=0.1d0

      nstcos=4

      nstcosold=4

c

c     pm=pmissing=-pb(nucleo residual)=p proton inicial (I.A.)

c



      w=eki-ekf



      Eb=w-epf+amassa

      print*,eb,amassb,eb-amassb



      akr=dsqrt(eb*eb-amassb*amassb)

      print*,eb*eb,amassb*amassb

      pm=akr



      pm2=pm*pm







      Tb=Eb-amassb



      Erec=Tb



c

c     calculo de w y Epi

c



      w2=w*w



      Epi=Epf-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     



      Epibar=dsqrt(pm2+amassp*amassp)



      wbar=(Epf-Epibar)



      w2bar=wbar*wbar



c

c     energia electron final

c



      ekf2=ekf*ekf





      print*,'epi=',epi*hc

      print*,'akr=',akr*hc





      aminfq=dabs(akr-akpf)



      amaxfq=akr+akpf











      aminiq=w



      amaxiq=eki+ekf



      qmin=max(aminiq,aminfq)



      qmax=min(amaxfq,amaxiq)



      sinphi=0.d0



      cosphi=1.d0



c      write(*,2393)eki*hc,ekf*hc,sekf*hc,epf*hc,

c     ,w*hc,pm*hc,eb*hc,tb*hc,epi*hc,wbar*hc,qmin*hc,

c     , qmax*hc,emiss*hc,cosin,cosfn,stcos







      if(qmax.le.qmin) then

      print*,'qmax,qmin ',qmax,qmin

      print*,aminiq,aminfq,amaxiq,amaxfq

c      goto 1

      endif





      cosin=(qmin*qmin-eki*eki-ekf*ekf)/2.d0/eki/ekf



      cosfn=(qmax*qmax-eki*eki-ekf*ekf)/2.d0/eki/ekf



 

c      print*,cosin,cosfn





c-----------------------------------------------

c    tick is the maximun step in dcos, so the step

c    is the minor of tick and (cosf-cosi)/nstcos

c    where nstcos is the minimum number of steps

c-------------------------------------------------     







      stcos=(cosfn-cosin)/nstcos



      if(stcos.eq.0.d0) goto 1 



      if(stcos.gt.tick)then



      astcos=stcos/tick

      nstcos=idint(astcos+0.5)

      nstcos=2*(nstcos/2+1)

      nstcos=min(nstcos,nstcosold)

      stcos=stcos/nstcos



      endif





c      print*,stcos

c







c      print*,'cosenos y steps',cosin,cosfn,stcos

c      print*,'qs y demas',qmin,qmax,akr





      facos=4.0d0

      sigacos=0.0d0

      ichcos=0







      write(*,2393)eki*hc,ekf*hc,sekf*hc,epf*hc,

     ,w*hc,pm*hc,eb*hc,tb*hc,epi*hc,wbar*hc,qmin*hc,

     , qmax*hc,emiss*hc,cosin,cosfn,stcos





2393  format(/' e lept i=',f9.3,' e lept f=',f9.3,'  e step=',f9.3,

     , '   MeV'/

     ,       '  Epf     =',f9.3,'   w     =',f9.3,'  pm    =',f9.3,

     , '   MeV'/

     ,       '  Eb      =',f9.3,'   Tb    =',f9.3,'  Epi   =',f9.3,

     , '   MeV'/

     ,       '  Wbar    =',f9.3,'  qmin   =',f9.3,'  qmax  =',f9.3,

     , '   MeV'/

     ,      ' Emiss     =',f9.3,'  cosin  =',f9.6,'  cosfn =',f9.6,

     ,      '  stcos    =',f9.6 /)



      cosfn=cosfn+stcos/2.d0

 





      do 1300 cosnu=cosin,cosfn,stcos



      q2=eki*eki+ekf*ekf+2.d0*eki*ekf*cosnu



c      print*,q2,w2



      cosif=cosnu



      q=dsqrt(q2)

      cosf=q2+akpf*akpf-akr*akr

      cosf=cosf/2.d0/q/akpf      

      cosm=akpf/akr*cosf-q/akr















      iysh=1





      cosgam=cosf







      cos2gam=cosgam*cosgam

      singam=dsqrt(dabs(1.d0-cos2gam))      







      pmx=akpf*singam*cosphi

      pmy=akpf*singam*sinphi

      pmz=akpf*cosgam-q







      pfx=pmx

      pfy=pmy

      pfz=akpf*cosgam











      q2=q*q

      Epi2=Epi*Epi

      Epf2=Epf*Epf

      cosif=cosnu



      qnu2=w2-q2



      qnu2bar=w2bar-q2



      sin2if2=-qnu2/4.d0/eki/ekf



      cos2if2=1.d0-sin2if2



      if(cos2if2.lt.0.and.dabs(cos2if2).lt.1.d-10)then

      cos2if2=-cos2if2

      endif



      tg2if2=sin2if2/cos2if2







      tgif2=dsqrt((tg2if2))







c

c     impresion de las variable cinematicas del electron y el proton...

c



c      thet=180.d0/pi*dacos(cosif)

c      gamma=180.d0/pi*dasin(singam)

c      phi=180.d0/pi*dacos(cosphi)

c      costh=qnu2/2/eki/ekf+1.d0

c      th=dacos(costh)



c      ekix=eki*ekf*dsin(th)/q

c      ekiz=(eki2-eki*ekf*dcos(th))/q

c      ekfx=ekix

c      ekfz=(-ekf2+eki*ekf*dcos(th))/q

c      ekiy=0.d0

c      ekfy=0.d0







c      thi=datan2(ekix,ekiz)

c      thf=datan2(ekfx,ekfz)















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



 





      if (qnu2.ne.0.d0)then



      sigmott=4.d0*alpha*alpha*ekf2*cos2if2/qnu2/qnu2

      sm=sigmott

      else



      sigmott=0

      sm=0

      endif

      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





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



      write(26,'(3f9.4,e13.5)')f1,f2,f3,sigmott

      write(26,3664)pm*hc,gamma,phi,thet,sigcc1,scc1ro,

     , 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

c      write(*,*)tg2if2,qnu2/q2,qnu2,q2,sin2if2,cos2if2

      Slt=-dsqrt((tg2if2-qnu2/q2))*2.d0

      Sl0=w2/q2

      Sl1=w/q*2.d0

      Att=asig*2.d0*dsqrt((tg2if2-qnu2/q2))*tgif2 

      Stt1=qnu2/q2

      Stt2=stt1/2.d0









1300  continue





      facekf=6.d0-facekf



      if(ichekf.eq.0)then

      sigekfi=sigcos

      endif 



      write(38,'(3f10.5,e13.5)')ekf*hc,q*hc,pm*hc,sigcos

      write(38,*)stcos,sig,factor,waverage



      

    1 continue





      sigtol=(sigaekf-sigcos)*sekf/3.d0





      write(37,'(f10.5,e13.5)')tpf*hc,sigtol





      close(37)

      close(38)

      close(46)

      close(47)

      close(13)

      close(57)

      close(67)

      close(68)

      close(11)

      close(13)

      close(26)

      close(24)

      close(25)

      close(1)

      end



