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.3d0/hc
c      amassp=0.1d0
      alph=137.d0

      cpmfac=1.d0

      ua=931.4944/hc
      sqr2=dsqrt(2.d0)
      sqr6=dsqrt(6.d0)
      fourpi=4.d0*pi
c      call cputim(it1)

      call consea

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

       if(icharg.eq.2)   then
c      neutron
    
       amassp=939.6/hc     
       anom=-1.91304 

       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





      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,*)Tpfin,tpffin,tpfstep,pmin,pmfin,dpm
      write(13,9991)tpfin,tpffin,tpfstep,pmin,pmfin,dpm
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,*)eki,Emiss,tick
      write(13,9991)eki,Emiss,tick
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,123)infile
      write(13,123)infile
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,1234)ipivot
      write(13,*)ipivot
      read(1,*)step,rfin,nstcos,cosmin,cosmax,imp
      write(13,1111)step,rfin,nstcos,cosmin,cosmax,imp
 1111 format(2e13.6,i5.3,2e13.5,i5.3)

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

c     Fin de la entrada de datos

c

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

c      open(24,file='ancrs3.out')
c      open(25,file='ancrs4.out')
      open(26,file='sigcc1.out')
      open(11,file=infile,status='old')
      open(37,file='2int.out')
      open(38,file='0int.out')
c      open(67,file='fact1.out')
c      open(68,file='nofact1.out')
      open(39,file='1int.out')
c      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
      sumu=0.d0
      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)
      fxu=ux(i)*ux(i)*xx(i)*xx(i)
      fx1=fx1*xx(i)*xx(i)
      sum2=sum2+fx1*xx(i)*xx(i)*fac
      sum=sum+fx1*fac
      sumu=sumu+fxu*fac

      fac=6.d0-fac

  99  continue

      nlines=nlin

      print*,'nlines=',nlines

      anorm=(sum-fx1)*(xx(2)-xx(1))/3.d0
      unorm=(sumu-fxu)*(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,*)'unorm=',unorm

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

c      endif

      write(13,*)'RMS RADIUS=',
     ,  dsqrt((sum2-fx1*xx(i)*xx(i))*(xx(2)-xx(1))/3.d0),anorm



      call spline(xx,ux,nlines,0.d0,0.d0,uy)
      call spline(xx,wx,nlines,0.d0,0.d0,wy)

      r=0.d0 

      do 451 i=1,kmax

      r=r+step

      call splint(xx,ux,uy,nlines,r,ub(i))
      call splint(xx,wx,wy,nlines,r,wb(i))

      if(ipt.eq.4)      wb(i)=-wb(i)
      if(ipt.eq.2)      wb(i)=0.d0
      if(ipt.eq.1)      ub(i)=0.d0

      if(imp.eq.1)print978,i,r,ub(i),wb(i)

 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
 

      emiss=emiss/hc
      amassb=amassb+emiss
      eki=eki/hc

      sekf=dpm/hc
      pmfmin=pmfin/hc

      eki2=eki*eki

      do 12345 tpfdo=tpfin,tpffin,tpfstep

      Tpf=tpfdo/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





c
c     energia electron incidente
c


c

c
c     calcular q y demas (la cinematica)
c

c
c     step en ekf
c




      sigaekf=0.d0
      ichekf=0
      facekf=2.0d0
      sigaekftt=0.0d0
      sigaekfts=0.0d0
      sigaekfst=0.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
      wlow=wmin
      wext=dsqrt(pmfmin**2+amassb**2)+epf-amassa

      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
      print*,'wext=',wext*hc,pmfmin*hc

      iw=0

      wlow=max(wplus,weki,wmin)

         wh1=wpm
         wl1=wlow
      print*,' case 1 ',wl1*hc,wh1*hc


      wl0=min(wpm,weki,wmins)

      if (wl0.gt.wmin) then


        if (wl0.lt.akpf) then

          wl2=wmin
          wh2=wl0
          print*,' case 2',wl2*hc,wh2*hc 
        else

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

        endif

      endif
 

      wl2=wmin
      wh2=wl0
      print*,' case 2 no if',wl2*hc,wh2*hc


      i1=0
      i2=0
      i3=0


      if(wl1.lt.wh1)then
      i1=1
      print*,'case 1 is possible'
      endif


      if(wl2.lt.wh2)then
      i2=1
      print*,'case 2 is possible'
      endif


      wl3=max(wpm,wmin)
      wh3=min(wmins,akpf)
      print*,' case 3 ',wl3*hc,wh3*hc



      if(wl3.lt.wh3)then
      i3=1
      print*,'case 3 is possible'
      endif

     
      wlow=min(wl1+999*(1-i1),wl2+999*(1-i2),wl3+999*(1-i3))
      whigh=max(wh1*i1,wh2*i2,wh3*i3)


      print*,' case ',i1,i2,i3,wlow*hc,whigh*hc




      whigh=min(whigh,wext)

      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


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





      isimps=0

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

      ekfmin=1.00000001*ekfmin
      ekfmax=ekfmax*0.9999999d0

      pmin=dsqrt((amassa+eki-ekfmax-epf)**2-amassb**2)
      pmax=dsqrt((amassa+eki-ekfmin-epf)**2-amassb**2)

      pscal=eki
      amin=dsqrt(pmin/pscal)
      amax=dsqrt(pmax/pscal)

      sekf=(amax-amin)/dpm
      nsimps=(amax-amin)/sekf
      nsimps=2*(nsimps/2) 
      sekf=(amax-amin)/nsimps

      amax=amax+sekf/30.d0

      print*,'amin, amax, sa',amin,amax,sekf


      do 1 alp=amin,amax,sekf

      ekf=eki-epf+amassa-dsqrt(alp**4*pscal**2+amassb**2)


      isimps=isimps+1


