      real(8) function fcn(xcos)
      implicit none
      real(8) :: xcos
      fcn = xcos
      end function fcn
!     -------------------------------------------------------------------------
!     differential cross section definition
!     -------------------------------------------------------------------------
!     spin-3/2 resonance model
!     -------------------------------------------------------------------------
!     Differential
!     -------------------------------------------------------------------------
      real(8) function dxs_spin32(egam,tman,xGamr,xMassr)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,tmannew,xGamr,xMassr
      real(8) :: a0,ap,Anorm,b0,st,xBrn
      real(8) :: a12flag,a32flag
!     Internal Variables
      real(8) :: dxsnorm,dxs_spin32temp,dfunc,deltafunc,jr,xcos
      real(8), dimension (24) :: listdfunc,listdelfunc
      integer :: i,angmom
!     External Variables
      real (8),external :: ampsqr
      real (8), external :: d3212p32p,d3212p12p,d3212p12m,d3212p32m
      real (8), external :: d3212m32p,d3212m12p,d3212m12m,d3212m32m
      real (8), external :: d3232p32p,d3232p12p,d3232p12m,d3232p32m
      real (8), external :: d3232m32p,d3232m12p,d3232m12m,d3232m32m
      dxsnorm =  (6084.375d0*alpha)/(egam**2*xMpro**2)
      angmom = 0
      jr = 3.d0/2.d0
      if (tman.gt. 0.d0) then
         xcos = 1.d0
!         write(0,*) xcos,tman
         tmannew =  (-2.d0*egam**2*xmpro + xmjpsi**2*xmpro + 
     -    egam*(xmjpsi**2 + xcos*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -        dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2 - 
     -            4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -	(xmpro*(2.d0*egam + xmpro)))))/
     -  (2.d0*egam + xmpro) 
      else
         xcos =   (2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      a12flag    = choice(1)
      a32flag    = choice(2)
      listdelfunc(1)  = 0.d0
      listdelfunc(2)  = 1.d0
      listdelfunc(3)  = 0.d0
      listdelfunc(4)  = 0.d0
      listdelfunc(5)  = 0.d0
      listdelfunc(6)  = 0.d0
      listdelfunc(7)  = 0.d0
      listdelfunc(8)  = 0.d0
      listdelfunc(9)  = 0.d0
      listdelfunc(10) = 1.d0
      listdelfunc(11) = 0.d0
      listdelfunc(12) = 0.d0
      listdelfunc(13) = 1.d0
      listdelfunc(14) = 0.d0
      listdelfunc(15) = 0.d0
      listdelfunc(16) = 0.d0
      listdelfunc(17) = 0.d0
      listdelfunc(18) = 0.d0
      listdelfunc(19) = 0.d0
      listdelfunc(20) = 0.d0
      listdelfunc(21) = 0.d0
      listdelfunc(22) = 0.d0
      listdelfunc(23) = 0.d0
      listdelfunc(24) = 1.d0
      listdfunc(1)  = d3212p32p(xcos)*a12flag
      listdfunc(2)  = d3212p12p(xcos)*a12flag
      listdfunc(3)  = d3212p12p(xcos)*a12flag
      listdfunc(4)  = d3212p12m(xcos)*a12flag
      listdfunc(5)  = d3212p12m(xcos)*a12flag
      listdfunc(6)  = d3212p32m(xcos)*a12flag
      listdfunc(7)  = d3212m32p(xcos)*a12flag
      listdfunc(8)  = d3212m12p(xcos)*a12flag
      listdfunc(9)  = d3212m12p(xcos)*a12flag
      listdfunc(10) = d3212m12m(xcos)*a12flag
      listdfunc(11) = d3212m12m(xcos)*a12flag
      listdfunc(12) = d3212m32m(xcos)*a12flag
      listdfunc(13) = d3232p32p(xcos)*a32flag
      listdfunc(14) = d3232p12p(xcos)*a32flag
      listdfunc(15) = d3232p12p(xcos)*a32flag
      listdfunc(16) = d3232p12m(xcos)*a32flag
      listdfunc(17) = d3232p12m(xcos)*a32flag
      listdfunc(18) = d3232p32m(xcos)*a32flag
      listdfunc(19) = d3232m32p(xcos)*a32flag
      listdfunc(20) = d3232m12p(xcos)*a32flag
      listdfunc(21) = d3232m12p(xcos)*a32flag
      listdfunc(22) = d3232m12m(xcos)*a32flag
      listdfunc(23) = d3232m12m(xcos)*a32flag
      listdfunc(24) = d3232m32m(xcos)*a32flag
!     Differential cross section calculation
      dxs_spin32temp = 0.d0
      if (tman .lt. 0.d0) then
         do i = 1,24
            dxs_spin32temp = dxs_spin32temp + dxsnorm*
     -	ampsqr(egam,tman,xGamr,xMassr,listdfunc(i)
     -	,listdelfunc(i),angmom,jr,xcos)
         enddo
      else
         do i = 1,24
            dxs_spin32temp = 
     -	dxs_spin32temp + dxsnorm*
     -	ampsqr(egam,tmannew,xGamr,xMassr,listdfunc(i)
     -	,listdelfunc(i),angmom,jr,xcos)
         enddo
      endif
      dxs_spin32 = dxs_spin32temp
      end function dxs_spin32
!     -------------------------------------------------------------------------
!     KLL
!     -------------------------------------------------------------------------
      real(8) function kll32(egam,tman,xGamr,xMassr)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,tmannew,xGamr,xMassr
      real(8) :: a0,ap,Anorm,b0,st,xBrn
!     Internal Variables
!     External Variables
      real (8), external :: ampsqr32pp,ampsqr32pm
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
!     Calculation
      kll32 = (ampsqr32pp(egam,tman,xGamr,xMassr)-
     -  ampsqr32pm(egam,tman,xGamr,xMassr))/
     -  (ampsqr32pp(egam,tman,xGamr,xMassr)+
     -  ampsqr32pm(egam,tman,xGamr,xMassr))
      end function kll32
!     -------------------------------------------------------------------------
!     Differential smear
!     -------------------------------------------------------------------------
      real(8) function dxs_spin32_smear(egam,tman,xGamr,xMassr)
      use math_constants
      use hadron_properties
      use storage
      use size
      implicit none
      real(8),dimension(20000) :: enerd
      real(8),dimension(20000) :: u1rd
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
!     Internal Variables
      real(8) :: gaussener,gausstman,smear,xcos,ethresh
      integer :: i1,n1p
      real(8) :: a0,ap,Anorm,b0,st,xBrn
!     External Variables
      real (8),external :: dxs_spin32
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      smear = choice(4)
      ethresh=(xmjpsi**2+2.d0*xmpro*xmjpsi)/(2.d0*xmpro)
      call dsg20r(ethresh,emax,nintegprec,enerd,n1p)
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      do i1=1,n1p
         gaussener=enerd(i1)
         gausstman=(-2.d0*gaussener**2*xmpro + xmjpsi**2*xmpro + 
     -    	gaussener*(xmjpsi**2 + xcos*dsqrt(xmpro*
     -		(2.d0*gaussener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*gaussener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(gaussener + xmpro))/
     -		(xmpro*(2.d0*gaussener + xmpro)))))/
     -  	(2.d0*gaussener + xmpro)
         u1rd(i1)=  dxs_spin32(gaussener,gausstman-.0000000001
     -	,xGamr,xMassr)*
     c		dexp(-(egam-gaussener)**2/(2.d0*smear**2))
     c		/(dsqrt(2.d0*pi)*smear)
      enddo
      call DRG20r(ethresh,emax,nintegprec,u1rd,dxs_spin32_smear)
      end function dxs_spin32_smear
!     -------------------------------------------------------------------------
!     KLL smear
!     -------------------------------------------------------------------------
      real(8) function kll32_smear(egam,tman,xGamr,xMassr)
      use math_constants
      use hadron_properties
      use storage
      use size
      implicit none
      real(8),dimension(20000) :: enerd
      real(8),dimension(20000) :: u1rd
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
!     Internal Variables
      real(8) :: gaussener,gausstman,smear,xcos,ethresh
      integer :: i1,n1p
      real(8) :: a0,ap,Anorm,b0,st,xBrn
!     External Variables
      real (8),external :: kll32
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      smear = choice(4)
      ethresh=(xmjpsi**2+2.d0*xmpro*xmjpsi)/(2.d0*xmpro)
      call dsg20r(ethresh,emax,nintegprec,enerd,n1p)
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      do i1=1,n1p
         gaussener=enerd(i1)
         gausstman=(-2.d0*gaussener**2*xmpro + xmjpsi**2*xmpro + 
     -    	gaussener*(xmjpsi**2 + xcos*dsqrt(xmpro*
     -		(2.d0*gaussener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*gaussener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(gaussener + xmpro))/
     -		(xmpro*(2.d0*gaussener + xmpro)))))/
     -  	(2.d0*gaussener + xmpro)
         u1rd(i1)=  kll32(gaussener,gausstman-.0000000001
     -	,xGamr,xMassr)*
     c		dexp(-(egam-gaussener)**2/(2.d0*smear**2))
     c		/(dsqrt(2.d0*pi)*smear)
      enddo
      call DRG20r(ethresh,emax,nintegprec,u1rd,kll32_smear)
      end function kll32_smear
!     -------------------------------------------------------------------------
!     SDM rho00 spin-3/2
!     -------------------------------------------------------------------------
      complex(8) function rho0032(egam,tman,xGamr,xMassr)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
      real(8) :: a0,ap,Anorm,b0,st,xBrn
      real(8) :: a12flag,a32flag
!     Internal Variables
      real(8) :: dxsnorm,dxs_spin32temp,dfunc,deltafunc,jr,xcos
      real(8) :: temp
      real(8), dimension (24) :: listdfunc,listdelfunc
      integer :: i1,angmom
!     External Variables
      complex (8),external :: amp
      real (8), external :: d3212p32p,d3212p12p,d3212p12m,d3212p32m
      real (8), external :: d3212m32p,d3212m12p,d3212m12m,d3212m32m
      real (8), external :: d3232p32p,d3232p12p,d3232p12m,d3232p32m
      real (8), external :: d3232m32p,d3232m12p,d3232m12m,d3232m32m
      jr = 3.d0/2.d0
      angmom = 0
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      a12flag    = choice(1)
      a32flag    = choice(2)
! Note that here I use a different order for the functions than usual! It makes it easier to follow compared to Mathematica.
      listdelfunc(1)  = 1.d0
      listdelfunc(2)  = 0.d0
      listdelfunc(3)  = 0.d0
      listdelfunc(4)  = 0.d0
      listdelfunc(5)  = 0.d0
      listdelfunc(6)  = 0.d0
      listdelfunc(7)  = 0.d0
      listdelfunc(8)  = 1.d0
      listdelfunc(9)  = 0.d0
      listdelfunc(10) = 0.d0
      listdelfunc(11) = 0.d0
      listdelfunc(12) = 0.d0
      listdelfunc(13) = 0.d0
      listdelfunc(14) = 0.d0
      listdelfunc(15) = 0.d0
      listdelfunc(16) = 0.d0
      listdelfunc(17) = 1.d0
      listdelfunc(18) = 0.d0
      listdelfunc(19) = 0.d0
      listdelfunc(20) = 0.d0
      listdelfunc(21) = 0.d0
      listdelfunc(22) = 0.d0
      listdelfunc(23) = 0.d0
      listdelfunc(24) = 1.d0
      listdfunc(1)  = d3212p12p(xcos)*a12flag
      listdfunc(2)  = d3212p32p(xcos)*a12flag
      listdfunc(3)  = d3212p12m(xcos)*a12flag
      listdfunc(4)  = d3212p12p(xcos)*a12flag
      listdfunc(5)  = d3212p32m(xcos)*a12flag
      listdfunc(6)  = d3212p12m(xcos)*a12flag
      listdfunc(7)  = d3232p12p(xcos)*a32flag
      listdfunc(8)  = d3232p32p(xcos)*a32flag
      listdfunc(9)  = d3232p12m(xcos)*a32flag
      listdfunc(10) = d3232p12p(xcos)*a32flag
      listdfunc(11) = d3232p32m(xcos)*a32flag
      listdfunc(12) = d3232p12m(xcos)*a32flag
      listdfunc(13) = d3232m12p(xcos)*a32flag
      listdfunc(14) = d3232m32p(xcos)*a32flag
      listdfunc(15) = d3232m12m(xcos)*a32flag
      listdfunc(16) = d3232m12p(xcos)*a32flag
      listdfunc(17) = d3232m32m(xcos)*a32flag
      listdfunc(18) = d3232m12m(xcos)*a32flag
      listdfunc(19) = d3212m12p(xcos)*a12flag
      listdfunc(20) = d3212m32p(xcos)*a12flag
      listdfunc(21) = d3212m12m(xcos)*a12flag
      listdfunc(22) = d3212m12p(xcos)*a12flag
      listdfunc(23) = d3212m32m(xcos)*a12flag
      listdfunc(24) = d3212m12m(xcos)*a12flag
      rho0032 = (0.d0,0.d0)
      temp = 0.d0
      do i1=1,24
        temp=temp+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos))
      enddo
      do i1=1,4
        rho0032=rho0032+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-2),listdelfunc(6*i1-2),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-2),listdelfunc(6*i1-2),
     -	angmom,jr,xcos))+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-3),listdelfunc(6*i1-3),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-3),listdelfunc(6*i1-3),
     -	angmom,jr,xcos))
      enddo
      rho0032=rho0032/temp
      end function rho0032

