!/ module sgvbt1
!  version 1.4    25 november 1991  for transp
!  from version 1.2 in sigmax.for (version 1.3 is in survey)
!
!.              %%%%%%%%%%%%%%%%%%%%%%%%%
!.              %                       %
!.              %       sgvbt1          %
!.              %                       %
!.              %%%%%%%%%%%%%%%%%%%%%%%%%
!
function sgvbt1(peb,pab,pti,pai,kzfast,katom,ksig)
  !  this function uses a multi interval gauss-legendre algorithm to
  !  evaluate the one dimensional integral for beam-target fusion reactivity.
  !  input index:
  !  peb (ev) the lab energy of the fast ion.
  !  pab (amu) the mass of the fast ion.
  !  pti (ev) the temperature of the target ions.
  !  pai (amu) the mass of the target ions.
  !  output index:
  !  sgvbt1  (cm**3/sec) the beam-target fusion reactivity
  !
  !============
  use iso_c_binding, only: fp => c_double
  use physconst_mod, only : two, kb_ev, amucgs=>amu_g !amucgs is one atomic mass unit in grams
  implicit none
  integer :: kzfast,katom,ksig,l,kgh,kgl,kdim,igh,igl,int,jlgl,jgl
  !============
  real(fp) :: peb,pab,pti,pai,sgvbt1,zrglx,zrglw,zrsig
  real(fp), parameter :: evs = kb_ev*1.0e+7_fp  !evs is one ev in ergs.
  real(fp), parameter :: evsinv = 1.0e-7_fp/kb_ev !evsinv is one erg in ev.
  real(fp), parameter :: zconv = amucgs*evsinv/two !zconv is 0.5*amucgs*evsinv; 
                                      !amlab*zconv*pvrel2 is the lab energy in ev.
  real(fp) :: zvbeam,aioniz,zsgion,zvlowr,zvuppr,zdvint
  real(fp) :: zsum,zvcntr,zvrel,zebrel,zsig,zexp,za2exp
  !============
  dimension zrglx(7,7),zrglw(7,7),zrsig(5)
  !  gauss-legendre evaluation points (see abramowitz and stegun).
  data zrglx(1,1)/0.57735027_fp/, &
       (zrglx(2,l),l=1,2)/0.86113631_fp, 0.33998104_fp/, &
       (zrglx(3,l),l=1,3)/0.93246951_fp, 0.66120939_fp, 0.23861919_fp/, &
       (zrglx(4,l),l=1,4)/0.96028986_fp, 0.79666648_fp, 0.52553241_fp, &
       0.18343464_fp/, &
       (zrglx(5,l),l=1,5)/0.97390653_fp, 0.86506337_fp, 0.67940957_fp, &
       0.43339539_fp, 0.14887434_fp/, &
       (zrglx(6,l),l=1,6)/0.98156063_fp, 0.90411726_fp, 0.76990267_fp, &
       0.58731795_fp, 0.36783150_fp, 0.12523341_fp/, &
       (zrglx(7,l),l=1,7)/0.98628381_fp, 0.92843488_fp, 0.82720132_fp, &
       0.68729290_fp, 0.51524864_fp, 0.31911237_fp, 0.10805495_fp/
  !  gauss-legendre weights (see abramowitz and stegun).
  data zrglw(1,1)/1._fp/, &
       (zrglw(2,l),l=1,2)/ 0.34785485_fp, 0.65214515_fp/, &
       (zrglw(3,l),l=1,3)/0.17132449_fp, 0.36076157_fp, 0.46791393_fp/, &
       (zrglw(4,l),l=1,4)/0.10122854_fp, 0.22238103_fp, 0.31370665_fp, &
       0.36268378_fp/, &
       (zrglw(5,l),l=1,5)/0.06667134_fp, 0.14945135_fp, 0.21908636_fp, &
       0.26926672_fp, 0.29552422_fp/, &
       (zrglw(6,l),l=1,6)/0.04717534_fp, 0.10693933_fp, 0.16007833_fp, &
       0.20316743_fp, 0.23349254_fp, 0.24914705_fp/, &
       (zrglw(7,l),l=1,7)/0.03511946_fp, 0.08015809_fp, 0.12151857_fp, &
       0.15720317_fp, 0.18553840_fp, 0.20519846_fp, 0.21526385_fp/
  !
  !  kgh determines # of major subdivisions of the range of integration.
  !  kgl is the number of evaluation points in a major subdivision;
  !.     the fit is exact for polynomials up to order 2*kgl-1.
  !.     kgh=8 and kgl=10 is very accurate (4 subdivisions, 5th order legendre).
  !  kdim =-5 the complete integral is evaluated, =-4 dominant term only.
  kgh=8
  kgl=10
  kdim=-5
  sgvbt1=0._fp
  zvbeam=sqrt(2._fp*peb*evs/(pab*amucgs))
  if(pti.le.0._fp) then
     sgvbt1=zvbeam*aioniz(kzfast,katom,peb,pab,ksig)
     return
  end if
  zsgion=sqrt(2._fp*pti*evs/(pai*amucgs))
  !  set the range of integration.
  zvlowr=max(0._fp,zvbeam-4._fp*zsgion)
  zvuppr=zvbeam+4._fp*zsgion
  zdvint=(zvuppr-zvlowr)/real(kgh,fp)
  zsum=0._fp
  igh=kgh/2
  igl=kgl/2
  !  loop over the major subdivisions.
  do int=1,igh
     !        set the center of a major subdivision of the range of integration.
     zvcntr=zvlowr+zdvint*real(2*int-1,fp)
     do jlgl=1,kgl
        !          evaluate the integrand at the points in a major subdivision.
        if(jlgl.le.igl) then
           jgl=jlgl
           zvrel=zvcntr-zdvint*zrglx(igl,jgl)
        else
           jgl=jlgl-igl
           zvrel=zvcntr+zdvint*zrglx(igl,jgl)
        end if
        !          zebrel is in ev (it is not the center of mass energy):
        zebrel=pab*amucgs*zvrel**2*evsinv*0.5_fp
        zsig=aioniz(kzfast,katom,zebrel,pab,ksig)
        !          evaluate the exponentials in the integrand.
        zexp=exp(-((zvbeam-zvrel)/zsgion)**2)
        za2exp=((zvbeam+zvrel)/zsgion)**2
        !          test and include the smaller exponential.
        if((za2exp.lt.15._fp).and.(kdim.eq.-5)) zexp= &
             zexp-exp(-za2exp)
        !          accumulate the integrand with its gauss-legendre weighting.
        zsum=zsum+zsig*(zvrel**2*zexp)*zrglw(igl,jgl)
     end do
  end do
  !  normalize and return with the reactivity.
  sgvbt1=zsum*zdvint/(zsgion*zvbeam*sqrt(3.141593_fp))
  return
end function sgvbt1
! 11sep2004 fgtok
