      Program chap1
c     Reads pi-nucleon phase shifts from file delal
c     S. Almehed, C Lovelace, Nucl Phys B40 (1972) 157
c     Produces disk file: chap1.dat with two columns
c     Center of mass energy and eta or d for selected nwave
c     1<=nwave <= 14 , eta :real part and d: imaginary part
c     This example is for nwave=1, and eta
c        nwave code
c     1-----s31     2-----p31     3-----p33     4-----s11
c     5-----p11     6-----p13     7-----d33     8-----d35
c     9-----f35     10----f37     11----d13     12----d15
c     13----f15     14----f17
c
c     gets rekf(14,100) aimtf(14,100), see routine argand
      Implicit Real*8 (a-h,o-z)
      Dimension ecm(100),eta(14,100), d(14,100),retf(14,100)
      Dimension aimtf(14,100),ceta(14,100),cd(14,100),cecm(100)
      pi=3.141592
      nwaves=14
      nes=71  !select nes=71 for No Rowe, Salomon, Landau phaseshifts
      nei=1
      Open(1,file='delal',status='unknown')
      Read (1,170) (ecm(ne),(eta(nwave,ne),d(nwave,ne),nwave=1,nwaves),
     1ne=nei,nes)
c     d tabulated in degrees, now convert to radians  
      Do  nwave=1,14
        Do  ne=1,nes
            d(nwave,ne)=d(nwave,ne)*pi/180.0
        End Do
      End Do
c     select i=0 no Rowe Salomon Landau phase shifts
c            i=1  yes R S L phase shifts
      i=1
      if(i.eq.0)Then
         Call argand(eta,d,retf,aimtf,nes)
         Call crossect(ecm,retf,aimtf,nes)  
      End If
      if(i.eq.1)Then 
        come=ecm(1)
        Call rosalan(come,ceta,cd,cecm)
c       output ceta is again in degrees convert to radians
        Do nwave=1,14
          Do ne=1,25
            cd(nwave,ne)=cd(nwave,ne)*pi/180.0 ! ro radians
          End Do
        End Do
        nes=68   ! number of energies considered 25+70=95 max
c       first 25  ecm eta d from RSL shifts, 26-90 as before
        Do nwave=1,14
          Do ne=26,nes
            cd(nwave,ne)=d(nwave,ne-24)
            ceta(nwave,ne)=eta(nwave,ne-24)
            cecm(ne)=ecm(ne-24)
          End Do
        End Do
        Call argand(ceta,cd,retf,aimtf,nes)
        Call crossect(cecm,retf,aimtf,nes)
 170    Format (13f6.1/8f6.1/8f6.1)
      End If
      Close(1)
      Stop
      End
      
      Subroutine argand(eta,d,retf,aimtf,nes)
c     ecm, eta, d, nes  are inputs to the routine
c     Returns args.  retf and aimtf and also
c     the output is a file "chap1ar.dat" with
c     the 2 columns real and imag. parts of partial
c     amplitudes f's. With this file you can 
c     plot Argand diagrams. The f's for each wave are
c     f=-0.5i(eta*exp(2id -1))
c     separating the real and imaginary parts
c     Compare with S.Almehed and C Lovelace p 163
c     i=0 :for Argand diagram
c     select below the nwave you like (from 1 to 14)    
      Implicit Real*8 (a-h,o-z)
      Real*8 mpi
      Dimension eta(14,100),d(14,100),retf(14,100)
      Dimension aimtf(14,100)
   
      pi=3.141592
      mpi=139.578
      Open(6,file='chap1ar.dat', status='unknown')
      Do  nwave=1,14
        Do  ne=1,nes
           twodel=2.0*d(nwave,ne)
           retf(nwave,ne)=0.5*eta(nwave,ne)*sin(twodel)
           aimtf(nwave,ne)=0.5*(1.0-eta(nwave,ne)*cos(twodel))
c          retf(nwave,ne)=0.5*sin(twodel) ! eta=1
c          aimtf(nwave,ne)=0.5*(1.0-cos(twodel)) !eta=1
        End Do
      End Do
      nwave=3
      Do  ne=1,nes-1
        Write(6,100)retf(nwave,ne),aimtf(nwave,ne)!argand dia
      End Do
 100  Format(1h ,e15.5,e15.4)
      Close(6)
      Return
      End

      Subroutine kinemat(come,kappa,pp)
