"""
    track_propagated_discontinuities!(integrator::DDEIntegrator)

Try to find a propagated discontinuity in the time interval `[integrator.t, integrator.t +
integrator.dt]` and add it to the set of discontinuities and grid points of the
`integrator`.
"""
function track_propagated_discontinuities!(integrator::DDEIntegrator)
  # calculate interpolation points
  interp_points = integrator.discontinuity_interp_points
  Θs = range(zero(integrator.t); stop = oneunit(integrator.t), length = interp_points)

  # for dependent lags and previous discontinuities
  for lag in integrator.sol.prob.dependent_lags, discontinuity in integrator.tracked_discontinuities
    # obtain time of previous discontinuity
    T = discontinuity.t

    # estimate subinterval of current integration step that contains a propagated
    # discontinuity induced by the lag and the previous discontinuity
    interval = discontinuity_interval(integrator, lag, T, Θs)

    # if a discontinuity exists in the current integration step
    if interval !== nothing
      # estimate time point of discontinuity
      t = discontinuity_time(integrator, lag, T, interval)

      # add new discontinuity of correct order at the estimated time point
      if integrator.sol.prob.neutral
        d = Discontinuity(t, discontinuity.order)
      else
        d = Discontinuity(t, discontinuity.order + 1)
      end
      push!(integrator.opts.d_discontinuities, d)
      push!(integrator.opts.tstops, t)

      # analogously to RADAR5 we do not strive for finding the first discontinuity
      break
    end
  end

  nothing
end

"""
    discontinuity_function(integrator::DDEIntegrator, lag, T, t)

Evaluate function ``f(x) = T + lag(u(x), p, x) - x`` at time point `t`, where `T` is time
point of a previous discontinuity and `lag` is a dependent delay.
"""
function discontinuity_function(integrator::DDEIntegrator, lag, T, t)
  tmp = get_tmp_cache(integrator)
  cache = tmp === nothing ? nothing : first(tmp)

  # estimate state at the given time
  if cache === nothing
    ut = integrator(t, Val{0})
  else
    integrator(cache, t, Val{0})
    ut = cache
  end

  T + lag(ut, integrator.p, t) - t
end

"""
    discontinuity_interval(integrator::DDEIntegrator, lag, T, Θs)

Return an estimated subinterval of the current integration step of the `integrator` that
contains a propagated discontinuity induced by the dependent delay `lag` and the
discontinuity at time point `T`, or `nothing`.

The interval is estimated by checking the signs of `T + lag(u(t), p, t) - t` for time points
`integrator.t .+ θs` in the interval `[integrator.t, integrator.t + integrator.dt]`.
"""
function discontinuity_interval(integrator::DDEIntegrator, lag, T, Θs)
  # use start and end point of last time interval to check for discontinuities
  previous_condition = discontinuity_function(integrator, lag, T, integrator.t)
  if isapprox(previous_condition, 0;
              rtol = integrator.discontinuity_reltol,
              atol = integrator.discontinuity_abstol)
    prev_sign = 0
  else
    prev_sign = cmp(previous_condition, zero(previous_condition))
  end
  new_condition = discontinuity_function(integrator, lag, T, integrator.t + integrator.dt)
  new_sign = cmp(new_condition, zero(new_condition))

  # if signs are different we already know that a discontinuity exists
  if prev_sign * new_sign < 0
    return (zero(eltype(Θs)), one(eltype(Θs)))
  end

  # recheck interpolation intervals if no discontinuity found yet
  prev_Θ = zero(eltype(Θs))

  for i in 2:length(Θs)
    # evaluate sign at next interpolation point
    new_Θ = Θs[i]
    new_t = integrator.t + new_Θ * integrator.dt
    new_condition = discontinuity_function(integrator, lag, T, new_t)
    new_sign = cmp(new_condition, zero(new_condition))

    # return estimated interval if we find a root or observe a switch of signs
    if new_sign == 0
      return (new_Θ, new_Θ)
    elseif prev_sign * new_sign < 0
      return (prev_Θ, new_Θ)
    else
      # otherwise update lower estimate of subinterval
      prev_sign = new_sign
      prev_Θ = new_Θ
    end
  end

  nothing
end

"""
    discontinuity_time(integrator::DDEIntegrator, lag, T, interval)

Estimate time point of the propagated discontinuity induced by the dependent delay
`lag` and the discontinuity at time point `T` inside the `interval` of the current
integration step of the `integrator`.
"""
function discontinuity_time(integrator::DDEIntegrator, lag, T, (bottom_Θ, top_Θ))
  if bottom_Θ == top_Θ
    # in that case we have already found the time point of a discontinuity
    Θ = top_Θ
  else
    # define function for root finding
    zero_func = let integrator = integrator, lag = lag, T = T, t = integrator.t, dt = integrator.dt
      (θ, p=nothing) -> discontinuity_function(integrator, lag, T, t + θ * dt)
    end

    Θ = NonlinearSolve.solve(
      NonlinearSolve.NonlinearProblem{false}(zero_func, (bottom_Θ, top_Θ)),
      NonlinearSolve.Falsi(),
    ).left
  end

  integrator.t + Θ * integrator.dt
end