c      tick=0.05d0
c      nstcos=4
c      nstcosold=4
c
c     pm=pmissing=-pb(nucleo residual)=p proton inicial (I.A.)
c

      w=eki-ekf

      Eb=w-epf+amassa

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

      pm=akr

      pm2=pm*pm


      ajacob=2.d0*alp**3*pscal**2/dsqrt(amassb**2+pm**2)

      print*,alp,ekf*hc,ajacob



      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



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

c
c     energia electron final
c

      ekf2=ekf*ekf





      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



      if(qmax.lt.qmin) goto 1



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

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

 
      cosin=min(cos1,cos2)
      cosfn=max(cos1,cos2)


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

      cosfn=min(cosmax,cosfn)
      cosin=max(cosmin,cosin)

      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)
c      nstcos=max(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=2.0d0
      sigacos=0.0d0
      ichcos=0
      sigacost=0.0d0


      write(*,2393)eki*hc,ekf*hc,akpf*hc,epf*hc,
     ,w*hc,pm*hc,epibar*hc,tb*hc,epi*hc,wbar*hc,qmin*hc,
     , qmax*hc,emiss*hc,cosin,cosfn,stcos
      ekfold=ekf

2393  format(/' e lept i=',f9.3,' e lept f=',f9.3,' akpf   =',f9.3,
     , '   MeV'/
     ,       '  Epf     =',f9.3,'   w     =',f9.3,'  pm    =',f9.3,
     , '   MeV'/
     ,       '  Epibar  =',f9.3,'   Tb    =',f9.3,'  Epi   =',f9.3,
     , '   MeV'/
     ,       '  Wbar    =',f9.3,'  qmin   =',f9.3,'  qmax  =',f9.3,
     , '   MeV'/
     ,      ' Emiss     =',f9.3,'  cosin  =',f9.6,'  cosfn =',f9.6,
     ,      '  stcos    =',f9.6 /)

c
c       las integrales analiticas....
c

      fac=4.d0
      l=lb2/2
      lbar=lb2bar/2
      r=0.d0
      sum1=0.d0
      sum2=0.d0



      apm=dabs(pm)



      do 551 i=1,kmax
      r=r+step
      xpm=r*apm



      hpl=djl(l,xpm)
      hplbar=djl(lbar,xpm)



c      call hl(xpm)





      DENS1=UB(I)*hpl*r*r
      dens2=WB(I)*hplbar*r*r
      SUM1=SUM1+DENS1*FAc
      sum2=sum2+dens2*fac
      FAC=6.D0-FAC

 551  CONTINUE

      SUM1=(SUM1-DENS1)*FCS
      sum2=(sum2-dens2)*fcs
      ubpm=sum1
      wbpm=-skb*sum2

c      print*,'integrales',pm*hc,ubpm,wbpm



c
c     un caso de test, wbpm=-pm/(Em+amassp)*ubpm
c

      epibar=dsqrt(pm*pm+amassp*amassp)

c      print*,'--------------- TEST---------------'

      if(ipt.eq.3)then

       wbpm=-dabs(pm)/(Epibar+amassp)*ubpm
       cpmfac=(epibar+amassp)/2.d0/epibar/unorm

      endif




      write(*,'(2i4.2,4e13.5)')l,lbar,pm*hc,
     , ubpm,wbpm,-dabs(pm)/(Epibar+amassp)*ubpm



c      cosin=.999999999*cosin
c      cosfn=.999999999*cosfn


      cosfn=cosfn+stcos/1.d10


 
      print*,cosin,cosfn

      incos=0

      do 1300 cosnu=cosin,cosfn,stcos

      incos=incos+1



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

c      print*,dsqrt(q2)*hc,w2,cosnu,ekf*hc

      cosif=cosnu

      q=dsqrt(q2)

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






      iysh=1

c

c     un caso de test

c    

      cosgam=cosf

c      print*,cosgam

      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



c

c
c
c
c
c     ahora se calculan los angulos respecto a los ejes en cuestion
c
c





c
c     construccion del Ylm(Pm) 
c      
c

c       calculamos los armonicos esfericos que necesitaremos despues
c       y los almacenamos, hay que colocar esta rutina 'antes' de los
c       otros armonicos.
c

      if(dabs(pm).gt.1.d-12)then

       sinpm=pmx/dabs(pm)
       cospm=pmz/dabs(pm)

       iysh=1

       if(sinpm.lt.0)then

         sinpm=dabs(sinpm)

         iysh=-1

       endif

       call harmsq(sinpm,cospm)

       do 7444 lll=0,npw

        do 7444 m=-npw,npw

         yspm(lll,m)=(iysh**(iabs(m)))*ysh(lll,m)

7444   continue

      else

       iysh=1

       if(pm.lt.0.d0)iysh=-1

       sinpm=0.d0

       cospm=1.d0

       do 7445 lll=0,npw

         yspm(lll,0)=iysh**lll*dsqrt(dfloat(lll+lll+1)/fourpi)

         do 7445 m=1,npw

          yspm(lll,m)=0.d0
          yspm(lll,-m)=0.d0

7445   continue

      endif

c      write(47,*)'thpm=',dasin(sinpm)*180.d0/pi
c      write(47,*)sinpm,cospm,l,lbar




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
      Epi2=Epi*Epi
      Epf2=Epf*Epf
      cosif=cosnu

      qnu2=w2-q2

      qnu2bar=w2bar-q2

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

      cos2if2=1.d0-sin2if2
c      if(cos2if2.eq.0.d0)cos2if2=1.d-15

c      print*,qnu2,4*eki*ekf,q2,w2,cosif,sin2if2,cos2if2

      tg2if2=sin2if2/cos2if2

c      print*,sin2if2,w2,q2,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      write(13,818)pm*hc,q*hc,eki*hc,ekf*hc,w*hc,
c     , thet,(Epi-amassp)*hc,(Epibar-amassp)*hc,akpf*hc,
c     , gamma,Emiss*hc

c      print818,pm*hc,q*hc,eki*hc,ekf*hc,w*hc,
c     , thet,(Epi-amassp)*hc,(Epibar-amassp)*hc,akpf*hc,
c     , gamma,Emiss*hc

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

 
      call ffactor(idip,icharg,ireac,qnu2,fff(1),eff,ggg,fcoup,f2s,g1s)

      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
c




C     MOTT'S CROSS SECTION
      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

      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

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



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



      ind=0



      Rlac=0.d0
      Rtac=0.d0
      Rttac=0.d0
      Rltac=0.d0
      Rattac=0.d0






      do 5999 mb2=-jb2,jb2,2

      do 5999 ispf2=-1,1,2

      ind=ind+1

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



      s0=cgc(lb2,mb2-ispf2,1,ispf2,jb2,mb2)*
     * yspm(l,(mb2-ispf2)/2)

      s0b=cgc(lb2bar,mb2-ispf2,1,ispf2,jb2,mb2)*
     * yspm(lbar,(mb2-ispf2)/2)


      sv=zeroz

      svb=zeroz


      do 7555 im2=-2,2,2

      m12=mb2-ispf2+im2

      f1b=cgc(lb2bar,m12,1,ispf2-im2,jb2,mb2)*yspm(lbar,m12/2)

      f1=cgc(lb2,m12,1,ispf2-im2,jb2,mb2)*  yspm(l,m12/2)






      f0=((-1)**(iabs(1-ispf2)/2))*sqr6*
     *     sex3j(2,1,1,im2,ispf2-im2,-ispf2)



      v(im2/2)= f0 * f1

      vb(im2/2)= f0* f1b



      sv=sv+zpf(-im2/2)*v(im2/2)*(-1)**(iabs(im2/2))

      svb=svb+zpf(-im2/2)*vb(im2/2)*(-1)**(iabs(im2/2))





7555  continue







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



      zcrosb(0)=  zi*(vb(-1)*zpf(1)-vb(1)*zpf(-1))
      zcrosb(1)=  zi*(vb(0)*zpf(1) -vb(1) *zpf(0))
      zcrosb(-1)= zi*(vb(-1)*zpf(0)-vb(0) *zpf(-1))





      do 7556 im2=-1,1,1

      vv(im2)=s0*zpf(im2)+zi*zcros(im2)

      vvb(im2)=s0b*zpf(im2)+zi*zcrosb(im2)

7556  continue

c
c     ahora construimos las funciones analiticas
c

c     primero unos factores

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

c

      srf=1.d0
      srb=1.d0/(Epf+amassp)

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

c     ahora los elementos matriciales

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

c

      zsuu=ubpm*s0
      zsdd=wbpm*srb*svb
      zsdu=ubpm*srb*sv
      zsud=wbpm*s0b




      do 7557 iv=-1,1,1

      zvdu(iv)=srb*ubpm*vv(iv)
      zvud(iv)=wbpm*vb(iv)
      zvuu(iv)=ubpm*v(iv)
      zvdd(iv)=srb*wbpm*vvb(iv)

7557  continue



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

       

       zalp(4)=(-zi)**(lb2/2)*srf*(zsuu+zsdd)

c
c      espacial  gamma_i
c

   

       do 558 im2=-2,2,2 

       zalp(im2/2)=(-zi)**(lb2/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)=-(-zi)**(lb2/2)*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)=(-zi)**(lb2/2)*srf*zbet(0)
      zbet(1)=(-zi)**(lb2/2)*srf*zbet(1)
      zbet(-1)=(-zi)**(lb2/2)*srf*zbet(-1)


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


      zgam(4)=(-zi)**(lb2/2)*srf*(zsuu-zsdd)*(Epibar+Epf)


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



      do 6666 i=-1,1

      zgam(i)=(-zi)**(lb2/2)*srf*(zsuu-zsdd)*(zpm(i)+zpf(i))