!————————————————————————————————
!     SDM rhop0 spin-3/2
!     -------------------------------------------------------------------------
      complex(8) function rhop032(egam,tman,xGamr,xMassr)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
      real(8) :: a0,ap,Anorm,b0,st,xBrn
      real(8) :: a12flag,a32flag
!     Internal Variables
      real(8) :: dxsnorm,dxs_spin32temp,dfunc,deltafunc,jr,xcos
      real(8) :: temp
      real(8), dimension (24) :: listdfunc,listdelfunc
      integer :: i1,angmom
!     External Variables
      complex (8),external :: amp
      real (8), external :: d3212p32p,d3212p12p,d3212p12m,d3212p32m
      real (8), external :: d3212m32p,d3212m12p,d3212m12m,d3212m32m
      real (8), external :: d3232p32p,d3232p12p,d3232p12m,d3232p32m
      real (8), external :: d3232m32p,d3232m12p,d3232m12m,d3232m32m
      jr = 3.d0/2.d0
      angmom = 0
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      a12flag    = choice(1)
      a32flag    = choice(2)
! Note that here I use a different order for the functions than usual! It makes it easier to follow compared to Mathematica.
      listdelfunc(1)  = 1.d0
      listdelfunc(2)  = 0.d0
      listdelfunc(3)  = 0.d0
      listdelfunc(4)  = 0.d0
      listdelfunc(5)  = 0.d0
      listdelfunc(6)  = 0.d0
      listdelfunc(7)  = 0.d0
      listdelfunc(8)  = 1.d0
      listdelfunc(9)  = 0.d0
      listdelfunc(10) = 0.d0
      listdelfunc(11) = 0.d0
      listdelfunc(12) = 0.d0
      listdelfunc(13) = 0.d0
      listdelfunc(14) = 0.d0
      listdelfunc(15) = 0.d0
      listdelfunc(16) = 0.d0
      listdelfunc(17) = 1.d0
      listdelfunc(18) = 0.d0
      listdelfunc(19) = 0.d0
      listdelfunc(20) = 0.d0
      listdelfunc(21) = 0.d0
      listdelfunc(22) = 0.d0
      listdelfunc(23) = 0.d0
      listdelfunc(24) = 1.d0
      listdfunc(1)  = d3212p12p(xcos)*a12flag
      listdfunc(2)  = d3212p32p(xcos)*a12flag
      listdfunc(3)  = d3212p12m(xcos)*a12flag
      listdfunc(4)  = d3212p12p(xcos)*a12flag
      listdfunc(5)  = d3212p32m(xcos)*a12flag
      listdfunc(6)  = d3212p12m(xcos)*a12flag
      listdfunc(7)  = d3232p12p(xcos)*a32flag
      listdfunc(8)  = d3232p32p(xcos)*a32flag
      listdfunc(9)  = d3232p12m(xcos)*a32flag
      listdfunc(10) = d3232p12p(xcos)*a32flag
      listdfunc(11) = d3232p32m(xcos)*a32flag
      listdfunc(12) = d3232p12m(xcos)*a32flag
      listdfunc(13) = d3232m12p(xcos)*a32flag
      listdfunc(14) = d3232m32p(xcos)*a32flag
      listdfunc(15) = d3232m12m(xcos)*a32flag
      listdfunc(16) = d3232m12p(xcos)*a32flag
      listdfunc(17) = d3232m32m(xcos)*a32flag
      listdfunc(18) = d3232m12m(xcos)*a32flag
      listdfunc(19) = d3212m12p(xcos)*a12flag
      listdfunc(20) = d3212m32p(xcos)*a12flag
      listdfunc(21) = d3212m12m(xcos)*a12flag
      listdfunc(22) = d3212m12p(xcos)*a12flag
      listdfunc(23) = d3212m32m(xcos)*a12flag
      listdfunc(24) = d3212m12m(xcos)*a12flag
      rhop032 = (0.d0,0.d0)
      temp = 0.d0
      do i1=1,24
        temp=temp+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos))
      enddo
      do i1=1,4
        rhop032=rhop032+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-5),listdelfunc(6*i1-5),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-3),listdelfunc(6*i1-3),
     -	angmom,jr,xcos))+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-4),listdelfunc(6*i1-4),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-2),listdelfunc(6*i1-2),
     -	angmom,jr,xcos))
      enddo
      rhop032=rhop032/temp
      end function rhop032
!-------------------------------------------------------------------------
!     SDM rhopm spin-3/2
!     -------------------------------------------------------------------------
      complex(8) function rhopm32(egam,tman,xGamr,xMassr)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
      real(8) :: a0,ap,Anorm,b0,st,xBrn
      real(8) :: a12flag,a32flag
!     Internal Variables
      real(8) :: dxsnorm,dxs_spin32temp,dfunc,deltafunc,jr,xcos
      real(8) :: temp
      real(8), dimension (24) :: listdfunc,listdelfunc
      integer :: i1,angmom
!     External Variables
      complex (8),external :: amp
      real (8), external :: d3212p32p,d3212p12p,d3212p12m,d3212p32m
      real (8), external :: d3212m32p,d3212m12p,d3212m12m,d3212m32m
      real (8), external :: d3232p32p,d3232p12p,d3232p12m,d3232p32m
      real (8), external :: d3232m32p,d3232m12p,d3232m12m,d3232m32m
      jr = 3.d0/2.d0
      angmom = 0
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      a12flag    = choice(1)
      a32flag    = choice(2)
! Note that here I use a different order for the functions than usual! It makes it easier to follow compared to Mathematica.
      listdelfunc(1)  = 1.d0
      listdelfunc(2)  = 0.d0
      listdelfunc(3)  = 0.d0
      listdelfunc(4)  = 0.d0
      listdelfunc(5)  = 0.d0
      listdelfunc(6)  = 0.d0
      listdelfunc(7)  = 0.d0
      listdelfunc(8)  = 1.d0
      listdelfunc(9)  = 0.d0
      listdelfunc(10) = 0.d0
      listdelfunc(11) = 0.d0
      listdelfunc(12) = 0.d0
      listdelfunc(13) = 0.d0
      listdelfunc(14) = 0.d0
      listdelfunc(15) = 0.d0
      listdelfunc(16) = 0.d0
      listdelfunc(17) = 1.d0
      listdelfunc(18) = 0.d0
      listdelfunc(19) = 0.d0
      listdelfunc(20) = 0.d0
      listdelfunc(21) = 0.d0
      listdelfunc(22) = 0.d0
      listdelfunc(23) = 0.d0
      listdelfunc(24) = 1.d0
      listdfunc(1)  = d3212p12p(xcos)*a12flag
      listdfunc(2)  = d3212p32p(xcos)*a12flag
      listdfunc(3)  = d3212p12m(xcos)*a12flag
      listdfunc(4)  = d3212p12p(xcos)*a12flag
      listdfunc(5)  = d3212p32m(xcos)*a12flag
      listdfunc(6)  = d3212p12m(xcos)*a12flag
      listdfunc(7)  = d3232p12p(xcos)*a32flag
      listdfunc(8)  = d3232p32p(xcos)*a32flag
      listdfunc(9)  = d3232p12m(xcos)*a32flag
      listdfunc(10) = d3232p12p(xcos)*a32flag
      listdfunc(11) = d3232p32m(xcos)*a32flag
      listdfunc(12) = d3232p12m(xcos)*a32flag
      listdfunc(13) = d3232m12p(xcos)*a32flag
      listdfunc(14) = d3232m32p(xcos)*a32flag
      listdfunc(15) = d3232m12m(xcos)*a32flag
      listdfunc(16) = d3232m12p(xcos)*a32flag
      listdfunc(17) = d3232m32m(xcos)*a32flag
      listdfunc(18) = d3232m12m(xcos)*a32flag
      listdfunc(19) = d3212m12p(xcos)*a12flag
      listdfunc(20) = d3212m32p(xcos)*a12flag
      listdfunc(21) = d3212m12m(xcos)*a12flag
      listdfunc(22) = d3212m12p(xcos)*a12flag
      listdfunc(23) = d3212m32m(xcos)*a12flag
      listdfunc(24) = d3212m12m(xcos)*a12flag
      rhopm32 = (0.d0,0.d0)
      temp = 0.d0
      do i1=1,24
        temp=temp+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos))
      enddo
      do i1=1,4
        rhopm32=rhopm32+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-5),listdelfunc(6*i1-5),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-1),listdelfunc(6*i1-1),
     -	angmom,jr,xcos))+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-4),listdelfunc(6*i1-4),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1),listdelfunc(6*i1),
     -	angmom,jr,xcos))
      enddo
      rhopm32=rhopm32/temp
      end function rhopm32