c     input: center of mass energy
c     output: kappa, pp :com pion momentum and lab pion momentum
c     Computes kappa: center of mass momentum from total ecm energy
c     computes lab pion energy from total com energy ecm
c     epi=(ecm**2-mpi**2-mn**2)/(2mn)
c     then finds lab pion momentum pp=sqrt(epi**2-mpi**2)
      Implicit Real*8 (a-h,o-z)
      Real*8 mn,mpi,kappa
      mpi=139.578 !pion mass
      mn=938.272 !proton mass
      epi=(come**2-mpi**2-mn**2)/(2*mn)
      pp=sqrt(epi**2-mpi**2)
      s = come**2
      kappa = Sqrt((s-(mn+mpi)*(mn+mpi))*(s-(mn-mpi)*(mn-mpi))/4./s)
      Return
      End

      Subroutine crossect(ecm,retf,aimtf,nener)  
c     input: ecm  aimtf and    
c     nener: number of com energies and aimmtf's needed
c     the output is a file "chap1st.dat"
c     that contains either the total cross sect. versus pi-lab momemt.
c     or the contribution of any partial wave to the cross section
c     or tpilab and reft for a wave
c     or tplilab and aimtf for a wase
c     depending of your choice of the i value
c      write(6,100) statement
      Implicit Real*8 (a-h,o-z)
      Real*8 kappa, mpi
      Dimension ecm(100),  retf(14,100),aimtf(14,100)
      pi=3.141592
      hbarc=197.3289
      mpi=139.578 ! pion mass      
c     fact : to give total cross section in mb
      Open(6,file='chap1st.dat', status='unknown')
c     call kinemat gives com enery (come) and
c     gets kappa: com pion momentum
c          pp:    lab. pion momentum  
c     tpilab: kinetic energy of the pion in lab system
c     i=0 tpilab and total cross section
c     i=1 tpilab and cntrbtn of any partial wave to crossection
c     i=2 tpilab and retf for a wave
c     i=3 tpilab and imtf for a wave
      i=0
      Do  ne=1,nener
        come=ecm(ne)
        Call kinemat(come,kappa,pp)
        tpilab=sqrt(pp*pp+mpi*mpi)-mpi
        fact=10.0*4.0*pi*(hbarc/kappa)**2
        swve=aimtf(1,ne)
        pwve=2.0*aimtf(3,ne)+aimtf(2,ne)
        dwve=3.0*aimtf(8,ne)+2.0*aimtf(7,ne)
        fwve=4.0*aimtf(10,ne)+3.0*aimtf(9,ne)
        sigtot=fact*(swve+pwve+dwve+fwve)
c       first write for sigma total
c       second write for contribution of a wave to sigma total
c       uncomment one and comment the other according to your
c       choice
        if(i.eq.0)Write(6,100) tpilab,sigtot
        if(i.eq.1)Write(6,100)tpilab,fact*2.0*aimtf(3,ne)
        if(i.eq.2)Write(6,100)tpilab,retf(3,ne)
        if(i.eq.3)Write(6,100)tpilab,aimtf(3,ne)
      End Do
      Close(6)
 100  Format(1h ,e15.5,e15.4)
      Return
      End
      
      Subroutine rosalan(come,ceta,cd,cecm)
c     Rowe, Salomon, Landau partial amplitudes
c     Phys.Rev C 18, (1978), 584-589
c     input: come=ecm(1)     
c     output: ceta(14,100), cd(14,100) and cecm(100)
c     eta, delta and ecm with its first 25 data including
c     the RSL phase shifts and corresponding com energies
c     parts taken from routine tpicm from lpott
c     the output is a file "chap1rsl.dat" with
c     first column: pp: lab pion momentum
c     second column: cd: phase shift in degrees for a
c     freely selected wave (1 to 14), 3 in this case 
      Implicit Real*8 (a-h,o-z)
      Real*8 kappa,kaph,dumkap
      Dimension ceta(14,100),cd(14,100),cecm(100)
      Open(6,file='chap1prsl.dat', status='unknown')      
      pi=3.141592
      Call kinemat(come,kappa,pp)
      If (kappa.le.5.) kappa = 6.
      kaph = (kappa-5.)/20.
      kappa = 0.
c      Do  nkap=1,25  ! as in tpicm
c       comment last line and uncomment next to plot phashifts
       Do  nkap=1,51  !for com pi momentum, up to 51 for P33graph
          kappa = kappa+kaph
          If (nkap.le.6) kappa = (nkap*1.5)
          Do  ni=1,14                   !for partial waves
             Call tpirsl (kappa,ni,ddum,w,tr,ti)