6666  continue


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



      zfive(4)=(-zi)**(lb2/2)*srf*(zsud+zsdu)



      do 9985 iv=-1,1,1

      zfive(iv)=(-zi)**(lb2/2)*srf*(zvuu(iv)+zvdd(iv))

9985  continue



      zrho=falp*zalp(4)+fbet*zbet(4)+fgam*zgam(4)+g1*zfive(4)

  







      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)))


      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)))   


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



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

c      write(47,'(f6.1,3i4.2,6f10.4)')pm*hc,
c     ,          ispf2,mb2,4,zalp(4),zbet(4),zfive(4)

c      do 1047 imq=-1,1,1
c      write(47,'(14x,i4.2,6f10.4)')imq,zalp(imq),zbet(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------------------------------------------------------

      Rl=dble(zrho*dconjg(zrho))+w2/q2*dble(zj1z*dconjg(zj1z))-
     -                           2.d0*w/q*dble(dconjg(zrho)*zj1z)

      Rt=dble(zj1x*dconjg(zj1x))+dble(zj1y*dconjg(zj1y))

      Rlt=dble(dconjg(zrho-w/q*zj1z)*(-sinphi*zj1y+cosphi*zj1x)) 

      Rtt=Stt1*sinphi2*dble(zj1y*dconjg(zj1x))-
     -    stt2*cosphi2*(dble(zj1x*dconjg(zj1x))-dble(zj1y*dconjg(zj1y)))

      Ratt=dimag(zj1x*dconjg(zj1y))

      Rlac=Rlac+Rl

      Rtac=Rtac+Rt

      Rltac=Rltac+Rlt

      Rttac=Rttac+Rtt

      Rattac=Rattac+Ratt



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



5999  continue



c  6999  format(/,'     Acumuladas       ',5f10.4/
c     ,   '------------------------')     

c      write(46,6999) rlac,rtac,rltac,rttac,rattac





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

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


      Waverage=Rlac*Sl+Rtac*St+Att*Rattac


      sinv=4.d0*ekf2*cos2if2

      FW=1.d0/dfloat(Jb2+1)*spft

c      factor=fcoup*fw*sinv *akpf/(2.d0*pi)**5*(4.d0*pi)**2*
c     *  Eb/amassb/2.d0 * amassp

      factor=fcoup*sinv*Epf*Eb/q*fw/(2.d0*pi)**3

      
      factor=cpmfac*factor

      fpwia=(4.d0*pi)**2*(Epf+amassp)/2.d0/Epf
  

      fconv=1.d16*(hc)**3


c
c     factor*Waverage es la seccion eficaz/S(R,Jr,c,j,Mr)
c     la funcion espectroscopica sera: Mb/Eb/Ef/Pf/SIGCC1
c

      fact2=factor*amassb/Eb/akpf/Epf*130.d0
       

c      spectf=fact2*Waverage/sigcc1



c
c     la parte axial
c


c      spa=fact2*att*rattac/sigcc1
c   
c     seccion eficaz en 10 -42 cm2/MeV

      esig=factor*fpwia*fconv*Waverage


      write(38,'(2f8.3,3e13.5,2e11.3)')ekf*hc,pm*hc,
     ,   esig, cosnu, q*hc,singam,cosgam


c      write(37,'(1f7.1,5e13.5)')pm*hc,spectf,spectro,spect,
c     , speclt

c      write(38,'(1f7.1,5e13.5)')pm*hc,ekf*hc,cosnu,esig,waverage,factor
      



      sigacos=sigacos+facos*esig

      sigacost=sigacost+esig

      if(ichcos.eq.0)then
      sigcosi=esig
      ichcos=1
      endif
	 

      facos=6.d0-facos

1300  continue

      sigcos=ajacob*(sigacos-sigcosi-esig)*stcos/3.d0
      sigcost=ajacob*(sigacost-esig/2.d0-sigcosi/2.d0)*stcos

      sigaekf=sigaekf+sigcos*facekf
      sigaekfts=sigaekfts+sigcost*facekf

      sigaekftt=sigaekftt+sigcost
      sigaekfst=sigaekfst+sigcos


      facekf=6.d0-facekf

      if(ichekf.eq.0)then
      sigekfti=sigcost
      sigekfi=sigcos
      ichekf=1
      endif 

      write(39,'(2f10.5,3e13.5,i5)')ekf*hc,pm*hc,sigcos/ajacob,
     ,sigcost/ajacob,ajacob,incos
      write(*,'(2e13.5,i5)')sigcos,sigcost,incos

c      write(38,*)stcos,esig,sigcos,factor,waverage

      
    1 continue

      print*,sigaekf,sigcos,sigekfi,sekf


      sigtol=(sigaekf-sigcos-sigekfi)*sekf/3.d0
      sigtolts=(sigaekfts-sigcost-sigekfti)*sekf/3.d0
      
      sigtoltt=(sigaekftt-sigcost/2.d0-sigekfti/2.d0)*sekf
      sigtolst=(sigaekfst-sigcos/2.d0-sigekfi/2.d0)*sekf


      write(37,'(f10.5,4e13.5,i5)')tpf*hc,sigtol,sigtolts,sigtoltt,
     , sigtolst,isimps
      write(*,'(4e13.5,i5)')sigtol,sigtolts,sigtoltt,sigtolst,isimps   

12345   continue

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



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

 

      akp=1.79285d0

      akn=-1.91304d0





      alph=1.d0/137.04

      pi=3.1415926    





      

      sin2w=0.2330

      Gf=1.1664e-11



      sqr3=sqrt(3.0)

      amass=939.d0

      hc=197.3d0



      am2v=840.d0**2

      am2a=1030.d0**2





c      f2s=-0.21

c      g1s=-0.19

c      f2s=0.

c      g1s=0.





      f1s=0.

      f1(0)=1.+f1s

      f2(0)=f2s+(akp+akn)

      g1(0)=g1s+0.55/3.



      f1(3)=0.5

      f2(3)=0.5*(akp-akn)

      g1(3)=0.63



      f1(8)=0.5*sqr3

      f2(8)=0.5*sqr3*(akp+akn)

      g1(8)=0.275/sqr3



      z0v=-0.5

      z0a=0.5

      z3v=1.-2*sin2w

      z3a=-1.

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

      z8a=-1./sqr3

c      print*,z0v,z3v,z8v,f2(0),f2(3),f2(8)

c      write(10,*)'proton-nu   proton anu  neutron nu  neutron-anu'



      Q2=-qnu2*hc*hc





      dip=1.+q2/am2v

      dip=1./dip/dip



      dips=1.+Q2/am2a

      dips=1./dips/dips



      tau=Q2/4./amass/amass







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

c

c     factores de forma EM del neutron y proton

c





      if (icharg.eq.1)then

        anom=akp

        ge0=1.0

      else

        anom=akn

        ge0=0.0

      endif









      Gm0=ge0+anom









      if(idip.eq.0)then

        dip=1.d0

        dips=1.d0

      endif



      gm=gm0*dip

      ge=ge0*dip



      if(icharg.eq.2)then

      ge=-anom*tau/(1.+5.6*tau)*dip

      endif







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

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







      gmt=eff1(icharg)+anom*eff2(icharg)

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

c      print*,'qnu2=',qnu2,ff1(1),ff2(1),gm/ge,gmt/get

c      if(ipl.eq.2)ff1=-anom*ff2

c      if(ipl.eq.3)ff2=0.0d0

c      if(ipl.eq.4)ff1=0.d0

c      falp=ff1

c      fbet=anom*ff2/2.d0/amassp



c      fff(1)=-anom*ff2/2.d0/amassp

c      fff(2)=ff1+anom*ff2

c      print*,fff(1),fff(2),'momentos'















      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





      fac=gf*gf/2.d0

c   /3.1415926*(197.3**2)*2.*amass*1e16



c     cross section in fm^2/MeV x 1e-16=cm2/MeV x 1.e-42)





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

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



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

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



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

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



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

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



      iii=icharg