!————————————————————————————————————
!     rho00 smear
!     -------------------------------------------------------------------------
      complex(8) function rho0032_smear(egam,tman,xGamr,xMassr)
      use math_constants
      use hadron_properties
      use storage
      use size
      implicit none
      real(8),dimension(20000) :: enerd
      complex(8),dimension(20000) :: u1rd
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
!     Internal Variables
      real(8) :: gaussener,gausstman,smear,xcos,ethresh
      integer :: i1,n1p
      real(8) :: a0,ap,Anorm,b0,st,xBrn
!     External Variables
      complex (8),external :: rho0032
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      smear = choice(4)
      ethresh=(xmjpsi**2+2.d0*xmpro*xmjpsi)/(2.d0*xmpro)
      call dsg20r(ethresh,emax,nintegprec,enerd,n1p)
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      do i1=1,n1p
         gaussener=enerd(i1)
         gausstman=(-2.d0*gaussener**2*xmpro + xmjpsi**2*xmpro + 
     -    	gaussener*(xmjpsi**2 + xcos*dsqrt(xmpro*
     -		(2.d0*gaussener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*gaussener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(gaussener + xmpro))/
     -		(xmpro*(2.d0*gaussener + xmpro)))))/
     -  	(2.d0*gaussener + xmpro)
         u1rd(i1)=  rho0032(gaussener,gausstman-.0000000001
     -	,xGamr,xMassr)*
     c		dexp(-(egam-gaussener)**2/(2.d0*smear**2))
     c		/(dsqrt(2.d0*pi)*smear)
      enddo
      call DRG20c(ethresh,emax,nintegprec,u1rd,rho0032_smear)
      end function rho0032_smear
!-------------------------------------------------------------------------
!     rhop0 smear
!     -------------------------------------------------------------------------
      complex(8) function rhop032_smear(egam,tman,xGamr,xMassr)
      use math_constants
      use hadron_properties
      use storage
      use size
      implicit none
      real(8),dimension(20000) :: enerd
      complex(8),dimension(20000) :: u1rd
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
!     Internal Variables
      real(8) :: gaussener,gausstman,smear,xcos,ethresh
      integer :: i1,n1p
      real(8) :: a0,ap,Anorm,b0,st,xBrn
!     External Variables
      complex (8),external :: rhop032
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      smear = choice(4)
      ethresh=(xmjpsi**2+2.d0*xmpro*xmjpsi)/(2.d0*xmpro)
      call dsg20r(ethresh,emax,nintegprec,enerd,n1p)
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      do i1=1,n1p
         gaussener=enerd(i1)
         gausstman=(-2.d0*gaussener**2*xmpro + xmjpsi**2*xmpro + 
     -    	gaussener*(xmjpsi**2 + xcos*dsqrt(xmpro*
     -		(2.d0*gaussener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*gaussener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(gaussener + xmpro))/
     -		(xmpro*(2.d0*gaussener + xmpro)))))/
     -  	(2.d0*gaussener + xmpro)
         u1rd(i1)=  rhop032(gaussener,gausstman-.0000000001
     -	,xGamr,xMassr)*
     c		dexp(-(egam-gaussener)**2/(2.d0*smear**2))
     c		/(dsqrt(2.d0*pi)*smear)
      enddo
      call DRG20c(ethresh,emax,nintegprec,u1rd,rhop032_smear)
      end function rhop032_smear


!————————————————————————————————————
!     rhopm smear
!     -------------------------------------------------------------------------
      complex(8) function rhopm32_smear(egam,tman,xGamr,xMassr)
      use math_constants
      use hadron_properties
      use storage
      use size
      implicit none
      real(8),dimension(20000) :: enerd
      complex(8),dimension(20000) :: u1rd
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
!     Internal Variables
      real(8) :: gaussener,gausstman,smear,xcos,ethresh
      integer :: i1,n1p
      real(8) :: a0,ap,Anorm,b0,st,xBrn
!     External Variables
      complex (8),external :: rhopm32
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      smear = choice(4)
      ethresh=(xmjpsi**2+2.d0*xmpro*xmjpsi)/(2.d0*xmpro)
      call dsg20r(ethresh,emax,nintegprec,enerd,n1p)
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      do i1=1,n1p
         gaussener=enerd(i1)
         gausstman=(-2.d0*gaussener**2*xmpro + xmjpsi**2*xmpro + 
     -    	gaussener*(xmjpsi**2 + xcos*dsqrt(xmpro*
     -		(2.d0*gaussener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*gaussener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(gaussener + xmpro))/
     -		(xmpro*(2.d0*gaussener + xmpro)))))/
     -  	(2.d0*gaussener + xmpro)
         u1rd(i1)=  rhopm32(gaussener,gausstman-.0000000001
     -	,xGamr,xMassr)*
     c		dexp(-(egam-gaussener)**2/(2.d0*smear**2))
     c		/(dsqrt(2.d0*pi)*smear)
      enddo
      call DRG20c(ethresh,emax,nintegprec,u1rd,rhopm32_smear)
      end function rhopm32_smear

!————————————————————————————————————
!     spin-5/2 resonance model
!     -------------------------------------------------------------------------
!     Differential
!     -------------------------------------------------------------------------
      real(8) function dxs_spin52(egam,tman,xGamr,xMassr)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,tmannew,xGamr,xMassr
      real(8) :: a0,ap,Anorm,b0,st,xBrn
      real(8) :: a12flag,a32flag
!     Internal Variables
      real(8) :: dxsnorm,dfunc,deltafunc,jr,xcos
      real(8), dimension (24) :: listdfunc,listdelfunc
      integer :: i,angmom
