detbal.f

c $Id: detbal.f,v 1.12 1999/01/18 09:57:00 ernst Exp $
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      subroutine detbal(sqrts,ityp1,ityp2,iso31,iso32,
     &                  em1,em2,itnew1,itnew2,dbfact)
c
c     Revision : 1.0
c
cinput sqrts   : sqrt(s)
cinput ityp1   : ityp of incoming particle 1
cinput ityp2   : ityp of incoming particle 2
cinput iso31   : 2*I3 of incoming particle 1
cinput iso32   : 2*I3 of incoming particle 2
cinput em1     : mass of incoming particle 1
cinput em2     : mass of incoming particle 2
cinput itnew1  : ityp of outgoing particle 1
cinput itnew2  : ityp of outgoing particle 2
c
coutput dbfact     : correction factor for cross section
c
c     This subroutine calculates a correction factor for the
c     partial crosssection based on the principle of detailed balance.
C
c     For {\tt CTOption(3)=0} a modified detailed balance (default) is used
c     which takes finite resonance widths into account. For
c     {\tt CTOption(3)=1} the old standard detailed balance relation is used.
c
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c
      implicit none

      include 'coms.f'
      include 'options.f'
      include 'comres.f'
      include 'newpart.f'
c

      real*8 sqrts, dbfact
      integer iso31, iso32
      real*8  em1, em2, clebweight
      integer ityp1, ityp2, itnew1, itnew2
c     local vars for integration of Breit-Wigner:
      real*8 mepsilon
      real*8 oq, q, minwid, factor
      integer inres, outres,idum1,idum2,idum3,idum4
c     called functions
      real*8 pcms, dbweight, massit
      real*8 pmean, widit, dgcgkfct
      integer isoit
      real*8 detbalin
      external detbalin

c
c     0.1 MeV shift for integrator-maxvalue
      parameter(mepsilon=0.0001)
c     minimal width for "unstable" particle
      parameter( minwid=1.d-4 )

ctp060202 to avoid warnings with gfortran compilation
      logical ctp060202
      ctp060202=.false.
      if(ctp060202)write(*,*)em1,em2
ctp060202 end

      idum1=0
      idum2=0
      idum3=0
      idum4=0
c
c     fix itypes, iso3 and  phase-space for outgoing particles
c
c     a) set up call to isocgk and getmass: determine outgoing isospins
c        and masses


c
c clebweight: actually areduction of given isospin_summed cross_section
c             to actual incoming channel - here it is used to probe
c             wether the process in question is isospin allowed or not
c
      clebweight=dbweight(isoit(ityp1),iso31,isoit(ityp2),iso32,
     &     isoit(itnew1),isoit(itnew2))
      if(clebweight.lt.0.00001) then
         dbfact=0.d0
         return
      endif


c
c     b) determine momenta
c
      pnnout=pcms(sqrts,massit(itnew1),massit(itnew2))
c
c     d) now calculate correction factor
c

c
c     call to dgcgkfct which calculates degeneracy factors and clebsches
c
      factor=dgcgkfct(ityp1,ityp2,iso31,iso32,itnew1,itnew2)


      if(CTOption(3).eq.0) then

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c modified detailed balance
c
c reference: Danielewicz and Bertsch: Nuclear Physics A533(1991) 712.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


c
c     count resonances in the incoming channel:
c
         inres=0
         if(widit(ityp1).gt.minwid) inres=inres+1
         if(widit(ityp2).gt.minwid) inres=inres+1

c
c     count resonances in the outgoing channel:
         outres=0
         if(widit(itnew1).gt.minwid) outres=outres+1
         if(widit(itnew2).gt.minwid) outres=outres+1

         if(inres.eq.0) then
c     in this case detbal without resonances, original prescription
            dbfact=factor*pnnout**2/(pnn**2)
c
cccccccccccccccccccccccccccccc
         elseif(inres.eq.1) then
c modified det-bal for one resonance
c

c     now generate the correction factor
         dbfact=factor*pnnout**2/
     &          pmean(sqrts,ityp1,iso31,ityp2,iso32,
     &                idum1,idum2,idum3,idum4,2)

