program srspectra

! Make plots of synchrotron radiation spectra
! Use BMAD PHOTON_ENERGY_INIT to throw values.
! Compare to TOMS757 routine SYNCH1 by Allan McLeod,
! obtained from 
! http://people.sc.fsu.edu/~jburkardt/f_src/toms757/toms757.html
! Integration routine CADRE obtained from the INTLIB library here:
! http://people.sc.fsu.edu/~jburkardt/f_src/intlib/intlib.html
!
! The Physics of Synchrotron Radiation by Albert Hoffmann (pp. 82, 89, 90, 111, 112)
! ----------------------------------------------------------------------------------
!
! The total power radiated by an electron of relativistic factor
! gamma while accelerated with bending radius rho is
!  Ps = 2 * r0 * c**3 * m_e * gamma**4 / ( 3 * rho**2 )
!    where r0 is the classical electron radius
!
! Defining the critical frequency omega_c: 
!        omega_c = ( 3 / 2 ) * gamma**3 * c / rho
! we have the critical energy
!        E_c = hbar * omega_c = ( hbar * c ) * ( 3 / 2 ) * gamma**3 / rho
! with hbar * c = 197.327 MeV * fm
!
! The power radiated per unit photon energy E_gamma = hbar*omega is
!  dP/dE_gamma = ( 9 * sqrt(3) / 8pi ) * ( Ps / E_c ) * SYNCH1 (omega/omega_c) 
! i.e. SYNCH1 is dimensionless.
!
! The total photon rate is
!   R = dN/dt = ( 15 * sqrt(3) / 8 ) * Ps / E_c
! We use this relationship in SYNRAD output to get the critical energy from the ratio
! of the power radiated in a segment to the photon rate from that segment.
!
! The average photon energy is
!     <E_gamma> = Ps / (dN/dt) = ( 8 * sqrt(3) / 45 ) * E_c = 0.3079 * E_c
!
! The photon rate spectrum, i.e. the photon rate into an energy interval dE_gamma, is
!   dR / dE_gamma =  dP/dE_gamma / E_gamma  
! which is proportional to 
!                           SYNCH1(x) / x
!
! 18 august 2010 jac

use photon_init_mod

implicit none

real*8 synch1
external synch1
real*8 drde
external drde

real*8 a, b, abserr, relerr, error, result 
integer ind 

real*8 synch1int, drdeint

real(rp) E_rel
real*8 eoverec
real*8 synch,drdeval
real*8 synchnorm,drdenorm

real*8 e_relmin /1.e-4/
integer i

! Calculate the power and rate integrals for normalization
a=e_relmin
! b=20 did not get the full integral, so there was a value
! drdenorm=1.001 among the 100000 examples.
! b=30 gave IND=3 (danger)
b=25.
abserr=0.
relerr = 1d-6
! Get normalization integral for TOMS747 power spectrum function
call cadre ( synch1, a, b, abserr, relerr, error, result, ind )
synch1int = result
write(6,'(a,e15.6,a,e15.6,a,i5)') &
  ' Rtn CADRE returns SYNCH1 integral ',result,' +- ',error,'   IND=',ind
! Get normalization integral for photon rate function
call cadre ( drde, a, b, abserr, relerr, error, result, ind )
drdeint = result
write(6,'(a,e15.6,a,e15.6,a,i5)') &
  ' Rtn CADRE returns DRDE integral ',result,' +- ',error,'   IND=',ind

open (50, file = 'srspectra.out', status = 'unknown') 

do i=1,100000
! Throw E/Ec according to BMAD library generation function
 call photon_energy_init(E_rel)
! Protect against numerical instability in Fcn DRDE
 eoverec = E_rel
 if(eoverec.lt.e_relmin)cycle
! Value of TOMS747 power spectrum for this value of E/Ec
 synch = synch1(eoverec)
! Value of photon rate function for this value of E/Ec
 drdeval = drde(eoverec)

 b = eoverec
 call cadre ( synch1, a, b, abserr, relerr, error, result, ind )
 synchnorm = result / synch1int 
 call cadre ( drde, a, b, abserr, relerr, error, result, ind )
 drdenorm = result / drdeint 
 if(i.le.20)write(6,'(5(e15.6,2x))')eoverec,synch,drdeval,synchnorm,drdenorm
 write(50,'(5(e15.6,2x))')eoverec,synch,drdeval,synchnorm,drdenorm