!     External Variables
      real (8), external :: ampsqr
      real (8), external :: d5212p32p,d5212p12p,d5212p12m,d5212p32m
      real (8), external :: d5212m32p,d5212m12p,d5212m12m,d5212m32m
      real (8), external :: d5232p32p,d5232p12p,d5232p12m,d5232p32m
      real (8), external :: d5232m32p,d5232m12p,d5232m12m,d5232m32m
      dxsnorm =  (6084.375d0*alpha)/(egam**2*xMpro**2)
      angmom = 1
      jr = 5.d0/2.d0
      if (tman.gt.0.d0) then
         xcos = 1.d0
         tmannew =  (-2.d0*egam**2*xmpro + xmjpsi**2*xmpro + 
     -    egam*(xmjpsi**2 + xcos*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -        dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2 - 
     -            4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -	(xmpro*(2.d0*egam + xmpro)))))/
     -  (2.d0*egam + xmpro)
      else
        xcos =   (2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      a12flag    = choice(1)
      a32flag    = choice(2)
      listdelfunc(1)  = 0.d0
      listdelfunc(2)  = 1.d0
      listdelfunc(3)  = 0.d0
      listdelfunc(4)  = 0.d0
      listdelfunc(5)  = 0.d0
      listdelfunc(6)  = 0.d0
      listdelfunc(7)  = 0.d0
      listdelfunc(8)  = 0.d0
      listdelfunc(9)  = 0.d0
      listdelfunc(10) = 1.d0
      listdelfunc(11) = 0.d0
      listdelfunc(12) = 0.d0
      listdelfunc(13) = 1.d0
      listdelfunc(14) = 0.d0
      listdelfunc(15) = 0.d0
      listdelfunc(16) = 0.d0
      listdelfunc(17) = 0.d0
      listdelfunc(18) = 0.d0
      listdelfunc(19) = 0.d0
      listdelfunc(20) = 0.d0
      listdelfunc(21) = 0.d0
      listdelfunc(22) = 0.d0
      listdelfunc(23) = 0.d0
      listdelfunc(24) = 1.d0
      listdfunc(1)  = d5212p32p(xcos)*a12flag
      listdfunc(2)  = d5212p12p(xcos)*a12flag
      listdfunc(3)  = d5212p12p(xcos)*a12flag
      listdfunc(4)  = d5212p12m(xcos)*a12flag
      listdfunc(5)  = d5212p12m(xcos)*a12flag
      listdfunc(6)  = d5212p32m(xcos)*a12flag
      listdfunc(7)  = d5212m32p(xcos)*a12flag
      listdfunc(8)  = d5212m12p(xcos)*a12flag
      listdfunc(9)  = d5212m12p(xcos)*a12flag
      listdfunc(10) = d5212m12m(xcos)*a12flag
      listdfunc(11) = d5212m12m(xcos)*a12flag
      listdfunc(12) = d5212m32m(xcos)*a12flag
      listdfunc(13) = d5232p32p(xcos)*a32flag
      listdfunc(14) = d5232p12p(xcos)*a32flag
      listdfunc(15) = d5232p12p(xcos)*a32flag
      listdfunc(16) = d5232p12m(xcos)*a32flag
      listdfunc(17) = d5232p12m(xcos)*a32flag
      listdfunc(18) = d5232p32m(xcos)*a32flag
      listdfunc(19) = d5232m32p(xcos)*a32flag
      listdfunc(20) = d5232m12p(xcos)*a32flag
      listdfunc(21) = d5232m12p(xcos)*a32flag
      listdfunc(22) = d5232m12m(xcos)*a32flag
      listdfunc(23) = d5232m12m(xcos)*a32flag
      listdfunc(24) = d5232m32m(xcos)*a32flag
!     Differential cross section calculation
      dxs_spin52 = 0.d0
      if (tman .lt. 0.d0) then
         do i = 1,24
            dxs_spin52 = 
     -	dxs_spin52 + dxsnorm*
     -	ampsqr(egam,tman,xGamr,xMassr,listdfunc(i)
     -	,listdelfunc(i),angmom,jr,xcos)
         enddo
      else
         do i = 1,24
            dxs_spin52 = 
     -	dxs_spin52 + dxsnorm*
     -	ampsqr(egam,tmannew,xGamr,xMassr,listdfunc(i)
     -	,listdelfunc(i),angmom,jr,xcos)
         enddo
      endif
      end function dxs_spin52
!-------------------------------------------------------------------------
!     KLL
!     -------------------------------------------------------------------------
      real(8) function kll52(egam,tman,xGamr,xMassr)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,tmannew,xGamr,xMassr
      real(8) :: a0,ap,Anorm,b0,st,xBrn
!     Internal Variables
!     External Variables
      real (8), external :: ampsqr52pp,ampsqr52pm
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
!     Calculation
      kll52 = (ampsqr52pp(egam,tman,xGamr,xMassr)-
     -  ampsqr52pm(egam,tman,xGamr,xMassr))/
     -  (ampsqr52pp(egam,tman,xGamr,xMassr)+
     -  ampsqr52pm(egam,tman,xGamr,xMassr))
      end function kll52
!-------------------------------------------------------------------------
!     Differential smear
!-------------------------------------------------------------------------
      real(8) function dxs_spin52_smear(egam,tman,xGamr,xMassr)
      use math_constants
      use hadron_properties
      use storage
      use size
      implicit none
      real(8),dimension(20000) :: enerd
      real(8),dimension(20000) :: u1rd
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
!     Internal Variables
      real(8) :: gaussener,gausstman,smear,xcos,ethresh
      integer :: i1,n1p
      real(8) :: a0,ap,Anorm,b0,st,xBrn
!     External Variables
      real (8),external :: dxs_spin52
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      smear = choice(4)
      ethresh=(xmjpsi**2+2.d0*xmpro*xmjpsi)/(2.d0*xmpro)
      call dsg20r(ethresh,emax,nintegprec,enerd,n1p)
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      do i1=1,n1p
         gaussener=enerd(i1)
         gausstman=(-2.d0*gaussener**2*xmpro + xmjpsi**2*xmpro + 
     -    	gaussener*(xmjpsi**2 + xcos*dsqrt(xmpro*
     -		(2.d0*gaussener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*gaussener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(gaussener + xmpro))/
     -		(xmpro*(2.d0*gaussener + xmpro)))))/
     -  	(2.d0*gaussener + xmpro)
         u1rd(i1)=  dxs_spin52(gaussener,gausstman-.0000000001
     -	,xGamr,xMassr)*
     c		dexp(-(egam-gaussener)**2/(2.d0*smear**2))
     c		/(dsqrt(2.d0*pi)*smear)
      enddo
      call DRG20r(ethresh,emax,nintegprec,u1rd,dxs_spin52_smear)
      end function dxs_spin52_smear
!     -------------------------------------------------------------------------
!     KLL smear
!     -------------------------------------------------------------------------
      real(8) function kll52_smear(egam,tman,xGamr,xMassr)
      use math_constants
      use hadron_properties
      use storage
      use size
      implicit none
      real(8),dimension(20000) :: enerd
      real(8),dimension(20000) :: u1rd
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
!     Internal Variables
      real(8) :: gaussener,gausstman,smear,xcos,ethresh
      integer :: i1,n1p
      real(8) :: a0,ap,Anorm,b0,st,xBrn
!     External Variables
      real (8),external :: kll52
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      smear = choice(4)
      ethresh=(xmjpsi**2+2.d0*xmpro*xmjpsi)/(2.d0*xmpro)
      call dsg20r(ethresh,emax,nintegprec,enerd,n1p)
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      do i1=1,n1p
         gaussener=enerd(i1)
         gausstman=(-2.d0*gaussener**2*xmpro + xmjpsi**2*xmpro + 
     -    	gaussener*(xmjpsi**2 + xcos*dsqrt(xmpro*
     -		(2.d0*gaussener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*gaussener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(gaussener + xmpro))/
     -		(xmpro*(2.d0*gaussener + xmpro)))))/
     -  	(2.d0*gaussener + xmpro)
         u1rd(i1)=  kll52(gaussener,gausstman-.0000000001
     -	,xGamr,xMassr)*
     c		dexp(-(egam-gaussener)**2/(2.d0*smear**2))
     c		/(dsqrt(2.d0*pi)*smear)
      enddo
      call DRG20r(ethresh,emax,nintegprec,u1rd,kll52_smear)
      end function kll52_smear



!     -------------------------------------------------------------------------
!     SDM rho00 spin-5/2
!     -------------------------------------------------------------------------
      complex(8) function rho0052(egam,tman,xGamr,xMassr)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
      real(8) :: a0,ap,Anorm,b0,st,xBrn
      real(8) :: a12flag,a32flag
!     Internal Variables
      real(8) :: dfunc,deltafunc,jr,xcos
      real(8) :: temp
      real(8), dimension (24) :: listdfunc,listdelfunc
      integer :: i1,angmom
!     External Variables
      complex (8),external :: amp
      real (8), external :: d5212p32p,d5212p12p,d5212p12m,d5212p32m
      real (8), external :: d5212m32p,d5212m12p,d5212m12m,d5212m32m
      real (8), external :: d5232p32p,d5232p12p,d5232p12m,d5232p32m
      real (8), external :: d5232m32p,d5232m12p,d5232m12m,d5232m32m
      jr = 5.d0/2.d0
      angmom = 1
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      a12flag    = choice(1)
      a32flag    = choice(2)
! Note that here I use a different order for the functions than usual! It makes it easier to follow compared to Mathematica.
      listdelfunc(1)  = 1.d0
      listdelfunc(2)  = 0.d0
      listdelfunc(3)  = 0.d0
      listdelfunc(4)  = 0.d0
      listdelfunc(5)  = 0.d0
      listdelfunc(6)  = 0.d0
      listdelfunc(7)  = 0.d0
      listdelfunc(8)  = 1.d0
      listdelfunc(9)  = 0.d0
      listdelfunc(10) = 0.d0
      listdelfunc(11) = 0.d0
      listdelfunc(12) = 0.d0
      listdelfunc(13) = 0.d0
      listdelfunc(14) = 0.d0
      listdelfunc(15) = 0.d0
      listdelfunc(16) = 0.d0
      listdelfunc(17) = 1.d0
      listdelfunc(18) = 0.d0
      listdelfunc(19) = 0.d0
      listdelfunc(20) = 0.d0
      listdelfunc(21) = 0.d0
      listdelfunc(22) = 0.d0
      listdelfunc(23) = 0.d0
      listdelfunc(24) = 1.d0
      listdfunc(1)  = d5212p12p(xcos)*a12flag
      listdfunc(2)  = d5212p32p(xcos)*a12flag
      listdfunc(3)  = d5212p12m(xcos)*a12flag
      listdfunc(4)  = d5212p12p(xcos)*a12flag
      listdfunc(5)  = d5212p32m(xcos)*a12flag
      listdfunc(6)  = d5212p12m(xcos)*a12flag
      listdfunc(7)  = d5232p12p(xcos)*a32flag
      listdfunc(8)  = d5232p32p(xcos)*a32flag
      listdfunc(9)  = d5232p12m(xcos)*a32flag
      listdfunc(10) = d5232p12p(xcos)*a32flag
      listdfunc(11) = d5232p32m(xcos)*a32flag
      listdfunc(12) = d5232p12m(xcos)*a32flag
      listdfunc(13) = d5232m12p(xcos)*a32flag
      listdfunc(14) = d5232m32p(xcos)*a32flag
      listdfunc(15) = d5232m12m(xcos)*a32flag
      listdfunc(16) = d5232m12p(xcos)*a32flag
      listdfunc(17) = d5232m32m(xcos)*a32flag
      listdfunc(18) = d5232m12m(xcos)*a32flag
      listdfunc(19) = d5212m12p(xcos)*a12flag
      listdfunc(20) = d5212m32p(xcos)*a12flag
      listdfunc(21) = d5212m12m(xcos)*a12flag
      listdfunc(22) = d5212m12p(xcos)*a12flag
      listdfunc(23) = d5212m32m(xcos)*a12flag
      listdfunc(24) = d5212m12m(xcos)*a12flag
      rho0052 = (0.d0,0.d0)
      temp = 0.d0
      do i1=1,24
        temp=temp+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos))
      enddo
      do i1=1,4
        rho0052=rho0052+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-3),listdelfunc(6*i1-3),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-3),listdelfunc(6*i1-3),
     -	angmom,jr,xcos))
        rho0052=rho0052+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-2),listdelfunc(6*i1-2),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-2),listdelfunc(6*i1-2),
     -	angmom,jr,xcos))
      enddo
      rho0052=rho0052/temp
      end function rho0052

!——————————————————————————————
!     SDM rhop0 spin-5/2
!     52-------------------------------------------------------------------------
      complex(8) function rhop052(egam,tman,xGamr,xMassr)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
      real(8) :: a0,ap,Anorm,b0,st,xBrn
      real(8) :: a12flag,a32flag
!     Internal Variables
      real(8) :: dxsnorm,dxs_spin32temp,dfunc,deltafunc,jr,xcos
      real(8) :: temp
      real(8), dimension (24) :: listdfunc,listdelfunc
      integer :: i1,angmom
!     External Variables
      complex (8),external :: amp
      real (8), external :: d5212p32p,d5212p12p,d5212p12m,d5212p32m
      real (8), external :: d5212m32p,d5212m12p,d5212m12m,d5212m32m
      real (8), external :: d5232p32p,d5232p12p,d5232p12m,d5232p32m
      real (8), external :: d5232m32p,d5232m12p,d5232m12m,d5232m32m
      jr = 5.d0/2.d0
      angmom = 1
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      a12flag    = choice(1)
      a32flag    = choice(2)