cccccccccccccccccccccccccc
      else
c     modified det-bal for two resonances
c     reference: S.A. Bass, private calculation
c
ccccccccccccccccccccccccccccccccc
         oq=0D0
         if(outres.gt.0) then

c here we have B* B* to B* N
            oq=pmean(sqrts,itnew1,-99,itnew2,-99,
     &               idum1,idum2,idum3,idum4,2)

         endif
ccccccccccccccccccccccccccccccccc

         q=pmean(sqrts,ityp1,iso31,ityp2,iso32,
     &           idum1,idum2,idum3,idum4,2)
c
c     now generate the correction factor
            if(outres.eq.0) then
               dbfact=factor*pnnout**2/(max(1.d-12,q))
            else
               dbfact=factor*oq/(max(1.d-12,q))
            endif
         endif

         return
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c original detailed balance
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      elseif(CTOption(3).eq.1) then
         dbfact=factor*pnnout**2/(pnn**2)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c error processing
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      else
         write(6,*)'undefined detailed balance mode in DETBAL'
         dbfact=1.
      endif
c
      return
      end


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      real*8 function ppiso(pid,ityp1,iso31,ityp2,iso32,itnew1,itnew2)
c
c     Revision : 1.0
c
cinput pid     : ID of process
cinput ityp1   : ityp of incoming particle 1
cinput ityp2   : ityp of incoming particle 2
cinput iso31   : 2*I3 of incoming particle 1
cinput iso32   : 2*I3 of incoming particle 2
cinput itnew1  : ityp of outgoing particle 1
cinput itnew2  : ityp of outgoing particle 2
c
coutput nniso     : isospin-factor for the reaction $p p \to B B$
c
c     This subroutine calculates the isospin-factor for resonance
c     excitation in inelastic proton proton collisions.
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c
      implicit none

      real*8 dbweight,factor,dgcgkfct
      integer isoit,pid,itag,iz1,iz2,i1,i2,im,jm,ityp1,ityp2
      integer iso31,iso32,itnew1,itnew2,itmp1,itmp2

      include 'comres.f'
      include 'newpart.f'


      i1=ityp1
      i2=ityp2
      iz1=iso31
      iz2=iso32
      im=itnew1
      jm=itnew2

      if(pid.gt.0) then
         ppiso=dbweight(i1,iz1,i2,iz2,isoit(im),isoit(jm))/
     /        dbweight(1,1,1,1,isoit(im),isoit(jm))
      else

         factor=dgcgkfct(i1,i2,iz1,iz2,nucleon,nucleon)
         if(factor.le.1.d-8) then
            ppiso=0.d0
            return
         endif

         nexit=2
         itot(1)=isoit(nucleon)
         itot(2)=isoit(nucleon)
         call isocgk4(isoit(i1),iz1,isoit(i2),iz2,itot,i3new,itag)
         itmp1=i3new(1)
         itmp2=i3new(2)
         ppiso=dbweight(nucleon,itmp1,nucleon,itmp2,
     &           isoit(im),isoit(jm))/
     /           dbweight(1,1,1,1,isoit(im),isoit(jm))
      endif

      return

      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      real*8 function dgcgkfct(ityp1,ityp2,iso31,iso32,itnew1,itnew2)
c
c     Revision : 1.0
c
cinput ityp1   : ityp of incoming particle 1
cinput ityp2   : ityp of incoming particle 2
cinput iso31   : 2*I3 of incoming particle 1
cinput iso32   : 2*I3 of incoming particle 2
cinput itnew1  : ityp of outgoing particle 1
cinput itnew2  : ityp of outgoing particle 2
c
coutput dgcgkfct     : product of degeneracy and cgk factor for detailed bal.
c
c     This subroutine calculates the product of the spin and isospin
c     degeneracy factors and the a isospin correction factor
c     (isospin dependence of cross section) for the detailed balance
c     cross section.
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c
      implicit none

      include 'comres.f'