enddo

close(50)

end program

real*8 function drde (x)
!
! Synchrotron radiation photon rate spectrum
!
! 20 Aug 2010 jac
!
real*8 synch1
external synch1
!
real*8 x,d0
!
d0=0.
if (x.gt.d0)then
 drde = synch1(x) / x
else
 drde = 0.
endif
!
end 

function synch1 ( xvalue )

!*****************************************************************************80
!
!! SYNCH1 calculates the synchrotron radiation function.
!
!  Discussion:
!
!    The function is defined by:
!
!      SYNCH1(x) = x * Integral ( x <= t < infinity ) K(5/3)(t) dt
!
!    where K(5/3) is a modified Bessel function of order 5/3.
!
!    The code uses Chebyshev expansions, the coefficients of which
!    are given to 20 decimal places.
!
!    This subroutine is set up to work on IEEE machines.
!
!  Modified:
!
!    07 September 2004
!
!  Author:
!
!    Allan McLeod,
!    Department of Mathematics and Statistics,
!    Paisley University, High Street, Paisley, Scotland, PA12BE
!    macl_ms0@paisley.ac.uk
!
!  Reference:
!
!    Allan McLeod,
!    Algorithm 757, MISCFUN: A software package to compute uncommon 
!      special functions,
!    ACM Transactions on Mathematical Software,
!    Volume 22, Number 3, September 1996, pages 288-301.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) XVALUE, the argument of the function.
!
!    Output, real ( kind = 8 ) SYNCH1, the value of the function.
!
  implicit none

  real ( kind = 8 ) cheval
  real ( kind = 8 ), parameter :: eight = 8.0D+00
  real ( kind = 8 ), parameter :: four = 4.0D+00
  real ( kind = 8 ), parameter :: half = 0.5D+00
  integer, parameter :: nterm1 = 12
  integer, parameter :: nterm2 = 10
  integer, parameter :: nterm3 = 21
  real ( kind = 8 ), parameter :: one = 1.0D+00
  real ( kind = 8 ) synch1
  real ( kind = 8 ), parameter :: three = 3.0D+00
  real ( kind = 8 ) x
  real ( kind = 8 ) xvalue
  real ( kind = 8 ), parameter :: zero = 0.0D+00

  real ( kind = 8 ) async1(0:13),async2(0:11),asynca(0:24), &
       cheb1,cheb2,conlow, &
       lnrtp2,pibrt3,t,twelve,xhigh1, &
       xhigh2,xlow,xpowth
  data twelve/ 12.0d0 /
  data conlow/2.14952824153447863671d0/
  data pibrt3/1.81379936423421785059d0/
  data lnrtp2/0.22579135264472743236d0/
  data async1/30.36468298250107627340d0, &
              17.07939527740839457449d0, &
               4.56013213354507288887d0, &
               0.54928124673041997963d0, &
               0.3729760750693011724d-1, &
               0.161362430201041242d-2, &
               0.4819167721203707d-4, &
               0.105124252889384d-5, &
               0.1746385046697d-7, &
               0.22815486544d-9, &
               0.240443082d-11, &
               0.2086588d-13, &
               0.15167d-15, &
               0.94d-18/
  data async2/0.44907216235326608443d0, &
              0.8983536779941872179d-1, &
              0.810445737721512894d-2, &
              0.42617169910891619d-3, &
              0.1476096312707460d-4, &
              0.36286336153998d-6, &
              0.666348074984d-8, &
              0.9490771655d-10, &
              0.107912491d-11, &
              0.1002201d-13, &
              0.7745d-16, &
              0.51d-18/
  data asynca(0)/ 2.13293051613550009848d0/
  data asynca(1)/ 0.7413528649542002401d-1/
  data asynca(2)/ 0.869680999099641978d-2/
  data asynca(3)/ 0.117038262487756921d-2/
  data asynca(4)/ 0.16451057986191915d-3/
  data asynca(5)/ 0.2402010214206403d-4/
  data asynca(6)/ 0.358277563893885d-5/
  data asynca(7)/ 0.54477476269837d-6/
  data asynca(8)/ 0.8388028561957d-7/
  data asynca(9)/ 0.1306988268416d-7/
  data asynca(10)/0.205309907144d-8/
  data asynca(11)/0.32518753688d-9/
  data asynca(12)/0.5179140412d-10/
  data asynca(13)/0.830029881d-11/
  data asynca(14)/0.133527277d-11/
  data asynca(15)/0.21591498d-12/
  data asynca(16)/0.3499673d-13/
  data asynca(17)/0.569942d-14/
  data asynca(18)/0.92906d-15/
  data asynca(19)/0.15222d-15/
  data asynca(20)/0.2491d-16/
  data asynca(21)/0.411d-17/
  data asynca(22)/0.67d-18/
  data asynca(23)/0.11d-18/
  data asynca(24)/0.2d-19/