c      write(13,'(i4,5e13.5)')iii,q2,ff1(iii),ff2(iii),gg1(iii)

c     i=1, proton

c     i=2, neutron



      if(ireac.ne.0)then



      ef1=ff1(icharg)

      ef2=ff2(icharg)

      eg1=gg1(icharg)



      fcoup=fac 



      else



      eg1=0.d0

      ef1=eff1(icharg)

      ef2=eff2(icharg)*anom

      fcoup=2.*(2.*pi)**2*alph*alph/qnu2/qnu2



      endif      

     

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



      end











      SUBROUTINE harmsq(spci,dpci)

      IMPLICIT real*8 (A-H,O-Z)

      include 'param.inc'

      common/harm/ysh(0:npw+1,-npw:npw)

*--------------------------------------------------*

*     CALCULATES spherical harmonics               *

*--------------------------------------------------*

c      print900,spci,dpci

c      write(13,900)spci,dpci

 900  format('sin in.=',f10.5,' cos in.=',f10.5,' sin fin.=',f10.5,

     ,' cos fin.=',f10.5)

      do 137 l=1,npw

      do 137 m=-npw,npw

      ysh(l,m)=0.d0

 137  continue

      PI=dacos(-1.d0)

      FOURPI=4.d0*PI

      spc=spci

      dpc=dpci

      if(dabs(spc).lt.1.d-7)then

        do 111 l=0,npw

        if=1

        if(dpc.lt.0.d0)if=(-1)**l

        ysh(l,0)=if*dsqrt((02.d0*l+1.d0)/fourpi)

 111  continue

      else



        DO 100 M=-npw,npw

