**************************************************
      subroutine phaseshifts(ich,q,s,z0,z1,z2)
**************************************************
! this routine returns the CM angular distribution parameters for pi-N
! elastic scattering.

! based on work by A. Shinozaki
! gh 02.06.20

! ich = Isospin channel (1: pure T=3/2,  2: mixed T=1/2,3/2)
! q  = pi-N CM momentum (GeV/c)
! s  = pi-N CM energy^2 (GeV)
! z0 = cos^0 parameter
! z1 = cos^1 parameter (in pi-N CM frame)
! z2 = cos^2 parameter

      implicit none
**** inputs
      integer ich
      real s,q
**** outputs
      real z0,z1,z2
*** local variables
      real formula,reck2,hbarc2
      real ss31,cs31,sp33,cp33,sp31,cp31,ss11,cs11,sp13,cp13,sp11,cp11
      complex fs31,fp33,fp31,fs11,fp13,fp11

      hbarc2=0.389379                  ! mb(GeV/c)^2
      reck2=hbarc2/q**2                ! mb
  
      ss31=formula(s,q,0,2)            ! S31
      cs31=sqrt(1.0-ss31**2)
      fs31=cmplx(ss31*cs31,ss31**2)

      sp33=formula(s,q,1,6)            ! P33
      cp33=sqrt(1.0-sp33**2)
      fp33=cmplx(sp33*cp33,sp33**2)

      sp31=formula(s,q,1,5)            ! P31
      cp31=sqrt(1.0-sp31**2)
      fp31=cmplx(sp31*cp31,sp31**2)

      if (ich.eq.1) then               ! pure isospin 3/2 channel (e.g. pi+p)
         fs11=cmplx(0.,0.)
         fp11=cmplx(0.,0.)
         fp13=cmplx(0.,0.)
         
! Verify against the I=3/2 equations of Lock & Measday (#5-32, p.120).
!         z0=reck2*( ss31**2 + (sp33*cp31-cp33*sp31)**2 )
!         z1=reck2*2.0*ss31*( 2.0*sp33*(cp33*cs31+sp33*ss31) +
!     *                         sp31*(cp31*cs31+sp31*ss31) )
!         z2=reck2*(3.0*(sp33**2) + 6.0*sp31*sp33*(cp33*cp31+sp33*sp31))
!      
!         write(6,*)' measday: ',z0,z1,z2

      else                             ! mixed isospin 1/2,3/2 channel (e.g. pi-p)

         ss11=formula(s,q,0,1)         ! S11
         cs11=sqrt(1.0-ss11**2)
         fs11=cmplx(ss11*cs11,ss11**2)

         sp13=formula(s,q,1,4)         ! P13
         cp13=sqrt(1.0-sp13**2)
         fp13=cmplx(sp13*cp13,sp13**2)

         sp11=formula(s,q,1,3)         ! P11
         cp11=sqrt(1.0-cp11**2)
         fp11=cmplx(sp11*cp11,sp11**2)    

      endif

! equations for A,B,C from Appendix F of Aki's thesis.

      z0=reck2*real(                                            ! cos^0
     *       (fs31+2.0*fs11)*conjg(fs31+2.0*fs11)
     *     +        ( (fp33+2.0*fp13)-(fp31+2.0*fp11) )
     *        *conjg( (fp33+2.0*fp13)-(fp31+2.0*fp11) ) )

      z1=reck2*2.0*real(                                        ! cos^1
     *     (fs31+2.0*fs11)*conjg( 2.0*(fp33+2.0*fp13)+(fp31+2.0*fp11) ))

      z2=reck2*3.0*real(                                        ! cos^2
     *       (fp33+2.0*fp13)*conjg(fp33+2.0*fp13)
     *     + 2.0*(fp33+2.0*fp13)*conjg(fp31+2.0*fp11) )

      if (ich.eq.2) then
         z0=z0/9.
         z1=z1/9.
         z2=z2/9.
      endif

      return
      end


**************************************************
      real function formula(s,k,l,ij)
**************************************************
! written by A. Shinozaki, Oct/00.

! this function subroutine returns equation 2.25 (p. 812) of Feshbach,
! "Theoretical Nuclear Physics".  This is the parameterization of the
! pi-N phase shifts for pion energies less than 400 MeV (LAB).

! s = pi-N CM energy^2
! k = pi-N CM momentum
! l = pi-N CM angular momentum
! ij= parameter selecting which set of phaseshift parameters to use

      implicit none
      real s,k,kom,kom2,ss,s0,tangent
      real x(6),ss0(6),k0(6),G0(6),b(6),c(6),d(6)
      real xmpi
      integer ij,l,l21

**********************************************************************
* Phase shift data are from G. Rowe, M. Solomon and R.H. Landau
! PRC 18 (1978) 584.
! L_2t2j=        S11      S31      P11      P13      P31      P33
      data x  /  0.44,    0.31,    0.61,    0.23,    0.22,    0.99   /
      data ss0/  1.550,   1.655,   1.435,   1.815,   1.850,   1.233  /
      data k0 /  0.477,   0.550,   0.393,   0.656,   0.678,   0.228  /
      data G0 /  0.105,   0.170,   0.230,   0.255,   0.200,   0.116  /
      data b  /  0.168,  -0.112,  -0.0571, -0.0131, -0.0291,  0.114  /
      data c  / -0.0354, -0.0307,  0.0254,  0.00122, 0.00345,-0.0154 /
      data d  / 27.E-4,  21.E-4, -29.E-4,  -0.4E-4, -1.5E-4,  7.2E-4 /
**********************************************************************

      xmpi = 0.13957
      kom  = k/xmpi
      ss   = s
      if (k.gt.0.353) then        ! PSA limit @ Tpi=400 MeV
         kom = 0.353/xmpi         
         ss  = 1.383**2
      endif
      kom2 = kom**2
      l21  = 2*l+1

      s0=ss0(ij)**2
      if (abs(s0-ss) .lt. 1.e-10) then
         formula=1.0
      else
         tangent = (kom**l21)* ( b(ij) + (c(ij)+d(ij)*kom2)*kom2 )
     &        + x(ij)* ( (k/k0(ij))**l21 ) * G0(ij)*ss0(ij)/(s0-ss)
         formula=tangent/sqrt(1.+tangent**2)           ! sin(delta_Ltj)

! Note on quadrant corrections:
! The phaseshifts (delta_Ltj) are from 0-180 degrees.
! The tangent function is positive from 0-90 degrees and negative from 90-180 deg.
! Sine is positive from 0-180 deg while Cosine is negative from 90-180 deg.
! This means that "formula" will give the wrong sign (negative) for sin(delta)
! from 90-180 deg.  However, this works out okay, since the cosine calculated 
! in the calling routine is always positive, and the real part of the partial 
! wave amplitude (cos*sin) will have the correct sign for 0-180 deg.
! The imaginary part (sin**2) is always positive.  
! Thus, everything works out correctly.

      endif

      return
      end