c            w is given by last call is com energy
             ceta(ni,nkap) = 1.
             cd(ni,nkap) = ddum*180./pi
c            tangent jumps at 90 prom + to -, correct to sum phases
c            for ni=3 
             if((ni.eq.3).and.(cd(3,nkap).lt.0.0)) cd(3,nkap)=
     1        180.0+cd(3,nkap)
         End Do
         cecm(nkap) = w
c        input w in kinemat and use pp: lab pion kin. ener.
         Call kinemat(w,dumkap,pp)         
c        To plot  S31 P31,P33, P11 phaseshits in degrees        
         Write(6,*)pp,cd(3,nkap) !here 
      End Do
      Close(6)
      Return
      End

      Subroutine tpirsl (kappa,ni,ddum,w,tr,ti)
c *** see RH Landau, Program lpott, 1981
c     calc of pin phases  and amplitudes with rowe,salomon,Landau
c     formula, or Landau s for higher l
      Implicit Real*8 (a-h,i,k,m,o-z)
      Dimension a(14), b(14), c(14), xx(14), w0(14), q0(14), gam(14)
c       use rowe,salomon,Landau for low e s,p,d13,d15;Landau for others
c
      data mpi,mn,pi,kold/139.578e+0,938.903,3.141593e+0,0./,
     1 a/-11.2,-2.91,11.4,16.8,-5.7,-1.31,11.2e-13,-3.15e-12,
     2 8.31e-18,5.46e-17,0.109,0.112,-38.26e-18,7.43e-18/,
     3 b/-30.7,3.45,-15.4,-35.4,25.4,1.22,-0.369e-14,0.118e-13,
     4 -0.363e-19,-0.236e-18,-0.031,-0.270,0.205e-18,-0.222e-19/,
     5 c/21.,-1.5,7.2,27.,-29.,-0.4,0.,0.,0.,0.,0.003,0.19,0.,0./,
     6 xx/.31,.22,.99,.44,.61,.23,0.,0.,0.,0.,.54,.43,0.,0./,
     7 w0/1655.,1850.,1233.,1550.,1435.,1815.,0.,0.,0.,0.,
     8 1525.,1670.,0.,0./,
     9 q0/550.,678.,228.,477.,393.,656.,0.,0.,0.,0.,459.,560.,0.,0./,
     a gam/170.,200.,116.,105.,230.,255.,0.,0.,0.,0.,125.,155.,0.,0./
      If (abs(kappa-kold).lt.1.e-03) GoTo 10
      kold = kappa
      rhl = kappa*kappa
      If (kappa.eq.0.) kappa = 1.e-10
      epikap = Sqrt(mpi*mpi+rhl)
      enkap = Sqrt(mn*mn+rhl)
      w = epikap+enkap
      s = w*w
      factor = (-1./kappa)*w/4./pi/pi/epikap/enkap
      q = kappa/mpi
      q2 = q*q
      q4 = q2*q2
c     determine l value for formula
 10   If ((ni.eq.1).or.(ni.eq.4)) l = 0
      If ((ni.eq.2).or.(ni.eq.3).or.(ni.eq.5).or.(ni.eq.6)) l = 1
      If ((ni.eq.7).or.(ni.eq.8).or.(ni.eq.11).or.(ni.eq.12)) l = 2
      If ((ni.eq.9).or.(ni.eq.10).or.(ni.eq.13).or.(ni.eq.14)) l = 3
      jj = 2*l+1
      If ((ni.eq.7).or.(ni.eq.8).or.(ni.eq.9).or.(ni.eq.10).or.(ni.eq.13
     1).or.(ni.eq.14)) GoTo 20
      rhl = q0(ni)/mpi
      del = (q**jj)*(a(ni)*1.e-2+b(ni)*q2*1.e-3+c(ni)*q4*1.e-4+xx(ni)*
     1gam(ni)*w0(ni)/(rhl**jj)/(w0(ni)*w0(ni)-s))
      ddum = atan(del)
      GoTo 30
c     no salomon fit, use old rhl fit
 20   ddum = (a(ni)+b(ni)*kappa)*(kappa**jj)*pi/180.
 30   tr = 0.5*factor*sin(2.*ddum)
      ti = 0.5*factor*(1.-cos(2.*ddum))
      Return
      End