*-------------------------------------------------------------------*

*     RECURENCE FORMULA FOR THE SPHER. HARM. ,KEEPTS IN YSH(L,M)    *

*     SUM IS THE TOTAL AMPLITUDE                                    *

*-------------------------------------------------------------------*

        A1=1.0D00

        A3=3.0D00

        IM=IABS(M)

        DM=Dfloat(M)

        DIM=Dfloat(IM)

C     DO 101 L=0,IM-1

C     YSH(L,M)=0.0D00

C101  CONTINUE

      IF(IM.EQ.0)THEN

      AUX1=1.0D0

      ELSE

      AUX1=SPC**IM

      if=1

      if(aux1.lt.0.d0)if=-1

      aux1=aux1*if**(im)

      ENDIF

      is=(-1)**im

      if(m.lt.0)is=1

c      print*,m,im,is,aux1,spc,2.d0**im,ea(im)

      YSH(IM,M)=is*DSQRT((2*DIM+A1)*dexp(eA(2*IM)))/DSQRT(FOURPI)/

     1((2.d0**IM)*dexp(eA(IM)))*AUX1

      YSH(IM+1,M)=DSQRT((2*DIM+A1)*(2*DIM+A3)/(DIM+DM+A1)/

     1(DIM-DM+A1))*DPC*YSH(IM,M)

