positions Subroutine

    subroutine reserve(self, n, positions, attitudes, colors)
        class(Trajectory), intent(inout) :: self
        integer, intent(in) :: n                            !! number of points
        logical, intent(in), optional :: positions          !! reserve positions (default true)
        logical, intent(in), optional :: attitudes          !! reserve attitudes (default false)
        logical, intent(in), optional :: colors             !! reserve colors (default false)
        logical :: pos, att, col
        pos = .true.
        att = .false.
        col = .false.
        if (present(positions)) pos = positions
        if (present(attitudes)) att = attitudes
        if (present(colors)) col = colors
        call grow_times(self, n, exact=.true.)
        if (pos) call grow_positions(self, n, exact=.true.)
        if (att) call grow_attitudes(self, n, exact=.true.)
        if (col) call grow_colors(self, n, exact=.true.)
    end subroutine reserve

    !> Free the storage not used: each array is reallocated to the points it
    !> holds (one copy, an array at a time), and arrays holding none are
    !> deallocated. For a trajectory that is complete: growing doubles the
    !> storage, so up to about half of it can be unused. The data and what is
    !> drawn of them don't change; pointers from [[times]], [[positions]], ...
    !> must be fetched again. Points can still be added (storage grows again).
    subroutine compact(self)
        class(Trajectory), intent(inout) :: self
        real(wp), allocatable :: t(:), posopt(:,:,:), att(:,:)
        real(c_float), allocatable :: col(:,:)
        if (allocated(self%t)) then
            if (size(self%t) /= self%n_times) then
                allocate(t(self%n_times))
                t = self%t(1:self%n_times)
                call move_alloc(t, self%t)
            end if
            if (self%n_times == 0) deallocate(self%t)
        end if
        if (allocated(self%posopt)) then
            if (size(self%posopt, 3) /= self%n_pos) then
                allocate(posopt(3, 0:self%nopt_, self%n_pos))
                posopt = self%posopt(:, :, 1:self%n_pos)
                call move_alloc(posopt, self%posopt)
            end if
            if (self%n_pos == 0) deallocate(self%posopt)
        end if
        if (allocated(self%att)) then
            if (size(self%att, 2) /= self%n_att) then
                allocate(att(4, self%n_att))
                att = self%att(:, 1:self%n_att)
                call move_alloc(att, self%att)
            end if
            if (self%n_att == 0) deallocate(self%att)
        end if
        if (allocated(self%col)) then
            if (size(self%col, 2) /= self%n_col) then
                allocate(col(4, self%n_col))
                col = self%col(:, 1:self%n_col)
                call move_alloc(col, self%col)
            end if
            if (self%n_col == 0) deallocate(self%col)
        end if
    end subroutine compact

    !> The memory the trajectory's data take (bytes), used or not: its
    !> storage's capacity (see [[reserve]], [[compact]]).
    pure integer(int64) function storage_bytes(self) result(b)
        class(Trajectory), intent(in) :: self
        b = 0_int64
        if (allocated(self%t)) b = b + storage_size(self%t, int64) / 8 * size(self%t, kind=int64)
        if (allocated(self%posopt)) b = b + storage_size(self%posopt, int64) / 8 * size(self%posopt, kind=int64)
        if (allocated(self%att)) b = b + storage_size(self%att, int64) / 8 * size(self%att, kind=int64)
        if (allocated(self%col)) b = b + storage_size(self%col, int64) / 8 * size(self%col, kind=int64)
    end function storage_bytes

    !> Remove all points (keeping the memory).
    subroutine clear(self)
        class(Trajectory), intent(inout) :: self
        self%n_times = 0
        self%n_pos = 0
        self%n_att = 0
        self%n_col = 0
        self%generation_ = self%generation_ + 1
        self%revision_ = self%revision_ + 1
    end subroutine clear

    !> Set the degrees of freedom (2 or 3) of positions and optionals; clears the trajectory.
    subroutine set_dof(self, dof)
        class(Trajectory), intent(inout) :: self
        integer, intent(in) :: dof                          !! degrees of freedom (2 or 3)
        self%dof_ = max(2, min(3, dof))
        call self%clear()
    end subroutine set_dof

    !> Degrees of freedom (2 or 3) of positions and optionals.
    pure integer function dof(self)
        class(Trajectory), intent(in) :: self
        dof = self%dof_
    end function dof

    !> Set the number of optional vectors per position; clears the trajectory.
    subroutine set_num_optionals(self, nopt)
        class(Trajectory), intent(inout) :: self
        integer, intent(in) :: nopt                         !! optional vectors per position (>= 0)
        if (nopt /= self%nopt_ .and. allocated(self%posopt)) deallocate(self%posopt)
        self%nopt_ = max(0, nopt)
        call self%clear()
    end subroutine set_num_optionals

    !> Number of optional vectors per position.
    pure integer function num_optionals(self)
        class(Trajectory), intent(in) :: self
        num_optionals = self%nopt_
    end function num_optionals

    !> Number of times.
    pure integer function num_times(self)
        class(Trajectory), intent(in) :: self
        num_times = self%n_times
    end function num_times

    !> Number of positions (at most num_times).
    pure integer function num_positions(self)
        class(Trajectory), intent(in) :: self
        num_positions = self%n_pos
    end function num_positions

    !> Number of attitudes (at most num_times).
    pure integer function num_attitudes(self)
        class(Trajectory), intent(in) :: self
        num_attitudes = self%n_att
    end function num_attitudes

    !> Number of colors (at most num_times).
    pure integer function num_colors(self)
        class(Trajectory), intent(in) :: self
        num_colors = self%n_col
    end function num_colors

    !> True if the trajectory has no times.
    pure logical function is_empty(self)
        class(Trajectory), intent(in) :: self
        is_empty = self%n_times == 0
    end function is_empty

    !> Number of changes so far (any added, set or cleared data).
    pure integer function revision(self)
        class(Trajectory), intent(in) :: self
        revision = self%revision_
    end function revision

    !> Number of times the data were cleared (by [[clear]], [[set_dof]] or
    !> [[set_num_optionals]]): while it is unchanged, data were only appended.
    pure integer function generation(self)
        class(Trajectory), intent(in) :: self
        generation = self%generation_
    end function generation

    ! ==================================================================
    ! adding data

    !> Add a time at the end of the time list.
    subroutine add_time(self, t)
        class(Trajectory), intent(inout) :: self
        real(wp), intent(in) :: t                           !! the time
        call grow_times(self, self%n_times + 1)
        self%n_times = self%n_times + 1
        self%t(self%n_times) = t
        self%revision_ = self%revision_ + 1
    end subroutine add_time

    !> Add a position (the first [[dof]] elements of `pos` are used) for the next
    !> time that has none. `ok` is false if every time already has a position.
    subroutine add_position(self, pos, ok)
        class(Trajectory), intent(inout) :: self
        real(wp), intent(in) :: pos(:)          !! position (first [[dof]] elements used)
        logical, intent(out), optional :: ok    !! false if not added
        integer :: d
        if (present(ok)) ok = .false.
        if (self%n_pos >= self%n_times) return
        call grow_positions(self, self%n_pos + 1)
        self%n_pos = self%n_pos + 1
        d = min(self%dof_, size(pos))
        self%posopt(:, :, self%n_pos) = 0.0_wp
        self%posopt(1:d, 0, self%n_pos) = pos(1:d)
        self%revision_ = self%revision_ + 1
        if (present(ok)) ok = .true.
    end subroutine add_position

    !> Add an attitude quaternion [x, y, z, w] for the next time that has none.
    !> `ok` is false if every time already has an attitude.
    subroutine add_attitude(self, q, ok)
        class(Trajectory), intent(inout) :: self
        real(wp), intent(in) :: q(4)            !! quaternion [x, y, z, w]
        logical, intent(out), optional :: ok    !! false if not added
        if (present(ok)) ok = .false.
        if (self%n_att >= self%n_times) return
        call grow_attitudes(self, self%n_att + 1)
        self%n_att = self%n_att + 1
        self%att(:, self%n_att) = q
        self%revision_ = self%revision_ + 1
        if (present(ok)) ok = .true.
    end subroutine add_attitude

    !> Add a color ([r, g, b], alpha 1, or [r, g, b, a], 0 to 1) for the next
    !> time that has none. `ok` is false if every time already has a color.
    subroutine add_color(self, color, ok)
        class(Trajectory), intent(inout) :: self
        real(wp), intent(in) :: color(:)        !! [r, g, b] or [r, g, b, a] (0 to 1)
        logical, intent(out), optional :: ok    !! false if not added
        if (present(ok)) ok = .false.
        if (self%n_col >= self%n_times) return
        call grow_colors(self, self%n_col + 1)
        self%n_col = self%n_col + 1
        call to_float_rgba(self%col(:, self%n_col), color)
        self%revision_ = self%revision_ + 1
        if (present(ok)) ok = .true.
    end subroutine add_color

    !> rgba from [r, g, b] (alpha 1) or [r, g, b, a].
    pure subroutine to_float_rgba(rgba, color)
        real(c_float), intent(out) :: rgba(4)   !! [r, g, b, a]
        real(wp), intent(in) :: color(:)        !! [r, g, b] or [r, g, b, a] (0 to 1)
        integer :: n
        n = min(4, size(color))
        rgba = 1.0_c_float
        rgba(1:n) = real(color(1:n), c_float)
    end subroutine to_float_rgba

    !> Set optional vector `index` (1 to num_optionals) of the most recently
    !> added position. `ok` is false if there is no position or no such optional.
    subroutine set_optional(self, index, v, ok)
        class(Trajectory), intent(inout) :: self
        integer, intent(in) :: index            !! the optional (1 to num_optionals)
        real(wp), intent(in) :: v(:)            !! the vector (its first [[dof]] elements are used)
        logical, intent(out), optional :: ok    !! false if not set
        integer :: d
        if (present(ok)) ok = .false.
        if (index < 1 .or. index > self%nopt_ .or. self%n_pos == 0) return
        d = min(self%dof_, size(v))
        self%posopt(1:d, index, self%n_pos) = v(1:d)
        self%revision_ = self%revision_ + 1
        if (present(ok)) ok = .true.
    end subroutine set_optional

    !> Add a point: a time with an optional position, attitude and color
    !> ([r, g, b] or [r, g, b, a]).
    !>
    !> A circular orbit, a point a minute (with `orbit%add`, the generic name):
    !>
    !> ```fortran
    !> orbit = Trajectory()
    !> do k = 0, 360
    !>     t = 60.0_wp * k
    !>     a = 2.0_wp * pi * t / 5400.0_wp
    !>     call orbit%add(t, pos=7000.0_wp * [cos(a), sin(a), 0.0_wp])
    !> end do
    !> ```
    subroutine add_point(self, t, pos, att, color)
        class(Trajectory), intent(inout) :: self
        real(wp), intent(in) :: t                       !! time
        real(wp), intent(in), optional :: pos(:)        !! position (first [[dof]] elements used)
        real(wp), intent(in), optional :: att(4)        !! attitude quaternion [x, y, z, w]
        real(wp), intent(in), optional :: color(:)      !! [r, g, b] or [r, g, b, a] (0 to 1)
        call self%add_time(t)
        if (present(pos)) call self%add_position(pos)
        if (present(att)) call self%add_attitude(att)
        if (present(color)) call self%add_color(color)
    end subroutine add_point

    !> Add points in bulk: times `t(n)`, with optional positions `pos(dof, n)`,
    !> attitudes `att(4, n)` and colors `color(3 or 4, n)`. Times, positions,
    !> attitudes and colors must all have been complete before (as many of each
    !> as times) for the given ones to line up with `t`.
    !>
    !> Points computed elsewhere, all at once (`times(n)`, `positions(3, n)`):
    !>
    !> ```fortran
    !> call orbit%add(times, pos=positions)
    !> ```
    subroutine add_points(self, t, pos, att, color)
        class(Trajectory), intent(inout) :: self
        real(wp), intent(in) :: t(:)                    !! (n): times
        real(wp), intent(in), optional :: pos(:,:)      !! (dof, n): positions
        real(wp), intent(in), optional :: att(:,:)      !! (4, n): attitude quaternions [x, y, z, w]
        real(wp), intent(in), optional :: color(:,:)    !! (3 or 4, n): colors [r, g, b(, a)] (0 to 1)
        integer :: n, i0, d
        n = size(t)
        i0 = self%n_times
        call grow_times(self, i0 + n)
        self%t(i0+1:i0+n) = t
        self%n_times = i0 + n
        if (present(pos)) then
            n = min(size(pos, 2), self%n_times - self%n_pos)
            d = min(self%dof_, size(pos, 1))
            call grow_positions(self, self%n_pos + n)
            self%posopt(:, :, self%n_pos+1:self%n_pos+n) = 0.0_wp
            self%posopt(1:d, 0, self%n_pos+1:self%n_pos+n) = pos(1:d, 1:n)
            self%n_pos = self%n_pos + n
        end if
        if (present(att)) then
            n = min(size(att, 2), self%n_times - self%n_att)
            call grow_attitudes(self, self%n_att + n)
            self%att(:, self%n_att+1:self%n_att+n) = att(:, 1:n)
            self%n_att = self%n_att + n
        end if
        if (present(color)) then
            n = min(size(color, 2), self%n_times - self%n_col)
            d = min(4, size(color, 1))
            call grow_colors(self, self%n_col + n)
            self%col(:, self%n_col+1:self%n_col+n) = 1.0_c_float
            self%col(1:d, self%n_col+1:self%n_col+n) = real(color(1:d, 1:n), c_float)
            self%n_col = self%n_col + n
        end if
        self%revision_ = self%revision_ + 1
    end subroutine add_points

    ! ==================================================================
    ! reading data

    !> Time i.
    pure real(wp) function time(self, i)
        class(Trajectory), intent(in) :: self
        integer, intent(in) :: i                            !! the time (1 to num_times)
        time = self%t(i)
    end function time

    !> Position i (z = 0 for a 2-DOF trajectory).
    pure function position(self, i) result(p)
        class(Trajectory), intent(in) :: self
        integer, intent(in) :: i                            !! the position (1 to num_positions)
        real(wp) :: p(3)
        p = self%posopt(:, 0, i)
    end function position

    !> Attitude i (quaternion [x, y, z, w]).
    pure function attitude(self, i) result(q)
        class(Trajectory), intent(in) :: self
        integer, intent(in) :: i                            !! the attitude (1 to num_attitudes)
        real(wp) :: q(4)
        q = self%att(:, i)
    end function attitude

    !> Color i ([r, g, b, a]).
    pure function color(self, i) result(c)
        class(Trajectory), intent(in) :: self
        integer, intent(in) :: i                            !! the color (1 to num_colors)
        real(wp) :: c(4)
        c = real(self%col(:, i), wp)
    end function color

    !> Optional vector k (1 to num_optionals) of position i.
    pure function optional(self, i, k) result(v)
        class(Trajectory), intent(in) :: self
        integer, intent(in) :: i, k                         !! the position (1 to num_positions); the optional (1 to num_optionals)
        real(wp) :: v(3)
        v = self%posopt(:, k, i)
    end function optional

    !> Pointer to the times (num_times), without copying.
    subroutine times(self, t)
        class(Trajectory), intent(in), target :: self
        real(wp), dimension(:), pointer, intent(out) :: t   !! (num_times): times (null if none)
        if (self%n_times == 0) then
            t => null()
        else
            t => self%t(1:self%n_times)
        end if
    end subroutine times

    !> Pointer to the positions (3, num_positions), without copying.
    subroutine positions(self, p)
        class(Trajectory), intent(in), target :: self
        real(wp), dimension(:,:), pointer, intent(out) :: p !! (3, num_positions): positions (null if none)
        if (self%n_pos == 0) then
            p => null()
        else
            p => self%posopt(:, 0, 1:self%n_pos)
        end if
    end subroutine positions