c

      integer iso31,iso32
      real*8  clebweight
      integer ityp1, ityp2,stot(4),itnew1,itnew2
c     components of degeneracy factor
      integer gin1,gin2,gout1,gout2
      real*8 dgfact
c     called functions
      real*8 dbweight
      integer jit,isoit
c
c     a) set up call to isocgk and getmass: determine outgoing isospins
c        and masses

c
c clebweight: reduction of given isospin_summed cross_section to actual
c             incoming channel
c
c 1) reduction:
      clebweight=dbweight(isoit(ityp1),iso31,isoit(ityp2),iso32,
     &     isoit(itnew1),isoit(itnew2))
      if(clebweight.lt.0.00001) then
         dgcgkfct=0.d0
         return
      endif


c     c) calculate degeneracy factors
c        reference: S. Bass, GSI-Report 93-13 p. 25 and references therein
c
c     get spins: in-channel stot(1 and 2), out-channel stot(3 and 4)
      stot(1)=jit(ityp1)
      stot(2)=jit(ityp2)
      stot(3)=jit(itnew1)
      stot(4)=jit(itnew2)
c

      gout1=(stot(3)+1)
      gout2=(stot(4)+1)
      gin1=(stot(1)+1)
      gin2=(stot(2)+1)
c
c
c     the degeneracy factor is
      dgfact=dble(gout1*gout2)/dble(gin1*gin2)
c
c     d) now calculate correction factor
c
c
      dgcgkfct=dgfact*clebweight
c
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real*8 function pmean(sqrts,itp1,iz1,itp2,iz2,
     &                      itp3,iz3,itp4,iz4,ipwr)
c
c     Revision : 1.0
c
cinput sqrts  : $\sqrt{s}$
cinput itp1   : ityp of particle 1
cinput iz1    : $2 \cdot I_3$ of particle 1
cinput itp2   : ityp of particle 2
cinput iz2    : $2 \cdot I_3$ of particle 2
cinput itp3   : ityp of particle 3
cinput iz3    : $2 \cdot I_3$ of particle 3
cinput itp4   : ityp of particle 4
cinput iz4    : $2 \cdot I_3$ of particle 4
cinput ipwr   : power of $p_{mean}$ to integrate
c
c     This function returns the value of the following integral:
c     \begin{displaymath}
c     \int\limits_{m_1= {\tt mmin}}^{\tt mmax}
c      p_{CMS}^{\tt ipwr}(\sqrt{s},m_1,m_2) A_1(m_1) A_2(m_2) \; dm_1 dm_2
c     \end{displaymath}
c     with $A_r(m)$ being the spectral function of the resonance:
c     \begin{displaymath}
c      A(m) = \displaystyle\frac{\Gamma(m)/2}{(m-m_0)^2+\Gamma(m)^2/4}
c     \end{displaymath}
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none

c      include "comres.f"

      real*8 sqrts,minwid,q1,q2
      real*8 mmin1,mfix,mmin2,maxs1,maxs2,mepsilon,diverg1
      real*8 smass
      integer itp1,ipwr,iz1,itp2,iz2,itp3,iz3,itp4,iz4
      integer inres,izz1,ind1,ind2

cfunctions
      real*8 detbalin,widit,massit,detbalin2,mminit,pcms
c      integer
      external detbalin,detbalin2

c     minimal width for "unstable" particle
      parameter( minwid=1.d-3 )
c     0.1 MeV shift for integrator-maxvalue
      parameter(mepsilon=0.0001)
ctp060202 to avoid warnings with gfortran compilation
      logical ctp060202
      ctp060202=.false.
      if(ctp060202)write(*,*)iz3,iz4
ctp060202 end


      if(sqrts.le.mminit(itp1)+mminit(itp2)) then
         pmean=0.d0
         return
      endif

c     count broad particles and store number in inres
c     NOTE: only particles 1 and 2 may be broad!!!!
c
      inres=0
      if(widit(itp1).gt.minwid) inres=inres+1
      if(widit(itp2).gt.minwid) inres=inres+1

      if(inres.eq.0) then