c      print*,m,im,is,aux1,spc,2.d0**im,ea(im)

      DO 102 L=IM+1,npw

      DL=dfloat(L)

      CT=DSQRT((2*DL+A1)*(2*DL+A3)/(DL+A1+DM)/(DL+A1-DM))

      CI=DSQRT((DL+DM)*(DL-DM)/(2*DL+A1)/(2*DL-A1))

      YSH(L+1,M)=CT*(DPC*YSH(L,M)-CI*YSH(L-1,M))

 102  CONTINUE

 100  CONTINUE

      ENDIF

*-------------------------------------------------------------*

*     THIS PRINTS A TEST, Y10=DSQRT(3/4PI)*DCOS(THET)...       *

*-------------------------------------------------------------*

c      WRITE(*,99)YSHi(2,0)/DSQRT(5.D0/FOURPI),

c     ,YSHi(10,0)/DSQRT(21.D0/FOURPI),

c     ,YSHi(2,0),DSQRT(1.25D0/FOURPI)*(3.*DPC*DPC-1.0),

c     ,YSHi(2,-1),YSHi(2,1)

 99   FORMAT(2E12.5)

      return

      end

C

c

c

*----------------------------------------------------------------------*

*----------------------------------------------------------------------*

*  CLEBS-GORDAN COEF. <L1,L2,M1,M2,L,M>=CG(2*L1,2*L2,2*M1,2*M2,2*J,2*M)*

*----------------------------------------------------------------------*

      real*8 FUNCTION CGC(J1,M1,j2,M2,J,M)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      CNO=DSQRT(DFLOAT(J+1))

c      print*,j1,m1,j2,m2,j,m

      cno=cno*(-1)**(iabs(J1-J2+M)/2)

      CNO2=Sex3J(J1,J2,J,M1,M2,-M)

      CGC=CNO*CNO2

      RETURN

      END 

      real*8 FUNCTION CGC1(J1,M1,j2,M2,J,M)

      IMPLICIT real*8 (A-H,O-Z)

c      print*,j1,m1,j2,m2,j,m

      if((m1+m2-m).ne.0)then

      cgc1=0.d0

      return

      endif

      if((j.lt.iabs(j1-j2)).or.(j.gt.(j1+j2)))then

      cgc1=0.d0

      return

      endif

      if((iabs(m1).gt.j1).or.(iabs(m).gt.j))then

      cgc1=0.d0

      return

      endif 

      if((j-j1).eq.1)then

      if(m2.eq.-1)cgc1=dsqrt(dfloat(j1-m)/dfloat(j1+1)/2.d0)

      if(m2.eq.1)cgc1=dsqrt(dfloat(j1+m)/dfloat(j1+1)/2.d0)

c      print*,'caso 1',m2,cgc1

      return

      endif

      if((j-j1).eq.-1)then

      if(m2.eq.-1)cgc1=dsqrt(dfloat(j1+m)/dfloat(j1+1)/2.d0)

      if(m2.eq.1)cgc1=-dsqrt(dfloat(j1-m)/dfloat(j1+1)/2.d0)

c      print*,'caso2',m2,cgc1

      return

      endif

      print*,'algo falla'

      RETURN

      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





      real*8 function djl(l,xpm)

      implicit real*8 (a-h,o-z)





      if(l.eq.0) then

      djl=dsin(xpm)/xpm

      return

      endif



      if(l.eq.1) then

      djl=dsin(xpm)/xpm/xpm-dcos(xpm)/xpm

      return

      endif



      if(l.eq.2) then

      djl=(3.d0/xpm/xpm/xpm-1.d0/xpm)*dsin(xpm)-3.d0/xpm/xpm*dcos(xpm)

      return

      endif





      if(l.eq.3) then

      djl=dcos(xpm)-15.d0*dcos(xpm)/xpm/xpm+15.d0*dsin(xpm)/xpm**3-

     -    6.d0*dsin(xpm)/xpm

      djl=djl/xpm

      return

      endif



      if(l.eq.4) then

      djl=-105d0*dcos(xpm)/xpm**3+10.d0*dcos(xpm)/xpm+dsin(xpm)+

     +     105.d0*dsin(xpm)/xpm**4-45.d0*dsin(xpm)/xpm/xpm

      djl=djl/xpm

      return

      endif





      if(l.ge.5)then

      print*,' preparar version con hl(pm), esta ya no es eficiente'

      stop

      return

      endif



      end 







      SUBROUTINE HL(X)

      IMPLICIT  double precision(A-H,O-Z)

      include 'param.inc'

      COMMON/UNO/Hp(0:100)

      dimension ahl(0:120)

      if(x.lt.0.025d0)then

      hp(0)=1.d0

      do 10 l=1,npw+3