! Note that here I use a different order for the functions than usual! It makes it easier to follow compared to Mathematica.
      listdelfunc(1)  = 1.d0
      listdelfunc(2)  = 0.d0
      listdelfunc(3)  = 0.d0
      listdelfunc(4)  = 0.d0
      listdelfunc(5)  = 0.d0
      listdelfunc(6)  = 0.d0
      listdelfunc(7)  = 0.d0
      listdelfunc(8)  = 1.d0
      listdelfunc(9)  = 0.d0
      listdelfunc(10) = 0.d0
      listdelfunc(11) = 0.d0
      listdelfunc(12) = 0.d0
      listdelfunc(13) = 0.d0
      listdelfunc(14) = 0.d0
      listdelfunc(15) = 0.d0
      listdelfunc(16) = 0.d0
      listdelfunc(17) = 1.d0
      listdelfunc(18) = 0.d0
      listdelfunc(19) = 0.d0
      listdelfunc(20) = 0.d0
      listdelfunc(21) = 0.d0
      listdelfunc(22) = 0.d0
      listdelfunc(23) = 0.d0
      listdelfunc(24) = 1.d0
      listdfunc(1)  = d5212p12p(xcos)*a12flag
      listdfunc(2)  = d5212p32p(xcos)*a12flag
      listdfunc(3)  = d5212p12m(xcos)*a12flag
      listdfunc(4)  = d5212p12p(xcos)*a12flag
      listdfunc(5)  = d5212p32m(xcos)*a12flag
      listdfunc(6)  = d5212p12m(xcos)*a12flag
      listdfunc(7)  = d5232p12p(xcos)*a32flag
      listdfunc(8)  = d5232p32p(xcos)*a32flag
      listdfunc(9)  = d5232p12m(xcos)*a32flag
      listdfunc(10) = d5232p12p(xcos)*a32flag
      listdfunc(11) = d5232p32m(xcos)*a32flag
      listdfunc(12) = d5232p12m(xcos)*a32flag
      listdfunc(13) = d5232m12p(xcos)*a32flag
      listdfunc(14) = d5232m32p(xcos)*a32flag
      listdfunc(15) = d5232m12m(xcos)*a32flag
      listdfunc(16) = d5232m12p(xcos)*a32flag
      listdfunc(17) = d5232m32m(xcos)*a32flag
      listdfunc(18) = d5232m12m(xcos)*a32flag
      listdfunc(19) = d5212m12p(xcos)*a12flag
      listdfunc(20) = d5212m32p(xcos)*a12flag
      listdfunc(21) = d5212m12m(xcos)*a12flag
      listdfunc(22) = d5212m12p(xcos)*a12flag
      listdfunc(23) = d5212m32m(xcos)*a12flag
      listdfunc(24) = d5212m12m(xcos)*a12flag
      rhop052 = (0.d0,0.d0)
      temp = 0.d0
      do i1=1,24
        temp=temp+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos))
      enddo
      do i1=1,4
        rhop052=rhop052+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-5),listdelfunc(6*i1-5),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-3),listdelfunc(6*i1-3),
     -	angmom,jr,xcos))+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-4),listdelfunc(6*i1-4),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-2),listdelfunc(6*i1-2),
     -	angmom,jr,xcos))
      enddo
      rhop052=rhop052/temp
      end function rhop052
!-------------------------------------------------------------------------
!     SDM rhopm spin-5/2
!     -------------------------------------------------------------------------
      complex(8) function rhopm52(egam,tman,xGamr,xMassr)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
      real(8) :: a0,ap,Anorm,b0,st,xBrn
      real(8) :: a12flag,a32flag
!     Internal Variables
      real(8) :: dxsnorm,dxs_spin32temp,dfunc,deltafunc,jr,xcos
      real(8) :: temp
      real(8), dimension (24) :: listdfunc,listdelfunc
      integer :: i1,angmom
!     External Variables
      complex (8),external :: amp
      real (8), external :: d5212p32p,d5212p12p,d5212p12m,d5212p32m
      real (8), external :: d5212m32p,d5212m12p,d5212m12m,d5212m32m
      real (8), external :: d5232p32p,d5232p12p,d5232p12m,d5232p32m
      real (8), external :: d5232m32p,d5232m12p,d5232m12m,d5232m32m
      jr = 5.d0/2.d0
      angmom = 1
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      a12flag    = choice(1)
      a32flag    = choice(2)
! Note that here I use a different order for the functions than usual! It makes it easier to follow compared to Mathematica.
      listdelfunc(1)  = 1.d0
      listdelfunc(2)  = 0.d0
      listdelfunc(3)  = 0.d0
      listdelfunc(4)  = 0.d0
      listdelfunc(5)  = 0.d0
      listdelfunc(6)  = 0.d0
      listdelfunc(7)  = 0.d0
      listdelfunc(8)  = 1.d0
      listdelfunc(9)  = 0.d0
      listdelfunc(10) = 0.d0
      listdelfunc(11) = 0.d0
      listdelfunc(12) = 0.d0
      listdelfunc(13) = 0.d0
      listdelfunc(14) = 0.d0
      listdelfunc(15) = 0.d0
      listdelfunc(16) = 0.d0
      listdelfunc(17) = 1.d0
      listdelfunc(18) = 0.d0
      listdelfunc(19) = 0.d0
      listdelfunc(20) = 0.d0
      listdelfunc(21) = 0.d0
      listdelfunc(22) = 0.d0
      listdelfunc(23) = 0.d0
      listdelfunc(24) = 1.d0
      listdfunc(1)  = d5212p12p(xcos)*a12flag
      listdfunc(2)  = d5212p32p(xcos)*a12flag
      listdfunc(3)  = d5212p12m(xcos)*a12flag
      listdfunc(4)  = d5212p12p(xcos)*a12flag
      listdfunc(5)  = d5212p32m(xcos)*a12flag
      listdfunc(6)  = d5212p12m(xcos)*a12flag
      listdfunc(7)  = d5232p12p(xcos)*a32flag
      listdfunc(8)  = d5232p32p(xcos)*a32flag
      listdfunc(9)  = d5232p12m(xcos)*a32flag
      listdfunc(10) = d5232p12p(xcos)*a32flag
      listdfunc(11) = d5232p32m(xcos)*a32flag
      listdfunc(12) = d5232p12m(xcos)*a32flag
      listdfunc(13) = d5232m12p(xcos)*a32flag
      listdfunc(14) = d5232m32p(xcos)*a32flag
      listdfunc(15) = d5232m12m(xcos)*a32flag
      listdfunc(16) = d5232m12p(xcos)*a32flag
      listdfunc(17) = d5232m32m(xcos)*a32flag
      listdfunc(18) = d5232m12m(xcos)*a32flag
      listdfunc(19) = d5212m12p(xcos)*a12flag
      listdfunc(20) = d5212m32p(xcos)*a12flag
      listdfunc(21) = d5212m12m(xcos)*a12flag
      listdfunc(22) = d5212m12p(xcos)*a12flag
      listdfunc(23) = d5212m32m(xcos)*a12flag
      listdfunc(24) = d5212m12m(xcos)*a12flag
      rhopm52 = (0.d0,0.d0)
      temp = 0.d0
      do i1=1,24
        temp=temp+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos))
      enddo
      do i1=1,4
        rhopm52=rhopm52+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-5),listdelfunc(6*i1-5),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-1),listdelfunc(6*i1-1),
     -	angmom,jr,xcos))+amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1-4),listdelfunc(6*i1-4),
     -	angmom,jr,xcos)
     -  *conjg(amp(egam,tman,xGamr,xMassr,
     -  listdfunc(6*i1),listdelfunc(6*i1),
     -	angmom,jr,xcos))
      enddo
      rhopm52=rhopm52/temp
      end function rhopm52

!     -------------------------------------------------------------------------
!     rho00 smear
!     -------------------------------------------------------------------------
      complex(8) function rho0052_smear(egam,tman,xGamr,xMassr)
      use math_constants
      use hadron_properties
      use storage
      use size
      implicit none
      real(8),dimension(20000) :: enerd
      complex(8),dimension(20000) :: u1rd
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
!     Internal Variables
      real(8) :: gaussener,gausstman,smear,xcos,ethresh
      integer :: i1,n1p
      real(8) :: a0,ap,Anorm,b0,st,xBrn
!     External Variables
      complex (8),external :: rho0052
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      smear = choice(4)
      ethresh=(xmjpsi**2+2.d0*xmpro*xmjpsi)/(2.d0*xmpro)
      call dsg20r(ethresh,emax,nintegprec,enerd,n1p)
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      do i1=1,n1p
         gaussener=enerd(i1)
         gausstman=(-2.d0*gaussener**2*xmpro + xmjpsi**2*xmpro + 
     -    	gaussener*(xmjpsi**2 + xcos*dsqrt(xmpro*
     -		(2.d0*gaussener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*gaussener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(gaussener + xmpro))/
     -		(xmpro*(2.d0*gaussener + xmpro)))))/
     -  	(2.d0*gaussener + xmpro)
         u1rd(i1)=  rho0052(gaussener,gausstman-.0000000001
     -	,xGamr,xMassr)*
     c		dexp(-(egam-gaussener)**2/(2.d0*smear**2))
     c		/(dsqrt(2.d0*pi)*smear)
      enddo
      call DRG20c(ethresh,emax,nintegprec,u1rd,rho0052_smear)
      end function rho0052_smear
!     


!————————————————————————————————————
!     rhop0 smear
!     -------------------------------------------------------------------------
      complex(8) function rhop052_smear(egam,tman,xGamr,xMassr)
      use math_constants
      use hadron_properties
      use storage
      use size
      implicit none
      real(8),dimension(20000) :: enerd
      complex(8),dimension(20000) :: u1rd
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
!     Internal Variables
      real(8) :: gaussener,gausstman,smear,xcos,ethresh
      integer :: i1,n1p
      real(8) :: a0,ap,Anorm,b0,st,xBrn
!     External Variables
      complex (8),external :: rhop052
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      smear = choice(4)
      ethresh=(xmjpsi**2+2.d0*xmpro*xmjpsi)/(2.d0*xmpro)
      call dsg20r(ethresh,emax,nintegprec,enerd,n1p)
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      do i1=1,n1p
         gaussener=enerd(i1)
         gausstman=(-2.d0*gaussener**2*xmpro + xmjpsi**2*xmpro + 
     -    	gaussener*(xmjpsi**2 + xcos*dsqrt(xmpro*
     -		(2.d0*gaussener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*gaussener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(gaussener + xmpro))/
     -		(xmpro*(2.d0*gaussener + xmpro)))))/
     -  	(2.d0*gaussener + xmpro)
         u1rd(i1)=  rhop052(gaussener,gausstman-.0000000001
     -	,xGamr,xMassr)*
     c		dexp(-(egam-gaussener)**2/(2.d0*smear**2))
     c		/(dsqrt(2.d0*pi)*smear)
      enddo
      call DRG20c(ethresh,emax,nintegprec,u1rd,rhop052_smear)
      end function rhop052_smear