c     in this case the Breit-Wigner distributions are Delta functions,
c     no integrations necessary
         smass=massit(itp2)
         if(itp3.ne.0) smass=smass+massit(itp3)
         if(itp4.ne.0) smass=smass+massit(itp4)

         pmean=pcms(sqrts,massit(itp1),smass)**ipwr

         return
c
cccccccccccccccccccccccccccccc
      elseif(inres.eq.1) then
c modified det-bal for one resonance
c
c first determine which particle is the resonance and store ityp in ind1
         if(widit(itp1).gt.minwid) then
            ind2=itp2
            ind1=itp1
            izz1=iz1
         else
            ind2=itp1
            ind1=itp2
            izz1=iz2
         endif

c     now set integration boundaries
         mmin1=mminit(ind1)
         mfix=mminit(ind2)

         if(itp3.ne.0) mfix=mfix+massit(itp3)
         if(itp4.ne.0) mfix=mfix+massit(itp4)

         maxs1=sqrts-mfix-mepsilon
c     the integration might be divided by the pole of the Breit-Wigner
c     then two integrations with diverg1 as upper or lower boundary
c     respectively are necessary
         diverg1=massit(ind1)
c
c     now perform integration in function detbalin
c     integrate f(m1)=pcms(sqrts,m1,m2)**ipwr*fbwnorm(m1,ityp1)
         q1=0.d0
         q2=0.d0
         if(mmin1.le.diverg1) then
            if(maxs1.gt.diverg1) then
               call qsimp(detbalin,mmin1,diverg1,
     &              ind1,izz1,mfix,sqrts,ipwr,q1,-1)
               call qsimp(detbalin,diverg1,maxs1,
     &              ind1,izz1,mfix,sqrts,ipwr,q2,1)
            else
               call qsimp(detbalin,mmin1,maxs1,
     &              ind1,izz1,mfix,sqrts,ipwr,q1,-1)
            endif
         else
               call qsimp(detbalin,mmin1,maxs1,
     &              ind1,izz1,mfix,sqrts,ipwr,q2,1)
         endif

         pmean=(q1+q2)

         return
c
cccccccccccccccccccccccccc
      else
c 2 resonances to integrate over
c
c
         if(itp3.ne.0) then
            write(6,*) 'ERROR in pmean: only one broad particle allowed'
            write(6,*) '                in case of 3 or 4 body decays!!'
            stop
         endif


c     outer integration:
c     set integration boundaries
         mmin1=mminit(itp1)
         mmin2=mminit(itp2)
         maxs2=sqrts-mmin1
         diverg1=massit(itp2)
         q1=0.d0
         q2=0.d0
         if(mmin2.le.diverg1) then
            if(maxs2.gt.diverg1) then
               call qsimp2(detbalin2,mmin2,diverg1,
     &              itp2,iz2,mmin1,itp1,iz1,ipwr,sqrts,q1,-1)
               call qsimp2(detbalin2,diverg1,maxs2,
     &              itp2,iz2,mmin1,itp1,iz1,ipwr,sqrts,q1,1)
            else
               call qsimp2(detbalin2,mmin2,maxs2,
     &              itp2,iz2,mmin1,itp1,iz1,ipwr,sqrts,q1,-1)
            endif
         else
            call qsimp2(detbalin2,mmin2,maxs2,
     &           itp2,iz2,mmin1,itp1,iz1,ipwr,sqrts,q1,1)
         endif

         pmean=(q1+q2)

         return
      endif


      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      real*8 function detbalin(m1,ityp1,iz1,m2,sqrts,ipwr)
c
c     Revision : 1.0
c
cinput m1     : mass of resonance (integration variable)
cinput ityp1  : ityp of Delta/N* resonance for fbwnorm()
cinput iz1    : $2\cdot I_3$ of resonance
cinput m2     : second mass for call to pcms
cinput sqrts  : sqrt(s)
cinput ipwr   : power for $p_mean$
c
c     This function is an integrand for the modified detailed balance:
c     \begin{displaymath}
c     detbalin(m1)=\, p_{CMS}^{\tt ipwr}(\sqrt{s},m_1,m_2) A_1(m_1)
c     \end{displaymath}
c     with $A_r(M)$ being the spectral function of the resonance:
c     \begin{displaymath}
c      A(m) = \displaystyle\frac{\Gamma(m)/2}{(m-m_0)^2+\Gamma(m)^2/4}
c     \end{displaymath}
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     arguments
      real*8 m1,m2,sqrts
      integer ityp1,iz1,ipwr
