draw_trail Subroutine

    subroutine draw_trail(self, t)
        type(CurveArtist), intent(inout) :: self
        real(wp), intent(in) :: t                           !! the current time
        real(wp), parameter :: unlimited = 0.5_wp * huge(1.0_wp)
        real(c_float), parameter :: far = 1.0e30_c_float
        real(c_float) :: values(20), hi
        real(wp) :: sgn, lo_key, hi_key, s_now, fade
        integer :: n, i1, i2, k
        call check_status(vkf_line_set_range(self%tail, 1_c_int, 0_c_int), 'trajectory artists')
        self%tail_shown = .false.
        n = self%n_written
        self%n_drawn = 0
        self%first_drawn_ = 1
        values = 0.0_c_float
        if (.not. associated(self%traj) .or. n < 2) then
            call self%chain%draw_range(1, 0)
            call self%chain%set_trail(values)
            return
        end if
        sgn = merge(-1.0_wp, 1.0_wp, decreasing_times(self%traj, n))   ! (keys sgn * time increase with the index)
        lo_key = sgn * t - merge(huge(1.0_wp), self%trail_%past, self%trail_%past >= unlimited)
        hi_key = sgn * t + merge(huge(1.0_wp), self%trail_%future, self%trail_%future >= unlimited)
        if (sgn < 0.0_wp) then                              ! (the past is at the larger keys)
            lo_key = sgn * t - merge(huge(1.0_wp), self%trail_%future, self%trail_%future >= unlimited)
            hi_key = sgn * t + merge(huge(1.0_wp), self%trail_%past, self%trail_%past >= unlimited)
        end if
        if (hi_key < key(1) .or. lo_key > key(n)) then     ! (the window misses the trajectory)
            call self%chain%draw_range(1, 0)
            call self%chain%set_trail(values)
            return
        end if
        i1 = max(1, keys_up_to(lo_key, .true.))                   ! the last point before the window, or the first
        i2 = min(n, keys_up_to(hi_key, .false.) + 1)                 ! the first after it, or the last
        self%first_drawn_ = i1
        self%n_drawn = i2 - i1 + 1
        call self%chain%draw_range(i1, i2)
        s_now = s_at(t)
        hi = real(s_now, c_float)
        fade = 1.0_wp - self%trail_%fade
        values(1:4) = [hi, real(s_now - real(hi, wp), c_float), real(sgn, c_float), &
                       real(merge(1, 0, self%trail_%replace(1)) + merge(2, 0, self%trail_%replace(2)), c_float)]
        if (self%trail_%past >= unlimited) then
            values(5:6) = -far
        else
            values(5) = real(sgn * (s_at(t - self%trail_%past) - s_now), c_float)
            values(6) = real(sgn * (s_at(t - fade * self%trail_%past) - s_now), c_float)
        end if
        if (self%trail_%future >= unlimited) then
            values(7:8) = far
        else
            values(7) = real(sgn * (s_at(t + fade * self%trail_%future) - s_now), c_float)
            values(8) = real(sgn * (s_at(t + self%trail_%future) - s_now), c_float)
        end if
        do k = 1, 2                                         ! each part's color (or the line's: then its alpha only)
            values(9 + 4 * (k - 1):12 + 4 * (k - 1)) = real(self%trail_%color(:, k), c_float)
            values(12 + 4 * (k - 1)) = real(self%trail_%opacity(k) * merge(self%trail_%color(4, k), 1.0_wp, &
                                                                           self%trail_%replace(k)), c_float)
        end do
        values(17:20) = real([self%trail_%pattern(1), self%trail_%factor(1), self%trail_%pattern(2), &
                              self%trail_%factor(2)], c_float)
        call self%chain%set_trail(values)

    contains

        !> Point i's key: its time, increasing with the index.
        real(wp) function key(i)
            integer, intent(in) :: i                        !! the point
            key = sgn * self%traj%time(i)
        end function key

        !> The number of points with keys below x, or at most x if `inclusive`
        !> (bisection: the keys increase).
        integer function keys_up_to(x, inclusive) result(m)
            real(wp), intent(in) :: x                       !! the key
            logical, intent(in) :: inclusive                !! keys equal to x counted
            integer :: low, high, mid
            low = 0
            high = n + 1
            do while (high - low > 1)
                mid = (low + high) / 2
                if (key(mid) < x .or. (inclusive .and. key(mid) == x)) then
                    low = mid
                else
                    high = mid
                end if
            end do
            m = low
        end function keys_up_to

        !> The distance along the curve at time tt: interpolated between the
        !> points about it (extrapolated beyond the ends).
        real(wp) function s_at(tt)
            real(wp), intent(in) :: tt                      !! the time
            integer :: j
            real(wp) :: t0, t1
            j = max(1, min(n - 1, keys_up_to(sgn * tt, .true.)))
            t0 = self%traj%time(j)
            t1 = self%traj%time(j + 1)
            s_at = self%chain%distance(j)
            if (t1 /= t0) s_at = s_at + (tt - t0) / (t1 - t0) * (self%chain%distance(j + 1) - s_at)
        end function s_at

    end subroutine draw_trail