! fsi_test.f

! program to test whether the pi-N phase shifts are being calculated correctly.

      implicit none

      real p_cm,s,t_beam,p_beam,cos_cm
      real xmpi,xmp,epi,ep,phase,formula
      real d2sig,z0,z1,z2

      integer ichannel,l,i

      xmpi = 0.13958
      xmp  = 0.93828

      do i=1,6
         if (i.eq.1) then
            write(6,*)' T=3/2, S31 Phase Shifts '
            ichannel=2
            l=0
         elseif (i.eq.2) then
            write(6,*)' T=3/2, P33 Phase Shifts '
            ichannel=6
            l=1                 ! P33
         elseif (i.eq.3) then
            write(6,*)' T=3/2, P31 Phase Shifts '
            ichannel=5
            l=1                 ! P31
         elseif (i.eq.4) then
            write(6,*)' T=1/2, S11 Phase Shifts '
            ichannel=1
            l=0                 ! S11
         elseif (i.eq.5) then
            write(6,*)' T=3/2, P13 Phase Shifts '
            ichannel=4
            l=1                 ! P13
         elseif (i.eq.6) then
            write(6,*)' T=3/2, P11 Phase Shifts '
            ichannel=3
            l=1                 ! P11
         endif
         
         do p_cm=0.,0.300,0.050
            epi = sqrt(p_cm**2 + xmpi**2)
            ep  = sqrt(p_cm**2 + xmp**2)
            s= (epi+ep)**2
            phase = asind(formula(s,p_cm,l,ichannel))
            write(6,*)p_cm,phase
         enddo
      enddo

! now calculate differential cross sections

      do i=1,2
         if (i.eq.1) then       ! pi+p (pure T=3/2)
            ichannel=1   
            write(6,*)' pi+p '
         elseif (i.eq.2) then   ! pi-p (T=1/2 and 3/2)
            ichannel=2
            write(6,*)' pi-p '
         endif
         
         do t_beam=0.0001,1.00,0.001 ! Tpi
!           write(6,*)' Starting Tpi=',t_beam
            
            s = (xmpi+xmp)**2 + 2.*t_beam*xmp
            p_beam = sqrt( (t_beam+xmpi)**2 - xmpi**2 )
            p_cm = p_beam*xmp/sqrt(s)
            call phaseshifts(ichannel,p_cm,s,z0,z1,z2)
            
            do cos_cm=-1.0,1.0,0.5
               d2sig=z0 + z1*cos_cm + z2*(cos_cm**2)
               write(6,*)cos_cm,d2sig
            enddo
         enddo
      enddo

      end