c     called functions
      real*8 fbwnorm,pcms

      detbalin=pcms(sqrts,m1,m2)**ipwr*fbwnorm(m1,ityp1,iz1)

      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      real*8 function detbalin2(m2,ityp2,iz2,min1,
     &                          ityp1,iz1,ipwr,sqrts)
c
c     Revision : 1.0
c
c     This function represents the integrand
c     \begin{displaymath}
c     detbalin2\,=\,A_2(m2)\;
c     \int \limits_{m_N + m_{\pi}}^{\sqrt{s} - m_2}
c     p_{rel}(\sqrt{s},m_1,m_2) A_1(m_1)  \, d m_1
c     \end{displaymath}
c     for the modified detailed balance with two resonances in the
c     incoming channel.
c     $A_r(M)$ is the spectral function of the resonance (see {\tt detbalin}).
c
c
cinput m2     : mass of resonance2 (outer integration in detbal)
cinput ityp1  : ityp of resonance1
cinput ityp2  : ityp of resonance2
cinput min1   : lower boundary for integration via {\tt qsimp}
ccinput max1   : upper boundary for integration via {\tt qsimp}
cinput iz1    : $2\cdot I_3$ of resonance 1
cinput iz2    : $2\cdot I_3$ of resonance 2
cinput ipwr   : power for $p_mean$
cinput sqrts  : $\sqrt{s}$
c
c     output:
c             detbalin2  : value of function
c
c     function and subroutine calls:
c                      detbalin (referenced as external)
c                      fbwnorm
c                      qsimp
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer ityp1,ityp2,iz1,iz2,ipwr
      real*8 m2,min1,max1,sqrts,q1,q2,diverg1,fbwnorm,mepsilon
      real*8 massit
      real*8 detbalin
      external detbalin
c     0.1 MeV shift for integrator-maxvalue
      parameter(mepsilon=0.0001)

      max1=sqrts-m2-mepsilon

      diverg1=massit(ityp1)
      q1=0.d0
      q2=0.d0
      if(min1.le.diverg1) then
         if(max1.gt.diverg1) then
            call qsimp(detbalin,min1,diverg1,
     &           ityp1,iz1,m2,sqrts,ipwr,q1,-1)
            call qsimp(detbalin,diverg1,max1,
     &           ityp1,iz1,m2,sqrts,ipwr,q2,1)
         else
            call qsimp(detbalin,min1,max1,
     &           ityp1,iz1,m2,sqrts,ipwr,q1,-1)
         endif
      else
            call qsimp(detbalin,min1,max1,
     &           ityp1,iz1,m2,sqrts,ipwr,q2,1)
      endif

      detbalin2=fbwnorm(m2,ityp2,iz2)*(q1+q2)

      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real*8 function fbwnorm(m,ires,iz1)
c
c     Revision : 1.0
c
cinput m   : mass of resonance
cinput ires: ityp of resonance
cinput iz1 : $2\cdot I_3$ of resonance
c
c     {\tt fbwnorm} returns a Breit-Wigner distribution (non-relativistic)
c     which is normalized to 1 in the limit of mass-independent
c     decay widths. However this function uses mass-dependent decay
c     widths when available. The function only uses widths down to
c     a lower boundary of 1 MeV, smaller widths are automatically set
c     to 1 MeV. For {\tt iz=-99} fixed widths are used instead of
c     a call to {\tt fwidth}. You should use {\tt fbrwig} for standard purpose
c     since in case of mass dependent widths fbwnorm() is not very well
c     defined for widths smaller than 1 MeV.
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      implicit none

      real*8 m,gam2,mres,fwidth,massit,widit,gam,minwid
      integer ires,ires1,iz1
      include 'comres.f'
      include 'coms.f'
      include 'comwid.f'