!————————————————————————————————————
!     rhopm smear
!     -------------------------------------------------------------------------
      complex(8) function rhopm52_smear(egam,tman,xGamr,xMassr)
      use math_constants
      use hadron_properties
      use storage
      use size
      implicit none
      real(8),dimension(20000) :: enerd
      complex(8),dimension(20000) :: u1rd
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
!     Internal Variables
      real(8) :: gaussener,gausstman,smear,xcos,ethresh
      integer :: i1,n1p
      real(8) :: a0,ap,Anorm,b0,st,xBrn
!     External Variables
      complex (8),external :: rhopm52
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      smear = choice(4)
      ethresh=(xmjpsi**2+2.d0*xmpro*xmjpsi)/(2.d0*xmpro)
      call dsg20r(ethresh,emax,nintegprec,enerd,n1p)
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      do i1=1,n1p
         gaussener=enerd(i1)
         gausstman=(-2.d0*gaussener**2*xmpro + xmjpsi**2*xmpro + 
     -    	gaussener*(xmjpsi**2 + xcos*dsqrt(xmpro*
     -		(2.d0*gaussener + xmpro))
     -        	*dsqrt((xmjpsi**4 + 4.d0*gaussener**2*xmpro**2 - 
     -          4.d0*xmjpsi**2*xmpro*(gaussener + xmpro))/
     -		(xmpro*(2.d0*gaussener + xmpro)))))/
     -  	(2.d0*gaussener + xmpro)
         u1rd(i1)=  rhopm52(gaussener,gausstman-.0000000001
     -	,xGamr,xMassr)*
     c		dexp(-(egam-gaussener)**2/(2.d0*smear**2))
     c		/(dsqrt(2.d0*pi)*smear)
      enddo
      call DRG20c(ethresh,emax,nintegprec,u1rd,rhopm52_smear)
      end function rhopm52_smear
!————————————————————————————————————
!     -------------------------------------------------------------------------
!     Total
!     -------------------------------------------------------------------------
!     -------------------------------------------------------------------------
      real(8) function totxs(egam,xGamr,xMassr)
      use hadron_properties
      use storage
      use size
      implicit none
      real(8),dimension(20000) :: costhd
      real(8),dimension(20000) :: u1rd
!     Input Variables
      real(8) :: egam,xGamr,xMassr
!     Internal Variables
      real(8) :: xcos,tman,dtdcos
      integer :: i1,n1p
      real(8) :: a0,ap,Anorm,b0,st,xBrn,smear
      character(5) :: inputtype
!     External Variables
      real (8),external :: dxs_spin52,dxs_spin52_smear,
     c	dxs_spin32,dxs_spin32_smear
      inputtype    = strings(1)
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      smear  = choice(4)
      call dsg20r(-1.d0,1.d0,ntotintprec,costhd,n1p)
      do i1=1,n1p
         xcos=costhd(i1)
         tman =  (-2.d0*egam**2*xmpro + xmjpsi**2*xmpro + 
     -    egam*(xmjpsi**2 + xcos*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -        dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2 - 
     -            4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -	(xmpro*(2.d0*egam + xmpro)))))/
     -  (2.d0*egam + xmpro)
         dtdcos =  (egam*(dsqrt(xmpro*(2.d0*egam + xmpro))*
     -        dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2 - 
     -            4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -	(xmpro*(2.d0*egam + xmpro)))))/
     -  (2.d0*egam + xmpro)
         if (inputtype.eq.'52'.and.smear.eq.0.d0) then
            u1rd(i1)=dtdcos*dxs_spin52(egam,tman,xGamr,xMassr)
         elseif (inputtype.eq.'52'.and.smear.gt.0.d0) then
            u1rd(i1)=dtdcos*dxs_spin52_smear(egam,tman,xGamr,xMassr)
         elseif (inputtype.eq.'32'.and.smear.eq.0.d0) then
            u1rd(i1)=dtdcos*dxs_spin32(egam,tman,xGamr,xMassr)
         elseif (inputtype.eq.'32'.and.smear.gt.0.d0) then
            u1rd(i1)=dtdcos*dxs_spin32_smear(egam,tman,xGamr,xMassr)
         endif
      enddo
      call DRG20r(-1.d0,1.d0,ntotintprec,u1rd,totxs)
      end function totxs
!     -------------------------------------------------------------------------
!     amplitude definition
!     -------------------------------------------------------------------------
      complex(8) function amp(egam,tman,xGamr,xMassr,dfunc,
     -  deltafunc,
     -	angmom,jr,xcos)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr,dfunc,deltafunc,jr
!     Internal Variables
      real(8) :: a0,ap,Anorm,b0,st,xBrn,xcos
      integer :: angmom
!     Parameter assignment
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
!     Unsquared amplitude
      amp = (0,1)*anorm*deltafunc*
     -   dexp
     -    ((b0*(-1.d0 + xcos)*egam*xMpro*
     -        dSqrt((xMjpsi**4 + 
     -            4.d0*egam**2*xMpro**2 - 
     -            4.d0*xMjpsi**2*xMpro*(egam + xMpro))/
     -          (xMpro*(2.d0*egam + xMpro))))/
     -      dSqrt(xMpro*(2.d0*egam + xMpro)))*
     -   ((1,0)*(-st + xMpro*(2.d0*egam + xMpro)))**
     -    (a0 + (ap*(egam*xMjpsi**2 - 
     -           2.d0*egam**2*xMpro + 
     -           xMjpsi**2*xMpro + 
     -           xcos*egam*
     -            dSqrt(xMpro*(2.d0*egam + xMpro))*
     -            dSqrt((xMjpsi**4 + 
     -                4.d0*egam**2*xMpro**2 - 
     -                4.d0*xMjpsi**2*xMpro*
     -                 (egam + xMpro))/
     -              (xMpro*(2.d0*egam + xMpro)))))/
     -       (2.d0*egam + xMpro)) + 
     -  (2.d0*alpha*dsqrt(2.d0)*dfunc*egam*Pi*
     -     dSqrt((1.d0 + 2.d0*jr)*xBrn*xGamr)*xMpro**2*
     -     (-xMassr**2 + xMpro**2)*
     -     dSqrt((4.d0**(1 - angmom)*alpha*(1.d0 + 2.d0*jr)*
     -         xBrn*xFjpsi**2*xGamr*xMassr**3*
     -         ((xMassr**2 - xMpro**2)/
     -            (xMassr*
     -              dSqrt(-xMjpsi**2 + 
     -                (xMassr**2 + xMjpsi**2 - 
     -                  xMpro**2)**2/(4.d0*xMassr**2))
     -              ))**(1.d0 + 2.d0*angmom))/
     -       (xMjpsi**2*xMpro*
     -         (xMassr**2 - xMpro**2)**2)))/
     -   (3.d0*dSqrt(xMassr*xMpro)*
     -     dsqrt(((alpha*egam*xMpro)/
     -        dSqrt(xMpro*(2.d0*egam + xMpro)))**3)*
     -     dsqrt(dsqrt(((xMjpsi**4 + 4.d0*egam**2*xMpro**2 - 
     -          4.d0*xMjpsi**2*xMpro*(egam + xMpro))/
     -        (xMpro*(2.d0*egam + xMpro)))))*
     -     ((0,1)*xGamr*xMassr 
     -       -xMassr**2 + 
     -          xMpro*(2.d0*egam + xMpro)))
      end function amp

!     -------------------------------------------------------------------------
!     amplitude square definition
!     -------------------------------------------------------------------------
      real(8) function ampsqr(egam,tman,xGamr,xMassr,dfunc,deltafunc,
     -	angmom,jr,xcos)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr,dfunc,deltafunc,jr
!     Internal Variables
      real(8) :: a0,ap,Anorm,b0,st,xBrn,xcos
!     External Variables
      complex (8),external :: amp
!     The amplitude is complex
      integer :: angmom
!     Parameter assignment
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
!     Squared amplitude calculation
      ampsqr = real(amp(egam,tman,xGamr,xMassr,dfunc,deltafunc,
     -	angmom,jr,xcos)*
     -  conjg(amp(egam,tman,xGamr,xMassr,dfunc,deltafunc,
     -	angmom,jr,xcos)))
      end function ampsqr

!     -------------------------------------------------------------------------
!     Sum of amplitudes squared for helicity non-flip, spin-3/2
!     -------------------------------------------------------------------------
      real(8) function ampsqr32pp(egam,tman,xGamr,xMassr)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,tmannew,xGamr,xMassr
      real(8) :: a0,ap,Anorm,b0,st,xBrn
      real(8) :: a12flag,a32flag
!     Internal Variables
      real(8) :: dxsnorm,dxs_spin32temp,dfunc,deltafunc,jr,xcos
      real(8), dimension (24) :: listdfunc,listdelfunc
      integer :: i1,angmom
!     External Variables
      real (8),external :: ampsqr
      real (8), external :: d3212p32p,d3212p12p,d3212p12m,d3212p32m
      real (8), external :: d3212m32p,d3212m12p,d3212m12m,d3212m32m
      real (8), external :: d3232p32p,d3232p12p,d3232p12m,d3232p32m
      real (8), external :: d3232m32p,d3232m12p,d3232m12m,d3232m32m
      jr = 3.d0/2.d0
      angmom = 0
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      a12flag    = choice(1)
      a32flag    = choice(2)
      listdelfunc(1)  = 0.d0
      listdelfunc(2)  = 1.d0
      listdelfunc(3)  = 0.d0
      listdelfunc(4)  = 0.d0
      listdelfunc(5)  = 0.d0
      listdelfunc(6)  = 0.d0
      listdelfunc(7)  = 0.d0
      listdelfunc(8)  = 0.d0
      listdelfunc(9)  = 0.d0
      listdelfunc(10) = 1.d0
      listdelfunc(11) = 0.d0
      listdelfunc(12) = 0.d0
      listdelfunc(13) = 1.d0
      listdelfunc(14) = 0.d0
      listdelfunc(15) = 0.d0
      listdelfunc(16) = 0.d0
      listdelfunc(17) = 0.d0
      listdelfunc(18) = 0.d0
      listdelfunc(19) = 0.d0
      listdelfunc(20) = 0.d0
      listdelfunc(21) = 0.d0
      listdelfunc(22) = 0.d0
      listdelfunc(23) = 0.d0
      listdelfunc(24) = 1.d0
      listdfunc(1)  = d3212p32p(xcos)*a12flag
      listdfunc(2)  = d3212p12p(xcos)*a12flag
      listdfunc(3)  = d3212p12p(xcos)*a12flag
      listdfunc(4)  = d3212p12m(xcos)*a12flag
      listdfunc(5)  = d3212p12m(xcos)*a12flag
      listdfunc(6)  = d3212p32m(xcos)*a12flag
      listdfunc(7)  = d3212m32p(xcos)*a12flag
      listdfunc(8)  = d3212m12p(xcos)*a12flag
      listdfunc(9)  = d3212m12p(xcos)*a12flag
      listdfunc(10) = d3212m12m(xcos)*a12flag
      listdfunc(11) = d3212m12m(xcos)*a12flag
      listdfunc(12) = d3212m32m(xcos)*a12flag
      listdfunc(13) = d3232p32p(xcos)*a32flag
      listdfunc(14) = d3232p12p(xcos)*a32flag
      listdfunc(15) = d3232p12p(xcos)*a32flag
      listdfunc(16) = d3232p12m(xcos)*a32flag
      listdfunc(17) = d3232p12m(xcos)*a32flag
      listdfunc(18) = d3232p32m(xcos)*a32flag
      listdfunc(19) = d3232m32p(xcos)*a32flag
      listdfunc(20) = d3232m12p(xcos)*a32flag
      listdfunc(21) = d3232m12p(xcos)*a32flag
      listdfunc(22) = d3232m12m(xcos)*a32flag
      listdfunc(23) = d3232m12m(xcos)*a32flag
      listdfunc(24) = d3232m32m(xcos)*a32flag
      ampsqr32pp = 0.d0
      do i1=1,6
        ampsqr32pp=ampsqr32pp+ampsqr(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos)
      enddo
      end function ampsqr32pp