!
!   Machine-dependent constants (suitable for IEEE machines)
!
  data xlow/2.98023224d-8/
  data xhigh1,xhigh2/809.595907d0,-708.396418d0/

  x = xvalue

  if ( x < zero ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SYNCH1 - Fatal error!'
    write ( *, '(a)' ) '  Argument X < 0.'
    synch1 = zero
  else if ( x < xlow ) then
    xpowth = x ** ( one / three )
    synch1 = conlow * xpowth
  else if ( x <= four ) then
    xpowth = x ** ( one / three )
    t = ( x * x / eight - half ) - half
    cheb1 = cheval ( nterm1, async1, t )
    cheb2 = cheval ( nterm2, async2, t )
    t = xpowth * cheb1 - xpowth**11 * cheb2
    synch1 = t - pibrt3 * x
  else if ( x <= xhigh1 ) then
    t = ( twelve - x ) / ( x + four )
    cheb1 = cheval ( nterm3, asynca, t )
    t = lnrtp2 - x + log ( sqrt ( x ) * cheb1 )
    if ( t < xhigh2 ) then
      synch1 = zero
    else
      synch1 = exp ( t )
    end if
  else
    synch1 = zero
  end if

  return
end

function cheval ( n, a, t )

!*****************************************************************************80
!
!! CHEVAL evaluates a Chebyshev series.
!
!  Discussion:
!
!    This function evaluates a Chebyshev series, using the
!    Clenshaw method with Reinsch modification, as analysed
!    in the paper by Oliver.
!
!  Modified:
!
!    07 August 2004
!
!  Author:
!
!    Allan McLeod,
!    Department of Mathematics and Statistics,
!    Paisley University, High Street, Paisley, Scotland, PA12BE
!    macl_ms0@paisley.ac.uk
!
!  Reference:
!
!    Allan McLeod,
!    Algorithm 757, MISCFUN: A software package to compute uncommon 
!      special functions,
!    ACM Transactions on Mathematical Software,
!    Volume 22, Number 3, September 1996, pages 288-301.
!
!    J Oliver,
!    An error analysis of the modified Clenshaw method for
!    evaluating Chebyshev and Fourier series,
!    Journal of the IMA, 
!    Volume 20, 1977, pages 379-391.
!
!  Parameters:
!
!    Input, integer N, the number of terms in the sequence.
!
!    Input, real ( kind = 8 ) A(0:N), the coefficients of the Chebyshev series.
!
!    Input, real ( kind = 8 ) T, the value at which the series is 
!    to be evaluated.
!
!    Output, real ( kind = 8 ) CHEVAL, the value of the Chebyshev series at T.
!
  implicit none

  integer n

  real ( kind = 8 ) :: a(0:n)
  real ( kind = 8 ) :: cheval
  real ( kind = 8 ) :: d1
  real ( kind = 8 ) :: d2
  real ( kind = 8 ), parameter :: half = 0.5D+00
  integer i
  real ( kind = 8 ) :: t
  real ( kind = 8 ), parameter :: test = 0.6D+00
  real ( kind = 8 ) :: tt
  real ( kind = 8 ), parameter :: two = 2.0D+00
  real ( kind = 8 ) :: u0
  real ( kind = 8 ) :: u1
  real ( kind = 8 ) :: u2
  real ( kind = 8 ), parameter :: zero = 0.0D+00

  u1 = zero
!
!  T <= -0.6, Reinsch modification.
!
  if ( t <= -test ) then

    d1 = zero
    tt = ( t + half ) + half
    tt = tt + tt

    do i = n, 0, -1
      d2 = d1
      u2 = u1
      d1 = tt * u2 + a(i) - d2
      u1 = d1 - u2
    end do

    cheval = ( d1 - d2 ) / two
!
!  -0.6 < T < 0.6, Standard Clenshaw method.
!
  else if ( t < test ) then

    u0 = zero
    tt = t + t

    do i = n, 0, -1
      u2 = u1
      u1 = u0
      u0 = tt * u1 + a(i) - u2
    end do

    cheval = ( u0 - u2 ) / two
!
!  0.6 <= T, Reinsch modification.
!
  else

    d1 = zero
    tt = ( t - half ) - half
    tt = tt + tt

    do i = n, 0, -1
      d2 = d1
      u2 = u1
      d1 = tt * u2 + a(i) + d2
      u1 = d1 + u2
    end do

    cheval = ( d1 + d2 ) / two

  end if

  return
end

subroutine cadre ( func, a, b, abserr, relerr, error, result, ind )

!*****************************************************************************80
!
!! CADRE estimates the integral of F(X) from A to B.
!
!  Discussion:
!
!    CADRE is the Cautious Adaptive Romberg Extrapolator.
!
!  Modified:
!
!    30 October 2000
!
!  Reference:
!
!    Philip Davis, Philip Rabinowitz,
!    Methods of Numerical Integration,
!    Second Edition,
!    Dover, 2007,
!    ISBN: 0486453391,
!    LC: QA299.3.D28.
!
!    Carl DeBoor, John Rice,
!    CADRE: An algorithm for numerical quadrature,
!    in Mathematical Software,
!    edited by John Rice,
!    Academic Press, 1971.
!    ISBN: 012587250X,
!    LC: QA1.M766,
!
!  Parameters:
!
!    Input, real ( kind = 8 ), external FUNC, the name of the function to be
!    integrated.  The user must declare the name an external parameter in the
!    calling program, write a function routine of the form 
!      FUNCTION FUNC ( X ) 
!    which evaluates the function at X, and pass the name of the function
!    in FUNC.
!
!    Input, real ( kind = 8 ) A, the lower limit of integration.
!
!    Input, real ( kind = 8 ) B, the upper limit of integration.
!
!    Input, real ( kind = 8 ) ABSERR, the absolute error tolerance.
!
!    Input, real ( kind = 8 ) RELERR, the relative error tolerance.
!
!    Output, real ( kind = 8 ) ERROR, an estimate of the absolute error.
!
!    Output, real ( kind = 8 ) RESULT, the approximate value of the integral.
!
!    Output, integer IND, reliability indicator.
!    If IND <= 2, RESULT is very reliable.  Higher values of
!    IND indicate less reliable values of RESULT.
!
  implicit none

  integer, parameter :: mxstge = 30
  integer, parameter :: maxtbl = 10
  integer, parameter :: maxts = 2049

  real ( kind = 8 ) a
  real ( kind = 8 ) abserr
  real ( kind = 8 ) ait(maxtbl)
  logical aitken
  real ( kind = 8 ), parameter :: aitlow = 1.1D+00
  real ( kind = 8 ), parameter :: aittol = 0.1D+00
  real ( kind = 8 ) astep
  real ( kind = 8 ) b
  real ( kind = 8 ) beg
  real ( kind = 8 ) begin(mxstge)
  real ( kind = 8 ) bma
  real ( kind = 8 ) curest
  real ( kind = 8 ) dif(maxtbl)
  real ( kind = 8 ) diff
  real ( kind = 8 ) end
  real ( kind = 8 ) ergoal
  real ( kind = 8 ) erra
  real ( kind = 8 ) errer
  real ( kind = 8 ) error
  real ( kind = 8 ) errr
  real ( kind = 8 ) est(mxstge)
  real ( kind = 8 ) fbeg
  real ( kind = 8 ) fbeg2
  real ( kind = 8 ) fend
  real ( kind = 8 ) fextm1
  real ( kind = 8 ) fextrp
  real ( kind = 8 ) finis(mxstge)
  real ( kind = 8 ) fn
  real ( kind = 8 ) fnsize
  real ( kind = 8 ), external :: func
  logical h2conv
  real ( kind = 8 ) h2next
  real ( kind = 8 ) h2tfex
  real ( kind = 8 ), parameter :: h2tol = 0.15D+00
  real ( kind = 8 ) hovn
  integer i
  integer ibeg
  integer ibegs(mxstge)
  integer iend
  integer ii
  integer iii
  integer ind
  integer istage
  integer istep
  integer istep2
  integer it
  integer l
  integer lm1
  integer n
  integer n2
  integer nnleft
  real ( kind = 8 ) prever
  real ( kind = 8 ) r(maxtbl)
  logical reglar
  logical reglsv(mxstge)
  real ( kind = 8 ) relerr
  real ( kind = 8 ) result
  logical right
  real ( kind = 8 ) rn(4)
  real ( kind = 8 ) rnderr
  real ( kind = 8 ) sing
  real ( kind = 8 ) singnx
  real ( kind = 8 ) slope
  real ( kind = 8 ) stage
  real ( kind = 8 ) step
  real ( kind = 8 ) stepmn
  real ( kind = 8 ) sum1
  real ( kind = 8 ) sumabs
  real ( kind = 8 ) t(maxtbl,maxtbl)
  real ( kind = 8 ) tabs
  real ( kind = 8 ) tabtlm
  real ( kind = 8 ), parameter :: tljump = 0.01D+00
  real ( kind = 8 ) ts(2049)
  real ( kind = 8 ) vint

  if ( a == b ) then
    result = 0.0D+00
    return
  end if
 
  begin(1:mxstge) = 0.0D+00
  est(1:mxstge) = 0.0D+00
  finis(1:mxstge) = 0.0D+00
  ibegs(1:mxstge) = 0
  reglsv(1:mxstge) = .false.
 
  vint = 0.0D+00
 
  rn(1:4) = (/ 0.7142005D+00, 0.3466282D+00, 0.8437510D+00, 0.1263305D+00 /)
 
  rnderr = epsilon ( rnderr )
  result = 0.0D+00
  error = 0.0D+00
  ind = 1
  bma = abs ( b - a )
  errr = min ( 0.1D+00, max ( abs ( relerr ), 10.0D+00*rnderr) )
  erra = abs ( abserr )
  stepmn = max ( bma / 2**mxstge, max ( bma, abs ( a ), abs ( b ) ) * rnderr )
  stage = 0.5D+00
  istage = 1
  curest = 0.0D+00
  fnsize = 0.0D+00
  prever = 0.0D+00
  reglar = .false.
  beg = a
  fbeg = func ( beg ) / 2.0D+00
  ts(1) = fbeg
  ibeg = 1
  end = b
  fend = func ( end ) / 2.0D+00
  ts(2) = fend
  iend = 2
 
10 continue
 
  right = .false.
 
20 continue

  step = end - beg
  astep = abs ( step )
 
  if ( astep < stepmn ) then
    ind = 5
    result = curest + vint
    return
  end if
 
  t(1,1) = fbeg + fend
  tabs = abs ( fbeg ) + abs ( fend )
  l = 1
  n = 1
  h2conv = .false.
  aitken = .false.
  go to 40
 
30 continue
 
40 continue
 
  lm1 = l
  l = l + 1
  n2 = n * 2
  fn = n2
  istep = ( iend - ibeg ) / n

  if ( 1 < istep ) then
    go to 60
  end if

  ii = iend
  iend = iend + n

  if ( maxts < iend ) then
    go to 440
  end if

  hovn = step / fn
 
  iii = iend
  do i = 1, n2, 2
    ts(iii) = ts(ii)
    ts(iii-1) = func ( end - i * hovn )
    iii = iii-2
    ii = ii-1
  end do
 
  istep = 2
 
60 continue
 
  istep2 = ibeg + istep / 2
 
  sum1 = 0.0D+00
  sumabs = 0.0D+00
  do i = istep2, iend, istep
    sum1 = sum1 + ts(i)
    sumabs = sumabs + abs ( ts(i) )
  end do
 
  t(l,1) = t(l-1,1) / 2.0D+00 + sum1 / fn
  tabs = tabs / 2.0D+00 + sumabs / fn
 
  n = n2
  it = 1
  vint = step * t(l,1)
  tabtlm = tabs * rnderr
  fnsize = max ( fnsize, abs ( t(l,1) ) )
  ergoal = max ( astep * rnderr * fnsize, &
    stage * max ( erra , errr * abs ( curest + vint ) ) )
  fextrp = 1.0D+00
  do i = 1, lm1
    fextrp = fextrp * 4.0D+00
    t(i,l) = t(l,i) - t(l-1,i)
    t(l,i+1) = t(l,i) + t(i,l) / ( fextrp - 1.0D+00 )
  end do
 
  errer = astep * abs ( t(1,l) )

  if ( 2 < l ) then
    go to 90
  end if

  if ( abs ( t(1,2) ) <= tabtlm ) then
    go to 290
  end if

  go to 40
 
90 continue
 
  do i = 2, lm1

    if ( tabtlm < abs ( t(i-1,l) ) ) then
      diff = t(i-1,lm1) / t(i-1,l)
    else
      diff = 0.0D+00
    end if

    t(i-1,lm1) = diff

  end do
 
  if ( abs ( 4.0D+00 - t(1,lm1) ) <= h2tol ) then
    go to 130
  end if

  if ( t(1,lm1) == 0.0D+00 ) then
    go to 120
  end if

  if ( abs ( 2.0D+00 - abs ( t(1,lm1) ) ) < tljump ) then
    go to 280
  end if

  if ( l == 3 ) then
    go to 30
  end if

  h2conv = .false.

  if ( abs ( ( t(1,lm1) - t(1,l-2) ) / t(1,lm1) ) <= aittol ) then
    go to 160
  end if
  
  if ( .not. reglar .and. l == 4 ) then
    go to 30
  end if
 
120 continue
 
  if ( errer <= ergoal ) then
    go to 310
  end if

  go to 380

130 continue

  if ( .not. h2conv ) then
    aitken = .false.
    h2conv = .true.
  end if

140 continue

  fextrp = 4.0D+00

150 continue

  it = it + 1
  vint = step * t(l,it)
  errer = abs ( step / ( fextrp - 1.0D+00 ) * t(it-1,l))

  if ( errer <= ergoal ) then
    go to 340
  end if

  if ( it == lm1 ) then
    go to 270
  end if

  if ( t(it,lm1) == 0.0D+00 ) then
    go to 150
  end if

  if ( t(it,lm1) <= fextrp ) then
    go to 270
  end if

  if ( abs ( t(it,lm1) / 4.0D+00 - fextrp ) / fextrp < aittol ) then
    fextrp = fextrp * 4.0D+00
  end if

  go to 150
 
160 continue

  if ( t(1,lm1) < aitlow ) then
    go to 380
  end if
 
  if ( .not. aitken ) then
    h2conv = .false.
    aitken = .true.
  end if
 
170 continue

  fextrp = t(l-2,lm1)

  if ( 4.5D+00 < fextrp ) then
    go to 140
  end if

  if ( fextrp < aitlow ) then
    go to 380
  end if

  if ( h2tol < abs ( fextrp - t(l-3,lm1) ) / t(1,lm1) ) then
    go to 380
  end if

  sing = fextrp
  fextm1 = fextrp - 1.0D+00

  ait(1) = 0.0D+00
  do i = 2, l
    ait(i) = t(i,1) + (t(i,1)-t(i-1,1)) / fextm1
    r(i) = t(1,i-1)
    dif(i) = ait(i) - ait(i-1)
  end do

  it = 2

190 continue

  vint = step * ait(l)

200 continue

  errer = errer / fextm1
 
  if ( errer <= ergoal ) then
    ind = max ( ind, 2 )
    go to 340
  end if
 
210 continue

  it = it + 1

  if ( it == lm1 ) then
    go to 270
  end if

  if ( it <= 3 ) then
    h2next = 4.0D+00
    singnx = 2.0D+00 * sing
  end if

  if ( h2next < singnx ) then
    go to 230
  end if

  fextrp = singnx
  singnx = 2.0D+00 * singnx
  go to 240

230 continue

  fextrp = h2next
  h2next = 4.0D+00 * h2next

240 continue
 
  do i = it, lm1
    if ( tabtlm < abs ( dif(i+1) ) ) then
      r(i+1) = dif(i) / dif(i+1)
    else
      r(i+1) = 0.0D+00
    end if
  end do
 
  h2tfex = -h2tol * fextrp

  if ( r(l) - fextrp < h2tfex ) then
    go to 270
  end if

  if ( r(l-1) - fextrp < h2tfex ) then
    go to 270
  end if

  errer = astep * abs ( dif(l) )
  fextm1 = fextrp - 1.0D+00
  do i = it, l
    ait(i) = ait(i)+dif(i) / fextm1
    dif(i) = ait(i)-ait(i-1)
  end do
 
  go to 190
 
270 continue

  fextrp = max ( prever / errer, aitlow )
  prever = errer
  if ( l < 5 ) then
    go to 40
  end if

  if ( 2.0D+00 < l - it .and. istage < mxstge ) then
    go to 370
  end if

  if ( errer / fextrp**( maxtbl - l ) < ergoal ) then
    go to 40
  end if

  go to 370
 
280 continue

  if ( ergoal < errer ) then
    go to 370
  end if

  diff = abs ( t(1,l) ) * 2.0D+00 * fn
  go to 340
 
290 continue

  slope = ( fend - fbeg ) * 2.0D+00
  fbeg2 = fbeg * 2.0D+00
 
  do i = 1, 4
    diff = abs ( func ( beg + rn(i) * step ) - fbeg2 - rn(i) * slope )
    if ( tabtlm < diff ) then
      go to 330
    end if
  end do
 
  go to 340
 
310 continue

  slope = ( fend - fbeg ) * 2.0D+00
  fbeg2 = fbeg * 2.0D+00
  i = 1
 
320 continue

  diff = abs ( func ( beg + rn(i) * step ) - fbeg2 - rn(i) * slope )
 
330 continue

  errer = max ( errer, astep * diff )

  if ( ergoal < errer ) then
    go to 380
  end if

  i = i+1

  if ( i <= 4 ) then
    go to 320
  end if

  ind = 3
 
340 continue

  result = result + vint
  error = error + errer
 
350 continue

  if ( right ) then
    go to 360
  end if

  istage = istage - 1

  if ( istage == 0 ) then
    return
  end if

  reglar = reglsv(istage)
  beg = begin(istage)
  end = finis(istage)
  curest = curest - est(istage+1) + vint
  iend = ibeg - 1
  fend = ts(iend)
  ibeg = ibegs(istage)
  go to 400
 
360 continue

  curest = curest + vint
  stage = stage * 2.0D+00
  iend = ibeg
  ibeg = ibegs(istage)
  end = beg
  beg = begin(istage)
  fend = fbeg
  fbeg = ts(ibeg)
  go to 10
 
370 continue

  reglar = .true.
 
380 continue
 
  if ( istage == mxstge ) then
    ind = 5
    result = curest + vint
    return
  end if
 
390 continue

  if ( right ) then
    go to 410
  end if

  reglsv(istage+1) = reglar
  begin(istage) = beg
  ibegs(istage) = ibeg
  stage = stage / 2.0D+00

400 continue

  right = .true.
  beg = ( beg + end ) / 2.0D+00
  ibeg = ( ibeg + iend ) / 2
  ts(ibeg) = ts(ibeg) / 2.0D+00
  fbeg = ts(ibeg)
  go to 20

410 continue

  nnleft = ibeg - ibegs(istage)
  if ( maxts <= end + nnleft ) then
    go to 440
  end if

  iii = ibegs(istage)
  ii = iend
  do i = iii, ibeg
    ii = ii + 1
    ts(ii) = ts(i)
  end do
 
  do i = ibeg, ii
    ts(iii) = ts(i)
    iii = iii + 1
  end do
 
  iend = iend + 1
  ibeg = iend - nnleft
  fend = fbeg
  fbeg = ts(ibeg)
  finis(istage) = end
  end = beg
  beg = begin(istage)
  begin(istage) = end
  reglsv(istage) = reglar
  istage = istage + 1
  reglar = reglsv(istage)
  est(istage) = vint
  curest = curest + est(istage)
  go to 10

440 continue

  ind = 4

460 continue

  result = curest + vint

  return
end