c     minimal width for "unstable" particle
      parameter( minwid=1.d-3 )

      ires1 = ires
      mres = massit(ires1)
      if(iz1.eq.-99.or.wtabflg.eq.0)then
        gam = widit(ires1)
      else
        gam = fwidth(ires1,iz1,m)
      end if
c     cutoff for small widths
      gam=max(gam,minwid)
      gam2=gam**2
      fbwnorm = 0.5*gam/(pi*((m-mres)**2+gam2/4.0))!*norm
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c  here start the numerical receipies routines for numerical integration
c  no more physics beyond this point in the file!!!
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c(c) numerical receipies, adapted for f(x,idum,dum,dum)
      SUBROUTINE qsimp(func,a,b,idum1,idum2,dum2,dum3,idum3,s,flag)
c
c     Simpson integration via Numerical Receipies.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      include 'options.f'

      INTEGER JMAX,j,idum1,idum2,idum3,flag
      REAL*8 a,b,func,s,EPS
      REAL*8 os,ost,st,dum2,dum3
      EXTERNAL func
      PARAMETER (JMAX=20)

      if(b-a.le.1.d-4) then
         s=0.d0
         return
      endif

      EPS = 6.d-2
      if (CTOption(35).eq.1) EPS=6.d-3

      ost=-1.d30
      os= -1.d30
      do 11 j=1,JMAX
         if(flag.eq.-1) then
            call midsqu1(func,a,b,idum1,idum2,dum2,dum3,idum3,st,j)
         elseif(flag.eq.1) then
            call midsql1(func,a,b,idum1,idum2,dum2,dum3,idum3,st,j)
         endif
        s=(9.*st-ost)/8.
        if (abs(s-os).le.EPS*abs(os)) return
        os=s
        ost=st
11    continue
      write(6,*)  'too many steps in qsimp, increase JMAX!'

      return
      END

C  (C) Copr. 1986-92 Numerical Recipes Software.



      SUBROUTINE midsqu1(funk,aa,bb,idum1,idum2,dum2,dum3,idum3,s,n)
c     modified midpoint rule; allows singuarity at upper limit
      implicit none
      integer idum1,idum2,idum3
      real*8 dum2,dum3
      INTEGER n
      REAL*8 aa,bb,s,funk
      EXTERNAL funk
      INTEGER it,j
      REAL*8 ddel,del,sum,tnm,x,func,a,b,xx

      func(x)=2.*x*funk(bb-x**2,idum1,idum2,dum2,dum3,idum3)

      b=sqrt(bb-aa)
      a=0.d0
      if (n.eq.1) then
      xx=0.5d0*(a+b)

      s=(b-a)*func(0.5d0*(a+b))

      else
        it=3**(n-2)
        tnm=it
        del=(b-a)/(3.*tnm)
        ddel=del+del
        x=a+0.5*del
        sum=0.
        do 11 j=1,it
          sum=sum+func(x)
          x=x+ddel
          sum=sum+func(x)
          x=x+del
11      continue
        s=(s+(b-a)*sum/tnm)/3.
      endif
      return
      END

      SUBROUTINE midsql1(funk,aa,bb,idum1,idum2,dum2,dum3,idum3,s,n)
c     modified midpoint rule; allows singularity at lower limit
      implicit none
      integer idum1,idum2,idum3
      real*8 dum2,dum3
      INTEGER n
      REAL*8 aa,bb,s,funk
      EXTERNAL funk
      INTEGER it,j
      REAL*8 ddel,del,sum,tnm,x,func,a,b
      func(x)=2.*x*funk(aa+x**2,idum1,idum2,dum2,dum3,idum3)
      b=sqrt(bb-aa)
      a=0.
      if (n.eq.1) then
        s=(b-a)*func(0.5*(a+b))
      else
        it=3**(n-2)
        tnm=it
        del=(b-a)/(3.*tnm)
        ddel=del+del
        x=a+0.5*del
        sum=0.
        do 11 j=1,it
          sum=sum+func(x)
          x=x+ddel
          sum=sum+func(x)
          x=x+del