!     -------------------------------------------------------------------------
!     Sum of amplitudes squared for helicity flip, spin-3/2
!     -------------------------------------------------------------------------
      real(8) function ampsqr32pm(egam,tman,xGamr,xMassr)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
      real(8) :: a0,ap,Anorm,b0,st,xBrn
      real(8) :: a12flag,a32flag
!     Internal Variables
      real(8) :: dxsnorm,dxs_spin32temp,dfunc,deltafunc,jr,xcos
      real(8), dimension (24) :: listdfunc,listdelfunc
      integer :: i1,angmom
!     External Variables
      real (8),external :: ampsqr
      real (8), external :: d3212p32p,d3212p12p,d3212p12m,d3212p32m
      real (8), external :: d3212m32p,d3212m12p,d3212m12m,d3212m32m
      real (8), external :: d3232p32p,d3232p12p,d3232p12m,d3232p32m
      real (8), external :: d3232m32p,d3232m12p,d3232m12m,d3232m32m
      jr = 3.d0/2.d0
      angmom = 0
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      a12flag    = choice(1)
      a32flag    = choice(2)
      listdelfunc(1)  = 0.d0
      listdelfunc(2)  = 1.d0
      listdelfunc(3)  = 0.d0
      listdelfunc(4)  = 0.d0
      listdelfunc(5)  = 0.d0
      listdelfunc(6)  = 0.d0
      listdelfunc(7)  = 0.d0
      listdelfunc(8)  = 0.d0
      listdelfunc(9)  = 0.d0
      listdelfunc(10) = 1.d0
      listdelfunc(11) = 0.d0
      listdelfunc(12) = 0.d0
      listdelfunc(13) = 1.d0
      listdelfunc(14) = 0.d0
      listdelfunc(15) = 0.d0
      listdelfunc(16) = 0.d0
      listdelfunc(17) = 0.d0
      listdelfunc(18) = 0.d0
      listdelfunc(19) = 0.d0
      listdelfunc(20) = 0.d0
      listdelfunc(21) = 0.d0
      listdelfunc(22) = 0.d0
      listdelfunc(23) = 0.d0
      listdelfunc(24) = 1.d0
      listdfunc(1)  = d3212p32p(xcos)*a12flag
      listdfunc(2)  = d3212p12p(xcos)*a12flag
      listdfunc(3)  = d3212p12p(xcos)*a12flag
      listdfunc(4)  = d3212p12m(xcos)*a12flag
      listdfunc(5)  = d3212p12m(xcos)*a12flag
      listdfunc(6)  = d3212p32m(xcos)*a12flag
      listdfunc(7)  = d3212m32p(xcos)*a12flag
      listdfunc(8)  = d3212m12p(xcos)*a12flag
      listdfunc(9)  = d3212m12p(xcos)*a12flag
      listdfunc(10) = d3212m12m(xcos)*a12flag
      listdfunc(11) = d3212m12m(xcos)*a12flag
      listdfunc(12) = d3212m32m(xcos)*a12flag
      listdfunc(13) = d3232p32p(xcos)*a32flag
      listdfunc(14) = d3232p12p(xcos)*a32flag
      listdfunc(15) = d3232p12p(xcos)*a32flag
      listdfunc(16) = d3232p12m(xcos)*a32flag
      listdfunc(17) = d3232p12m(xcos)*a32flag
      listdfunc(18) = d3232p32m(xcos)*a32flag
      listdfunc(19) = d3232m32p(xcos)*a32flag
      listdfunc(20) = d3232m12p(xcos)*a32flag
      listdfunc(21) = d3232m12p(xcos)*a32flag
      listdfunc(22) = d3232m12m(xcos)*a32flag
      listdfunc(23) = d3232m12m(xcos)*a32flag
      listdfunc(24) = d3232m32m(xcos)*a32flag
      ampsqr32pm = 0.d0
      do i1=13,18
        ampsqr32pm=ampsqr32pm+ampsqr(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos)
      enddo
      end function ampsqr32pm



!     -------------------------------------------------------------------------
!     Sum of amplitudes squared for helicity non-flip, spin-5/2
!     -------------------------------------------------------------------------
      real(8) function ampsqr52pp(egam,tman,xGamr,xMassr)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
      real(8) :: a0,ap,Anorm,b0,st,xBrn
      real(8) :: a12flag,a32flag
!     Internal Variables
      real(8) :: dxsnorm,dxs_spin52temp,dfunc,deltafunc,jr,xcos
      real(8), dimension (24) :: listdfunc,listdelfunc
      integer :: i1,angmom
!     External Variables
      real (8),external :: ampsqr
      real (8), external :: d5212p32p,d5212p12p,d5212p12m,d5212p32m
      real (8), external :: d5212m32p,d5212m12p,d5212m12m,d5212m32m
      real (8), external :: d5232p32p,d5232p12p,d5232p12m,d5232p32m
      real (8), external :: d5232m32p,d5232m12p,d5232m12m,d5232m32m
      jr = 5.d0/2.d0
      angmom = 1
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      a12flag    = choice(1)
      a32flag    = choice(2)
      listdelfunc(1)  = 0.d0
      listdelfunc(2)  = 1.d0
      listdelfunc(3)  = 0.d0
      listdelfunc(4)  = 0.d0
      listdelfunc(5)  = 0.d0
      listdelfunc(6)  = 0.d0
      listdelfunc(7)  = 0.d0
      listdelfunc(8)  = 0.d0
      listdelfunc(9)  = 0.d0
      listdelfunc(10) = 1.d0
      listdelfunc(11) = 0.d0
      listdelfunc(12) = 0.d0
      listdelfunc(13) = 1.d0
      listdelfunc(14) = 0.d0
      listdelfunc(15) = 0.d0
      listdelfunc(16) = 0.d0
      listdelfunc(17) = 0.d0
      listdelfunc(18) = 0.d0
      listdelfunc(19) = 0.d0
      listdelfunc(20) = 0.d0
      listdelfunc(21) = 0.d0
      listdelfunc(22) = 0.d0
      listdelfunc(23) = 0.d0
      listdelfunc(24) = 1.d0
      listdfunc(1)  = d5212p32p(xcos)*a12flag
      listdfunc(2)  = d5212p12p(xcos)*a12flag
      listdfunc(3)  = d5212p12p(xcos)*a12flag
      listdfunc(4)  = d5212p12m(xcos)*a12flag
      listdfunc(5)  = d5212p12m(xcos)*a12flag
      listdfunc(6)  = d5212p32m(xcos)*a12flag
      listdfunc(7)  = d5212m32p(xcos)*a12flag
      listdfunc(8)  = d5212m12p(xcos)*a12flag
      listdfunc(9)  = d5212m12p(xcos)*a12flag
      listdfunc(10) = d5212m12m(xcos)*a12flag
      listdfunc(11) = d5212m12m(xcos)*a12flag
      listdfunc(12) = d5212m32m(xcos)*a12flag
      listdfunc(13) = d5232p32p(xcos)*a32flag
      listdfunc(14) = d5232p12p(xcos)*a32flag
      listdfunc(15) = d5232p12p(xcos)*a32flag
      listdfunc(16) = d5232p12m(xcos)*a32flag
      listdfunc(17) = d5232p12m(xcos)*a32flag
      listdfunc(18) = d5232p32m(xcos)*a32flag
      listdfunc(19) = d5232m32p(xcos)*a32flag
      listdfunc(20) = d5232m12p(xcos)*a32flag
      listdfunc(21) = d5232m12p(xcos)*a32flag
      listdfunc(22) = d5232m12m(xcos)*a32flag
      listdfunc(23) = d5232m12m(xcos)*a32flag
      listdfunc(24) = d5232m32m(xcos)*a32flag
      ampsqr52pp = 0.d0
      do i1=1,6
        ampsqr52pp=ampsqr52pp+ampsqr(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos)
      enddo
      end function ampsqr52pp

!     -------------------------------------------------------------------------
!     Sum of amplitudes squared for helicity flip, spin-5/2
!     -------------------------------------------------------------------------
      real(8) function ampsqr52pm(egam,tman,xGamr,xMassr)
      use math_constants
      use physics_constants
      use hadron_properties
      use storage
      use size
      implicit none
!     Input Variables
      real(8) :: egam,tman,xGamr,xMassr
      real(8) :: a0,ap,Anorm,b0,st,xBrn
      real(8) :: a12flag,a32flag
!     Internal Variables
      real(8) :: dxsnorm,dxs_spin52temp,dfunc,deltafunc,jr,xcos
      real(8), dimension (24) :: listdfunc,listdelfunc
      integer :: i1,angmom