c      print*,dbeac(l+l+1)

      hp(l)=x**l/(dbeac(l+l+1))

 10   continue

      return

      endif

      IMAX=50+idint(x)



c      imax=60

      aHL(IMAX)=1.D-280 

      aHL(IMAX-1)=0.0D0

      DO 100 I=IMAX-1,1,-1

      aHL(I-1)=(2.D0*Dfloat(I)+1.D0)/X*aHL(I)-aHL(I+1) 

 100  CONTINUE

      REF=DSIN(X)/X 

      ESCALA=REF/aHL(0)

      DO 200 I=0,npw+3

      Hp(I)=ESCALA*aHL(I)

 200  CONTINUE

      RETURN

      END 

      SUBROUTINE azh(L,X,ZH)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      COMPLEX*16 ZI,ZH,ZEROZ,coef

      ZI=dCMPLX(0.D0,1.D0)

      ZEROZ=dCMPLX(0.D0,0.D0)

      ZH=ZEROZ

      DO 100 N=1,l+1

      COEF=FA(L+N-1)/FA(N-1)/FA(1+L-N)*(2.D0**(1-N))

      ZH=ZH+COEF/X**N*(ZI**(2+N-L))

 100  CONTINUE

      ZH=ZH*CdEXP(ZI*X)

      RETURN

      END 

C

      SUBROUTINE SPLINE(X,Y,N,YP1,YPN,Y2)

c      PARAMETER (NMAX=npunt1)

      implicit  real*8 (a-h,o-z)

      include 'param.inc'

      DIMENSION X(npunt2),Y(npunt2),Y2(npunt2),U(npunt2)

      IF (YP1.GT..99E30) THEN

        Y2(1)=0.d0

        U(1)=0.d0

      ELSE

        Y2(1)=-0.5d0

        U(1)=(3.d0/(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)

      ENDIF

      DO 11 I=2,N-1

        SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1))

        P=SIG*Y2(I-1)+2.d0

        Y2(I)=(SIG-1.d0)/P

        U(I)=(6.d0*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))

     *      /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P

11    CONTINUE

      IF (YPN.GT..99d30) THEN

        QN=0.d0

        UN=0.d0

      ELSE

        QN=0.5d0

        UN=(3.d0/(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))

      ENDIF

      Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.d0)

      DO 12 K=N-1,1,-1

        Y2(K)=Y2(K)*Y2(K+1)+U(K)

12    CONTINUE

      RETURN

      END

 

      SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y)

      implicit real*8 (a-h,o-z)

      include 'param.inc'

      DIMENSION XA(npunt2),YA(npunt2),Y2A(npunt2)

      KLO=1

      KHI=N

1     IF (KHI-KLO.GT.1) THEN

        K=(KHI+KLO)/2

        IF(XA(K).GT.X)THEN

          KHI=K

        ELSE

          KLO=K

        ENDIF

      GOTO 1

      ENDIF

      H=XA(KHI)-XA(KLO)

      IF (H.EQ.0.d0) PAUSE 'Bad XA input.'

      A=(XA(KHI)-X)/H

      B=(X-XA(KLO))/H

      Y=A*YA(KLO)+B*YA(KHI)+

     *      ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.d0

      RETURN

      END

      subroutine zeul(zs,zt,j,zw,jm)     

      implicit real*8 (a-h,o-z)

      include 'param.inc'

      dimension rzw(0:npw),dizw(0:npw),zw(0:j)

      complex*16 zs,zt,zw

      do 100 k=0,j

c        print*,k,j

        rzw(k)=dble(zw(k))

        dizw(k)=dimag(zw(k))

 100  continue

      rterm=dble(zt)

      diter=dimag(zt)

      call eulsum(rsum,rterm,j,rzw)

      call eulsum(disum,diterm,j,dizw)

      zs=dcmplx(rsum,disum)                  

      do 101 k=0,j

        print*,j,k,npw,'    estamos aqui'

        zw(k)=dcmplx(rzw(k),dizw(k))

        print*,j,k,npw,zw(k),'    estamos aqui'

 101  continue

      return

      end

c

c

c

      SUBROUTINE EULSUM(SUM,TERM,JTERM,WKSP)

      implicit real*8 (a-h,o-z)

      include 'param.inc'

      DIMENSION WKSP(0:jterm+1)

c      print*,'jterm=',jterm

      IF(JTERM.EQ.0)THEN

        NTERM=1

        WKSP(0)=TERM

        SUM=0.5d0*TERM

      ELSE

        TMP=WKSP(0)

        WKSP(0)=TERM

        DO 11 J=0,NTERM

c

c          print*,jterm,j+1

c

          DUM=WKSP(J+1)

          WKSP(J+1)=0.5d0*(WKSP(J)+TMP)

          TMP=DUM