11      continue
        s=(s+(b-a)*sum/tnm)/3.
      endif
      return
      END


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c(c) numerical receipies, adapted for f(x,idum,dum,dum)
      SUBROUTINE qsimp2(func,a,b,idum1,idum2,dum1,idum3,idum4,
     &                  idum5,dum2,s,flag)
c
c     Simpson integration via Numerical Receipies.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none

      include 'options.f'

      INTEGER JMAX,j,idum1,idum2,idum3,idum4,idum5,flag
      REAL*8 a,b,s,EPS
      REAL*8 os,ost,st,dum1,dum2
      REAL*8 func
      PARAMETER (JMAX=10)
      external func

      if(b-a.le.1.d-4) then
         s=0.d0
         return
      endif

      EPS = 6.d-2
      if (CTOption(35).eq.1) EPS=6.d-3

      ost=-1.d30
      os= -1.d30
      do 11 j=1,JMAX
         if(flag.eq.-1) then
            call midsqu2(func,a,b,idum1,idum2,dum1,idum3,idum4,
     &                  idum5,dum2,st,j)
         elseif(flag.eq.1) then
            call midsql2(func,a,b,idum1,idum2,dum1,idum3,idum4,
     &                  idum5,dum2,st,j)
         endif
        s=(9.*st-ost)/8.

        if (abs(s-os).le.EPS*abs(os)) return
        os=s
        ost=st
11    continue

      write(6,*)  'too many steps in qsimp2, increase JMAX!'

      return
      END


      SUBROUTINE midsqu2(funk,aa,bb,idum1,idum2,dum1,idum3,idum4,
     &                  idum5,dum2,s,n)
c     modified midpoint rule; allows singuarity at upper limit
      implicit none
      integer idum1,idum2,idum3,idum4,idum5
      real*8 dum1,dum2
      INTEGER n
      REAL*8 aa,bb,s,funk
      EXTERNAL funk
      INTEGER it,j
      REAL*8 ddel,del,sum,tnm,x,func,a,b
      func(x)=2.*x*funk(bb-x**2,
     &                  idum1,idum2,dum1,idum3,idum4,idum5,dum2)
      b=sqrt(bb-aa)
      a=0.
      if (n.eq.1) then
        s=(b-a)*func(0.5*(a+b))
      else
        it=3**(n-2)
        tnm=it
        del=(b-a)/(3.*tnm)
        ddel=del+del
        x=a+0.5*del
        sum=0.
        do 11 j=1,it
          sum=sum+func(x)
          x=x+ddel
          sum=sum+func(x)
          x=x+del
11      continue
        s=(s+(b-a)*sum/tnm)/3.
      endif
      return
      END

      SUBROUTINE midsql2(funk,aa,bb,idum1,idum2,dum1,idum3,idum4,
     &                  idum5,dum2,s,n)
c     modified midpoint rule; allows singularity at lower limit
      implicit none
      integer idum1,idum2,idum3,idum4,idum5
      real*8 dum1,dum2
      INTEGER n
      REAL*8 aa,bb,s,funk
      EXTERNAL funk
      INTEGER it,j
      REAL*8 ddel,del,sum,tnm,x,func,a,b
      func(x)=2.*x*funk(aa+x**2,
     &                  idum1,idum2,dum1,idum3,idum4,idum5,dum2)
      b=sqrt(bb-aa)
      a=0.
      if (n.eq.1) then
        s=(b-a)*func(0.5*(a+b))
      else
        it=3**(n-2)
        tnm=it
        del=(b-a)/(3.*tnm)
        ddel=del+del
        x=a+0.5*del
        sum=0.
        do 11 j=1,it
          sum=sum+func(x)
          x=x+ddel
          sum=sum+func(x)
          x=x+del
11      continue
        s=(s+(b-a)*sum/tnm)/3.
      endif
      return
      END