!     External Variables
      real (8),external :: ampsqr
      real (8), external :: d5212p32p,d5212p12p,d5212p12m,d5212p32m
      real (8), external :: d5212m32p,d5212m12p,d5212m12m,d5212m32m
      real (8), external :: d5232p32p,d5232p12p,d5232p12m,d5232p32m
      real (8), external :: d5232m32p,d5232m12p,d5232m12m,d5232m32m
      jr = 5.d0/2.d0
      angmom = 1
      if (tman.gt. 0.d0) then
         xcos = 1.d0
      else
         xcos =(2.d0*egam*tman - egam*xmjpsi**2 + 2.d0*egam**2*xmpro
     -	 + tman*xmpro - 
     -    xmjpsi**2*xmpro)/
     -  (egam*dsqrt(xmpro*(2.d0*egam + xmpro))*
     -    dsqrt((xmjpsi**4 + 4.d0*egam**2*xmpro**2
     -	 - 4.d0*xmjpsi**2*xmpro*(egam + xmpro))/
     -      (xmpro*(2.d0*egam + xmpro))))
      endif
      a0    = coef(1)
      ap    = coef(2)
      Anorm = coef(3)
      st    = coef(4)
      b0    = coef(5)
      xBrn  = coef(6)
      a12flag    = choice(1)
      a32flag    = choice(2)
      listdelfunc(1)  = 0.d0
      listdelfunc(2)  = 1.d0
      listdelfunc(3)  = 0.d0
      listdelfunc(4)  = 0.d0
      listdelfunc(5)  = 0.d0
      listdelfunc(6)  = 0.d0
      listdelfunc(7)  = 0.d0
      listdelfunc(8)  = 0.d0
      listdelfunc(9)  = 0.d0
      listdelfunc(10) = 1.d0
      listdelfunc(11) = 0.d0
      listdelfunc(12) = 0.d0
      listdelfunc(13) = 1.d0
      listdelfunc(14) = 0.d0
      listdelfunc(15) = 0.d0
      listdelfunc(16) = 0.d0
      listdelfunc(17) = 0.d0
      listdelfunc(18) = 0.d0
      listdelfunc(19) = 0.d0
      listdelfunc(20) = 0.d0
      listdelfunc(21) = 0.d0
      listdelfunc(22) = 0.d0
      listdelfunc(23) = 0.d0
      listdelfunc(24) = 1.d0
      listdfunc(1)  = d5212p32p(xcos)*a12flag
      listdfunc(2)  = d5212p12p(xcos)*a12flag
      listdfunc(3)  = d5212p12p(xcos)*a12flag
      listdfunc(4)  = d5212p12m(xcos)*a12flag
      listdfunc(5)  = d5212p12m(xcos)*a12flag
      listdfunc(6)  = d5212p32m(xcos)*a12flag
      listdfunc(7)  = d5212m32p(xcos)*a12flag
      listdfunc(8)  = d5212m12p(xcos)*a12flag
      listdfunc(9)  = d5212m12p(xcos)*a12flag
      listdfunc(10) = d5212m12m(xcos)*a12flag
      listdfunc(11) = d5212m12m(xcos)*a12flag
      listdfunc(12) = d5212m32m(xcos)*a12flag
      listdfunc(13) = d5232p32p(xcos)*a32flag
      listdfunc(14) = d5232p12p(xcos)*a32flag
      listdfunc(15) = d5232p12p(xcos)*a32flag
      listdfunc(16) = d5232p12m(xcos)*a32flag
      listdfunc(17) = d5232p12m(xcos)*a32flag
      listdfunc(18) = d5232p32m(xcos)*a32flag
      listdfunc(19) = d5232m32p(xcos)*a32flag
      listdfunc(20) = d5232m12p(xcos)*a32flag
      listdfunc(21) = d5232m12p(xcos)*a32flag
      listdfunc(22) = d5232m12m(xcos)*a32flag
      listdfunc(23) = d5232m12m(xcos)*a32flag
      listdfunc(24) = d5232m32m(xcos)*a32flag
      ampsqr52pm = 0.d0
      do i1=13,18
        ampsqr52pm=ampsqr52pm+ampsqr(egam,tman,xGamr,xMassr,
     -  listdfunc(i1),listdelfunc(i1),
     -	angmom,jr,xcos)
      enddo
      end function ampsqr52pm

!————————————————————————————————————
!     d-functions
!     -------------------------------------------------------------------------
      real(8) function d3212p32p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3212p32p = (dsqrt(1.5d0)*dsqrt(1.d0 - xcos)*(1.d0 + xcos))/2.d0
      end function d3212p32p      

      real(8) function d3212p12p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3212p12p = (dsqrt(1.d0 + xcos)*(-1.d0 + 3.d0*xcos))
     -	/(2.d0*dsqrt(2.d0))
      end function d3212p12p  

      real(8) function d3212p12m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3212p12m = -(dsqrt(2.d0 - 2.d0*xcos)*(1.d0 + 3.d0*xcos))/4.d0
      end function d3212p12m      

      real(8) function d3212p32m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3212p32m = -(dsqrt(1.5d0)*(-1.d0 + xcos)
     -	*dsqrt(1.d0 + xcos))/2.d0
      end function d3212p32m    

      real(8) function d3212m32p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3212m32p = -(dsqrt(1.5d0)*(-1.d0 + xcos)
     -	*dsqrt(1.d0 + xcos))/2.d0
      end function d3212m32p      

      real(8) function d3212m12p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3212m12p = (dsqrt(2.d0 - 2.d0*xcos)*(1.d0 + 3.d0*xcos))/4.d0
      end function d3212m12p  

      real(8) function d3212m12m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3212m12m = (dsqrt(1.d0 + xcos)*(-1.d0 + 3.d0*xcos))
     -	/(2.d0*dsqrt(2.d0))
      end function d3212m12m      

      real(8) function d3212m32m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3212m32m = -(dsqrt(1.5d0)*dsqrt(1.d0 - xcos)*(1.d0 + xcos))/2.d0
      end function d3212m32m    

      real(8) function d3232p32p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3232p32p = (1.d0 + xcos)**1.5d0/(2.d0*dsqrt(2.d0))
      end function d3232p32p      

      real(8) function d3232p12p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3232p12p = -(dsqrt(1.5d0)*dsqrt(1.d0 - xcos)*(1.d0 + xcos))/2.d0
      end function d3232p12p  

      real(8) function d3232p12m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3232p12m = -(dsqrt(1.5d0)*(-1.d0 + xcos)
     -	*dsqrt(1.d0 + xcos))/2.d0
      end function d3232p12m      

      real(8) function d3232p32m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3232p32m = -(1.d0 - xcos)**1.5d0/(2.d0*dsqrt(2.d0))
      end function d3232p32m    

      real(8) function d3232m32p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3232m32p = (1.d0 - xcos)**1.5/(2.d0*dsqrt(2.d0))
      end function d3232m32p      

      real(8) function d3232m12p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3232m12p = -(dsqrt(1.5d0)*(-1.d0 + xcos)
     -	*dsqrt(1.d0 + xcos))/2.d0
      end function d3232m12p  

      real(8) function d3232m12m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3232m12m = (dsqrt(1.5d0)*dsqrt(1.d0 - xcos)*(1.d0 + xcos))/2.d0
      end function d3232m12m      

      real(8) function d3232m32m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d3232m32m = (1.d0 + xcos)**1.5/(2.d0*dsqrt(2.d0))
      end function d3232m32m    

      real(8) function d5212p32p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5212p32p = (dsqrt(1.d0 - xcos)*(1.d0 + xcos)
     -	*(-1.d0 + 5.d0*xcos))/4.d0
      end function d5212p32p      

      real(8) function d5212p12p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5212p12p = (dsqrt(1.d0 + xcos)*(-1.d0 - 2.d0*xcos
     -	 + 5.d0*xcos**2))/(2.d0*dsqrt(2.d0))
      end function d5212p12p  

      real(8) function d5212p12m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5212p12m = (dsqrt(2.d0 - 2.d0*xcos)*(1.d0 - 2.d0*xcos
     -	 - 5.d0*xcos**2))/4.d0
      end function d5212p12m      

      real(8) function d5212p32m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5212p32m = -((-1.d0 + xcos)*dsqrt(1.d0 + xcos)
     -	*(1.d0 + 5.d0*xcos))/4.d0
      end function d5212p32m    

      real(8) function d5212m32p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5212m32p = -((-1.d0 + xcos)*dsqrt(1.d0 + xcos)
     -	*(1.d0 + 5.d0*xcos))/4.d0
      end function d5212m32p      

      real(8) function d5212m12p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5212m12p = (dsqrt(2.d0 - 2.d0*xcos)
     -	*(-1.d0 + 2.d0*xcos + 5.d0*xcos**2))/4.d0
      end function d5212m12p  

      real(8) function d5212m12m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5212m12m = (dsqrt(1.d0 + xcos)*(-1.d0 - 2.d0*xcos
     -	 + 5.d0*xcos**2))/(2.d0*dsqrt(2.d0))
      end function d5212m12m      

      real(8) function d5212m32m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5212m32m = -(dsqrt(1.d0 - xcos)*(1.d0 + xcos)
     -	*(-1.d0 + 5.d0*xcos))/4.d0
      end function d5212m32m    

      real(8) function d5232p32p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5232p32p = ((1.d0 + xcos)**1.5*(-3.d0 + 5.d0*xcos))
     -	/(4.d0*dsqrt(2.d0))
      end function d5232p32p      

      real(8) function d5232p12p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5232p12p = -(dsqrt(1.d0 - xcos)*(1.d0 + xcos)*
     -	(-1.d0 + 5.d0*xcos))/4.d0
      end function d5232p12p  

      real(8) function d5232p12m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5232p12m = -((-1.d0 + xcos)*dsqrt(1.d0 + xcos)
     -	*(1.d0 + 5.d0*xcos))/4.d0
      end function d5232p12m      

      real(8) function d5232p32m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5232p32m = -((1.d0 - xcos)**1.5*(3.d0 + 5.d0*xcos))
     -	/(4.d0*dsqrt(2.d0))
      end function d5232p32m    

      real(8) function d5232m32p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5232m32p =  ((1.d0 - xcos)**1.5*(3.d0 + 5.d0*xcos))
     -	/(4.d0*dsqrt(2.d0))
      end function d5232m32p      

      real(8) function d5232m12p(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5232m12p = -((-1.d0 + xcos)*dsqrt(1.d0 + xcos)
     -	*(1.d0 + 5.d0*xcos))/4.d0
      end function d5232m12p  

      real(8) function d5232m12m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5232m12m = (dsqrt(1.d0 - xcos)*(1.d0 + xcos)
     -	*(-1.d0 + 5.d0*xcos))/4.d0
      end function d5232m12m      

      real(8) function d5232m32m(xcos)
      implicit none
!     Input Variables
      real(8) :: xcos
      d5232m32m = ((1.d0 + xcos)**1.5*(-3.d0 + 5.d0*xcos))
     -	/(4.d0*dsqrt(2.d0))
      end function d5232m32m