11      CONTINUE

        IF(dABS(WKSP(NTERM+1)).LE.dABS(WKSP(NTERM)))THEN

c        print*,nterm+1

          SUM=SUM+0.5d0*WKSP(NTERM+1)

          NTERM=NTERM+1

        ELSE

          SUM=SUM+WKSP(NTERM+1)

        ENDIF

      ENDIF

      RETURN

      END

C 

C     CLEBS-GORDAN  <L1,L,M1,M2,L,M>=CG(2*L1,2*L2,2*M1,2*M2,2*J,2*M)   

C                                                                        

      DOUBLE PRECISION FUNCTION CG(J1,m1,j2,m2,J,M)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      CNO=DSQRT(DBLE(J+1))*DBLE((-1)**((J1-J2+M)/2))

      CG=CNO*Sex3J(J1,J2,J,M1,M2,-M)

      RETURN

      END 

C

C     3 J   S I M B O L S, L A , L B , L C = 2 * ( J 1 , J 2 , J 3 )

C

C                    N L F , N B T , N G A = 2 * ( M 1 , M 2 , M 3 )

C

      DOUBLE PRECISION FUNCTION SIM3J(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

      SIM3J=0.0D00

      ELSE

      A1=FA((LA+NLF)/2)

      A1=DSQRT(A1)

      A2=FA((LB+NBT)/2)

      A2=DSQRT(A2)

      A3=FA((LB-NBT)/2)

      A3=DSQRT(A3)

      A4=FA((LC+NGA)/2)

      A4=DSQRT(A4)

      A5=FA((LC-NGA)/2)

      A5=DSQRT(A5)

      A6=FA((LA-NLF)/2)

      A6=DSQRT(A6)

      D1=DELTA(LA,LB,LC)

      isg=1

      if(mod((la-lb-nga)/2,2).ne.0)isg=-1

      CNO=isg*DSQRT(D1)

      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=FA(N)/A1

      B2=FA((LC-LB+NLF)/2+N)/A2

      B3=FA((LC-LA-NBT)/2+N)/A3

      B4=FA((LA+LB-LC)/2-N)/A4

      B5=FA((LA-NLF)/2-N)/A5

      B6=FA((LB+NBT)/2-N)/A6

      CNTER=DBLE((-1)**N)*B1*B2*B3*B4*B5*B6 

      IF(CNTER.NE.0.0D00)SUM=SUM+1.D0/CNTER

100   CONTINUE

      SIM3J=CNO*SUM 

      ENDIF

      RETURN

      END 

C

C     F A C T O R I A L ( N )  

C



      SUBROUTINE CONSFA

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      COMMON/FACT/AFA(350)

      FB=1.0D00

      DO 10 I=1,250 

      FB=FB*DBLE(I) 

      AFA(I)=FB

 10   CONTINUE

      RETURN

      END 

C

      DOUBLE PRECISION FUNCTION FA(N)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      COMMON/FACT/AFA(350)

      IF(N.EQ.0)THEN

      FA=1.0D0

      ELSE

      FA=AFA(N)

      ENDIF

      RETURN

      END 

      DOUBLE PRECISION FUNCTION SIM6J(J1,J2,J3,J4,J5,J6)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      SUM=0.D0

      DO 100 N=MAX0(J1+J2+J3,J1+J5+J6,J4+J2+J6,J3+J4+J5),

     ,MIN0(J1+J2+J4+J5,J2+J3+J5+J6,J3+J1+J6+J4),2

      C1=dble((-1)**(N/2))*FA((N+2)/2)/FA((N-J1-J2-J3)/2)

      C2=FA((N-J1-J5-J6)/2)*FA((N-J4-J2-J6)/2)

      C3=FA((N-J4-J5-J3)/2)*FA((J1+J2+J4+J5-N)/2)

      C4=FA((J2+J3+J5+J6-N)/2)*FA((J3+J1+J6+J4-N)/2)

      SUM=SUM+C1/C2/C3/C4

  100 CONTINUE

      IF(SUM.NE.0.D0)THEN

      SIM6J=SUM*DSQRT(delta(j1,j2,j3)*delta(j1,j5,j6)*delta(j4,j2,j6)*

     *delta(j4,j5,j3))

      else

      sim6j=0.0d0

      endif

      return

      end

C 

C     R a c a h   W - c o e f f i c i e n t

C

C         a s   a   f u n c t i o n   

C

C        o f   6 - j   S y m b o l s

C

      function w(j1,j2,l2,l1,j3,l3)

      implicit double precision(a-h,o-z)

      fase=dble((-1)**((j1+j2+l1+l2)/2))

      w=fase*sim6j(j1,j2,j3,l1,l2,l3)

      return

      end

      DOUBLE PRECISION FUNCTION DELTA(LA,LB,LC)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      IF((LA+LB+LC).LT.0)THEN 

      DELTA=0.0D00

      RETURN

      ELSE

      DELTA=FA((LA+LB-LC)/2)*FA((LB+LC-LA)/2)*FA((LC+LA-LB)/2)/

     1FA((LA+LB+LC)/2+1)

      RETURN

      ENDIF

      END

