two_sided Function

    function surface_node(xyz, idx, rgba, normals, colors, two_sided, sphere) result(node)
        real(wp), intent(in) :: xyz(:,:)                    !! (3, n) vertices
        integer, intent(in) :: idx(:,:)                     !! (3, m) triangles, vertex numbers (1-based)
        real(wp), intent(in) :: rgba(4)                     !! color [r, g, b, a] (0 to 1)
        real(wp), intent(in), optional :: normals(:,:)      !! (3, n): lit
        real(wp), intent(in), optional :: colors(:,:)       !! (4, n)
        logical, intent(in), optional :: two_sided          !! (lit) both faces (default: front only)
        real(wp), intent(in), optional :: sphere(4)         !! (lit, translucent) center and radius
        type(c_ptr) :: node
        real(c_float), dimension(:,:), pointer :: v, n, uv, c
        integer(c_int32_t), dimension(:), pointer :: indices
        integer :: nv, ni, flags, i
        logical :: lit, translucent
        node = c_null_ptr
        nv = size(xyz, 2)
        ni = size(idx)
        if (nv < 3 .or. ni < 3) return
        lit = present(normals)
        flags = 1
        if (lit) then
            translucent = rgba(4) < 1.0_wp
            if (present(colors)) translucent = translucent .or. any(colors(4, :) < 1.0_wp)
            flags = 0
            if (translucent) flags = flags + 8
            if (present(two_sided)) then
                if (two_sided) flags = flags + 16
            end if
        end if
        if (present(colors)) flags = flags + 32
        node = vkf_mesh_create(int(nv, c_int), int(ni, c_int), real(rgba, c_float), c_null_ptr, c_null_ptr, &
                               int(flags, c_int))
        if (.not. c_associated(node)) return
        call mesh_arrays(node, nv, ni, v, n, uv, indices)
        v = real(xyz, c_float)
        if (lit) then
            n = real(normals, c_float)
        else
            n = 0.0_c_float
            n(3, :) = 1.0_c_float
        end if
        uv = 0.0_c_float
        indices = int(reshape(idx, [ni]) - 1, c_int32_t)
        if (present(colors)) then
            call c_f_pointer(vkf_mesh_array(node, 5_c_int), c, [4, nv])
            do i = 1, nv
                c(:, i) = real(colors(:, i) * rgba, c_float)
            end do
        end if
        if (lit .and. iand(flags, 8) /= 0 .and. present(sphere)) &
            call check_status(vkf_mesh_set_bound(node, real(sphere(1), c_double), real(sphere(2), c_double), &
                                                 real(sphere(3), c_double), real(sphere(4), c_double)), 'shapes')
    end function surface_node

    !> Points along a closed curve, as separate segments: pairs (p(i), p(i+1)).
    pure function closed_segments(p) result(s)
        real(wp), intent(in) :: p(:,:)                      !! (3, n) points of the curve
        real(wp) :: s(3, 2 * size(p, 2))
        integer :: i, n
        n = size(p, 2)
        do i = 1, n
            s(:, 2*i-1) = p(:, i)
            s(:, 2*i) = p(:, modulo(i, n) + 1)
        end do
    end function closed_segments

    !> Points along an open curve, as separate segments.
    pure function open_segments(p) result(s)
        real(wp), intent(in) :: p(:,:)                      !! (3, n) points of the curve
        real(wp) :: s(3, 2 * (size(p, 2) - 1))
        integer :: i
        do i = 1, size(p, 2) - 1
            s(:, 2*i-1) = p(:, i)
            s(:, 2*i) = p(:, i + 1)
        end do
    end function open_segments

    !> Append points b to a (allocatable, (3, n)).
    pure subroutine append(a, b)
        real(wp), allocatable, intent(inout) :: a(:,:)      !! (3, n) points (unallocated: none)
        real(wp), intent(in) :: b(:,:)                      !! (3, m) points to append
        real(wp), allocatable :: tmp(:,:)
        if (.not. allocated(a)) allocate(a(3, 0))
        allocate(tmp(3, size(a, 2) + size(b, 2)))
        tmp(:, 1:size(a, 2)) = a
        tmp(:, size(a, 2)+1:) = b
        call move_alloc(tmp, a)
    end subroutine append

    ! ==================================================================
    ! CoordinateAxes

    !> Coordinate axes named `name`, of color `color` ([r, g, b] or [r, g, b,
    !> a]), 1 long, with major ticks every 0.5 and minor every 0.25 (OpenFrames'
    !> defaults).
    function new_coordinate_axes(name, color) result(a)
        character(len=*), intent(in) :: name                !! the frame's name
        real(wp), intent(in), optional :: color(:)          !! [r, g, b] or [r, g, b, a] (0 to 1)
        type(CoordinateAxes) :: a
        a%ReferenceFrame = ReferenceFrame(name, color)
        a%image_file = ''
        call place_axes(a, a%length_)
    end function new_coordinate_axes

    !> The frame's own axes (arrows) at `length` along x, y and z, half as long.
    subroutine place_axes(self, length)
        class(ReferenceFrame), intent(inout) :: self
        real(wp), intent(in) :: length                      !! scene units
        call self%move_x_axis([length, 0.0_wp, 0.0_wp], 0.5_wp * length)
        call self%move_y_axis([0.0_wp, length, 0.0_wp], 0.5_wp * length)
        call self%move_z_axis([0.0_wp, 0.0_wp, length], 0.5_wp * length)
    end subroutine place_axes

    !> The axes' length (scene units, >= 0; default 1); the frame's own axes are
    !> moved to the lines' ends.
    subroutine set_axis_length(self, length)
        class(CoordinateAxes), intent(inout) :: self
        real(wp), intent(in) :: length                      !! scene units (< 0: ignored)
        if (length < 0.0_wp) return
        self%length_ = length
        call place_axes(self, length)
        call make_axes(self)
    end subroutine set_axis_length

    !> The axes' length (scene units).
    pure real(wp) function axis_length(self)
        class(CoordinateAxes), intent(in) :: self
        axis_length = self%length_
    end function axis_length

    !> Which axes to draw: an ior of X_AXIS, Y_AXIS, Z_AXIS.
    subroutine set_draw_axes(self, axes)
        class(CoordinateAxes), intent(inout) :: self
        integer, intent(in) :: axes                         !! an ior of X_AXIS, Y_AXIS, Z_AXIS
        self%drawn = axes
        call make_axes(self)
    end subroutine set_draw_axes

    !> Which axes are drawn: an ior of X_AXIS, Y_AXIS, Z_AXIS.
    pure integer function draw_axes(self)
        class(CoordinateAxes), intent(in) :: self
        draw_axes = self%drawn
    end function draw_axes

    !> The axis lines' width (window units; default 1).
    subroutine set_axis_width(self, width)
        class(CoordinateAxes), intent(inout) :: self
        real(wp), intent(in) :: width                       !! window units (<= 0: ignored)
        if (width <= 0.0_wp) return
        self%width_ = width
        call make_axes(self)
    end subroutine set_axis_width

    !> The axis lines' width (window units).
    pure real(wp) function axis_width(self)
        class(CoordinateAxes), intent(in) :: self
        axis_width = self%width_
    end function axis_width

    !> Major and minor tick spacings (>= 0; 0 for none).
    subroutine set_tick_spacing(self, major, minor)
        class(CoordinateAxes), intent(inout) :: self
        real(wp), intent(in) :: major, minor                !! scene units (< 0: unchanged)
        if (major >= 0.0_wp) self%major = major
        if (minor >= 0.0_wp) self%minor = minor
        call make_axes(self)
    end subroutine set_tick_spacing

    !> Major and minor tick spacings.
    pure subroutine tick_spacing(self, major, minor)
        class(CoordinateAxes), intent(in) :: self
        real(wp), intent(out) :: major, minor               !! scene units
        major = self%major
        minor = self%minor
    end subroutine tick_spacing

    !> Major and minor tick sizes (window units; default 8 and 4).
    subroutine set_tick_size(self, major, minor)
        class(CoordinateAxes), intent(inout) :: self
        real(wp), intent(in) :: major, minor                !! window units (<= 0: unchanged)
        if (major > 0.0_wp) self%major_size = major
        if (minor > 0.0_wp) self%minor_size = minor
        call make_axes(self)
    end subroutine set_tick_size

    !> Major and minor tick sizes (window units).
    pure subroutine tick_size(self, major, minor)
        class(CoordinateAxes), intent(in) :: self
        real(wp), intent(out) :: major, minor               !! window units
        major = self%major_size
        minor = self%minor_size
    end subroutine tick_size

    !> Tick shape: a MARKER_* shape (default MARKER_DISK; OpenFrames'
    !> setTickShader with its marker shaders).
    subroutine set_tick_shape(self, shape)
        class(CoordinateAxes), intent(inout) :: self
        integer, intent(in) :: shape                        !! a MARKER_* shape
        self%shape_ = shape
        call make_axes(self)
    end subroutine set_tick_shape

    !> Tick shape: a MARKER_* shape.
    pure integer function tick_shape(self)
        class(CoordinateAxes), intent(in) :: self
        tick_shape = self%shape_
    end function tick_shape

    !> Draw the ticks as an image (OpenFrames' setTickImage; found with
    !> find_data_file; '' for disks again). False if it cannot be read.
    logical function set_tick_image(self, file) result(ok)
        class(CoordinateAxes), intent(inout) :: self
        character(len=*), intent(in) :: file                !! image file ('' for disks)
        character(len=:), allocatable :: path
        type(c_ptr) :: image
        ok = .true.
        image = c_null_ptr
        if (len_trim(file) > 0) then
            path = find_data_file(file)
            if (len(path) > 0) image = vkf_image_read(path, len(path, kind=c_int))
            ok = c_associated(image)
            if (.not. ok) return
        end if
        if (c_associated(self%image)) call vkf_image_delete(self%image)
        self%image = image
        self%image_file = trim(file)
        self%shape_ = merge(MARKER_IMAGE, MARKER_DISK, c_associated(image))
        call make_axes(self)
    end function set_tick_image

    !> The tick image file ('' for none).
    function tick_image(self) result(file)
        class(CoordinateAxes), intent(in) :: self
        character(len=:), allocatable :: file
        file = self%image_file
    end function tick_image

    !> The axes' color (lines and ticks).
    subroutine axes_set_color(self, color)
        class(CoordinateAxes), intent(inout) :: self
        real(wp), intent(in) :: color(:)                    !! [r, g, b] or [r, g, b, a] (0 to 1)
        call self%ReferenceFrame%set_color(color)
        call make_axes(self)
    end subroutine axes_set_color

    !> What kind of frame this is: 'CoordinateAxes'.
    function axes_frame_info(self) result(info)
        class(CoordinateAxes), intent(in) :: self
        character(len=:), allocatable :: info
        info = 'CoordinateAxes'
    end function axes_frame_info

    !> The bounding sphere: the origin, radius the axis length.
    subroutine axes_bound(self, center, radius)
        class(CoordinateAxes), intent(in) :: self
        real(wp), intent(out) :: center(3), radius          !! the bounding sphere, frame coordinates (scene units)
        center = 0.0_wp
        radius = self%length_
    end subroutine axes_bound

    !> Draw the axes under `group` (the contents); with a window's `viewer`, later
    !> drawings are compiled.
    subroutine axes_build_contents(self, viewer, group)
        class(CoordinateAxes), intent(inout) :: self
        type(c_ptr), intent(in) :: viewer                   !! a window's viewer (null: headless)
        type(c_ptr), intent(in) :: group                    !! the frame's contents group
        self%nodes%contents = group
        call make_axes(self)
        self%nodes%shown = c_associated(viewer)
    end subroutine axes_build_contents

    !> Make the axes' lines and tick marks, in place of the old (nothing before
    !> build_contents).
    subroutine make_axes(self)
        class(CoordinateAxes), intent(inout) :: self
        real(wp), allocatable :: lines(:,:), major(:,:), minor(:,:)
        real(wp) :: e(3), mark, f
        type(c_ptr) :: group
        integer :: k, i
        if (.not. c_associated(self%nodes%contents)) return
        allocate(lines(3, 0), major(3, 0), minor(3, 0))
        do k = 1, 3
            if (iand(self%drawn, 2**(k-1)) == 0) cycle
            e = 0.0_wp
            e(k) = 1.0_wp
            call append(lines, reshape([0.0_wp, 0.0_wp, 0.0_wp, self%length_ * e], [3, 2]))
            if (self%major > 0.0_wp) then
                do i = 1, int(self%length_ / self%major + 1.0e-9_wp)
                    call append(major, reshape(i * self%major * e, [3, 1]))
                end do
            end if
            if (self%minor > 0.0_wp) then
                do i = 1, int(self%length_ / self%minor + 1.0e-9_wp)
                    mark = i * self%minor
                    if (self%major > 0.0_wp) then                   ! not where a major tick is
                        f = mark / self%major
                        if (abs(f - nint(f)) < 1.0e-6_wp) cycle
                    end if
                    call append(minor, reshape(mark * e, [3, 1]))
                end do
            end if
        end do
        group = vkf_switch_create()
        call add(group, line_node(lines, self%color(), self%width_, .true.))
        call add(group, marker_node(major, self%color(), self%major_size, self%shape_, self%length_, self%image))
        call add(group, marker_node(minor, self%color(), self%minor_size, self%shape_, self%length_, self%image))
        call self%nodes%replace(group)
    end subroutine make_axes

    ! ==================================================================
    ! LatLonGrid

    !> A latitude/longitude grid named `name`, of color `color`, on a unit
    !> sphere, every 30 degrees (OpenFrames' defaults).
    function new_lat_lon_grid(name, color) result(g)
        character(len=*), intent(in) :: name                !! the frame's name
        real(wp), intent(in), optional :: color(:)          !! [r, g, b] or [r, g, b, a] (0 to 1)
        type(LatLonGrid) :: g
        g%ReferenceFrame = ReferenceFrame(name, color)
        call place_grid_axes(g)
    end function new_lat_lon_grid

    !> The frame's own axes (arrows) at the ellipsoid's semi-axes, half as long.
    subroutine place_grid_axes(self)
        type(LatLonGrid), intent(inout) :: self
        call self%move_x_axis([self%radii_(1), 0.0_wp, 0.0_wp], 0.5_wp * self%radii_(1))
        call self%move_y_axis([0.0_wp, self%radii_(2), 0.0_wp], 0.5_wp * self%radii_(2))
        call self%move_z_axis([0.0_wp, 0.0_wp, self%radii_(3)], 0.5_wp * self%radii_(3))
    end subroutine place_grid_axes

    !> The ellipsoid's semi-axes (> 0; e.g. a Sphere's radii, a little larger
    !> to draw above its surface) and the line spacings (radians; 0 for no
    !> latitude or longitude lines), as OpenFrames' setParameters.
    subroutine grid_set_parameters(self, radius_x, radius_y, radius_z, lat_spacing, lon_spacing)
        class(LatLonGrid), intent(inout) :: self
        real(wp), intent(in) :: radius_x, radius_y, radius_z  !! scene units (<= 0: unchanged)
        real(wp), intent(in) :: lat_spacing, lon_spacing    !! radians (< 0: unchanged)
        if (radius_x > 0.0_wp) self%radii_(1) = radius_x
        if (radius_y > 0.0_wp) self%radii_(2) = radius_y
        if (radius_z > 0.0_wp) self%radii_(3) = radius_z
        if (lat_spacing >= 0.0_wp) self%lat_spacing = lat_spacing
        if (lon_spacing >= 0.0_wp) self%lon_spacing = lon_spacing
        call place_grid_axes(self)
        call make_grid(self)
    end subroutine grid_set_parameters

    !> The ellipsoid's semi-axes and the line spacings, as OpenFrames' getParameters.
    pure subroutine grid_parameters(self, radius_x, radius_y, radius_z, lat_spacing, lon_spacing)
        class(LatLonGrid), intent(in) :: self
        real(wp), intent(out) :: radius_x, radius_y, radius_z  !! scene units
        real(wp), intent(out) :: lat_spacing, lon_spacing   !! radians
        radius_x = self%radii_(1)
        radius_y = self%radii_(2)
        radius_z = self%radii_(3)
        lat_spacing = self%lat_spacing
        lon_spacing = self%lon_spacing
    end subroutine grid_parameters

    !> Widths of the grid lines and of the equator and prime meridian (window
    !> units; defaults 1 and 3).
    subroutine set_line_widths(self, width, main_width)
        class(LatLonGrid), intent(inout) :: self
        real(wp), intent(in) :: width, main_width           !! window units (<= 0: unchanged)
        if (width > 0.0_wp) self%width_ = width
        if (main_width > 0.0_wp) self%main_width = main_width
        call make_grid(self)
    end subroutine set_line_widths

    !> Widths of the grid lines and of the equator and prime meridian (window
    !> units).
    pure subroutine line_widths(self, width, main_width)
        class(LatLonGrid), intent(in) :: self
        real(wp), intent(out) :: width, main_width          !! window units
        width = self%width_
        main_width = self%main_width
    end subroutine line_widths

    !> The lines' color.
    subroutine grid_set_color(self, color)
        class(LatLonGrid), intent(inout) :: self
        real(wp), intent(in) :: color(:)                    !! [r, g, b] or [r, g, b, a] (0 to 1)
        call self%ReferenceFrame%set_color(color)
        call make_grid(self)
    end subroutine grid_set_color

    !> What kind of frame this is: 'LatLonGrid'.
    function grid_frame_info(self) result(info)
        class(LatLonGrid), intent(in) :: self
        character(len=:), allocatable :: info
        info = 'LatLonGrid'
    end function grid_frame_info

    !> The bounding sphere: the origin, radius the largest semi-axis.
    subroutine grid_bound(self, center, radius)
        class(LatLonGrid), intent(in) :: self
        real(wp), intent(out) :: center(3), radius          !! the bounding sphere, frame coordinates (scene units)
        center = 0.0_wp
        radius = maxval(self%radii_)
    end subroutine grid_bound

    !> Draw the grid under `group` (the contents); with a window's `viewer`, later
    !> drawings are compiled.
    subroutine grid_build_contents(self, viewer, group)
        class(LatLonGrid), intent(inout) :: self
        type(c_ptr), intent(in) :: viewer                   !! a window's viewer (null: headless)
        type(c_ptr), intent(in) :: group                    !! the frame's contents group
        self%nodes%contents = group
        call make_grid(self)
        self%nodes%shown = c_associated(viewer)
    end subroutine grid_build_contents

    !> The point at geodetic latitude lat and longitude lon on the ellipsoid
    !> (OpenFrames' LatLonGrid formula: the surface point whose normal has
    !> those angles).
    pure function surface_point(r, lat, lon) result(p)
        real(wp), intent(in) :: r(3)                        !! semi-axes (scene units)
        real(wp), intent(in) :: lat, lon                    !! radians
        real(wp) :: p(3), n(3), r2(3)
        n = [cos(lat) * cos(lon), cos(lat) * sin(lon), sin(lat)]
        r2 = r**2
        p = r2 * n / sqrt(dot_product(r2, n**2))
    end function surface_point

    !> Make the grid's lines (circles in 1 degree steps, meridians in 0.5 degree
    !> steps), in place of the old.
    subroutine make_grid(self)
        type(LatLonGrid), intent(inout) :: self
        integer, parameter :: n_around = 360, n_meridian = 360          ! 1 and 0.5 degree steps
        real(wp), allocatable :: grid(:,:), main(:,:)
        real(wp) :: circle(3, n_around), meridian(3, n_meridian + 1), lat, lon
        type(c_ptr) :: group
        integer :: i, k, sgn
        if (.not. c_associated(self%nodes%contents)) return
        allocate(grid(3, 0), main(3, 0))
        if (self%lat_spacing > 0.0_wp) then                    ! the equator, then pairs north and south
            k = 0
            do while (k * self%lat_spacing < 0.5_wp * pi - 1.0e-6_wp)
                do sgn = 1, merge(1, 2, k == 0)
                    lat = merge(1.0_wp, -1.0_wp, sgn == 1) * k * self%lat_spacing
                    do i = 1, n_around
                        circle(:, i) = surface_point(self%radii_, lat, 2.0_wp * pi * (i - 1) / n_around)
                    end do
                    if (k == 0) then
                        call append(main, closed_segments(circle))
                    else
                        call append(grid, closed_segments(circle))
                    end if
                end do
                k = k + 1
            end do
        end if
        if (self%lon_spacing > 0.0_wp) then                    ! the prime meridian, then the others
            k = 0
            do while (k * self%lon_spacing < 2.0_wp * pi - 1.0e-6_wp)
                lon = k * self%lon_spacing
                do i = 1, n_meridian + 1
                    meridian(:, i) = surface_point(self%radii_, -0.5_wp * pi + pi * (i - 1) / n_meridian, lon)
                end do
                if (k == 0) then
                    call append(main, open_segments(meridian))
                else
                    call append(grid, open_segments(meridian))
                end if
                k = k + 1
            end do
        end if
        group = vkf_switch_create()
        call add(group, line_node(grid, self%color(), self%width_, .true.))
        call add(group, line_node(main, self%color(), self%main_width, .true.))
        call self%nodes%replace(group)
    end subroutine make_grid

    ! ==================================================================
    ! RadialPlane

    !> A radial plane named `name` (decorations of color `color`): radius 1,
    !> circles every 0.5, radial lines every 30 degrees, light gray and
    !> transparent (OpenFrames' defaults).
    function new_radial_plane(name, color) result(p)
        character(len=*), intent(in) :: name                !! the frame's name
        real(wp), intent(in), optional :: color(:)          !! [r, g, b] or [r, g, b, a] (0 to 1)
        type(RadialPlane) :: p
        p%ReferenceFrame = ReferenceFrame(name, color)
        call place_axes(p, p%radius_)
    end function new_radial_plane

    !> The plane's radius (> 0), circle spacing and radial line spacing
    !> (radians; 0 for none), as OpenFrames' setParameters.
    subroutine plane_set_parameters(self, radius, radial_spacing, lon_spacing)
        class(RadialPlane), intent(inout) :: self
        real(wp), intent(in) :: radius                      !! scene units (<= 0: unchanged)
        real(wp), intent(in) :: radial_spacing              !! scene units (< 0: unchanged)
        real(wp), intent(in) :: lon_spacing                 !! radians (< 0: unchanged)
        if (radius > 0.0_wp) self%radius_ = radius
        if (radial_spacing >= 0.0_wp) self%radial_spacing = radial_spacing
        if (lon_spacing >= 0.0_wp) self%lon_spacing = lon_spacing
        call place_axes(self, self%radius_)
        call make_plane(self)
    end subroutine plane_set_parameters

    !> The plane's radius, circle spacing and radial line spacing, as OpenFrames'
    !> getParameters.
    pure subroutine plane_parameters(self, radius, radial_spacing, lon_spacing)
        class(RadialPlane), intent(in) :: self
        real(wp), intent(out) :: radius, radial_spacing     !! scene units
        real(wp), intent(out) :: lon_spacing                !! radians
        radius = self%radius_
        radial_spacing = self%radial_spacing
        lon_spacing = self%lon_spacing
    end subroutine plane_parameters

    !> The plane's color ([r, g, b] or [r, g, b, a]; default 0.8 gray, 0.2 opaque).
    subroutine set_plane_color(self, color)
        class(RadialPlane), intent(inout) :: self
        real(wp), intent(in) :: color(:)                    !! [r, g, b] or [r, g, b, a] (0 to 1)
        call set_rgba(self%plane_color_, color)
        call make_plane(self)
    end subroutine set_plane_color

    !> The plane's color [r, g, b, a].
    pure function plane_color(self) result(c)
        class(RadialPlane), intent(in) :: self
        real(wp) :: c(4)
        c = self%plane_color_
    end function plane_color

    !> The lines' color (default white, 0.2 opaque).
    subroutine plane_set_line_color(self, color)
        class(RadialPlane), intent(inout) :: self
        real(wp), intent(in) :: color(:)                    !! [r, g, b] or [r, g, b, a] (0 to 1)
        call set_rgba(self%line_color_, color)
        call make_plane(self)
    end subroutine plane_set_line_color

    !> The lines' color [r, g, b, a].
    pure function plane_line_color(self) result(c)
        class(RadialPlane), intent(in) :: self
        real(wp) :: c(4)
        c = self%line_color_
    end function plane_line_color

    !> What kind of frame this is: 'RadialPlane'.
    function plane_frame_info(self) result(info)
        class(RadialPlane), intent(in) :: self
        character(len=:), allocatable :: info
        info = 'RadialPlane'
    end function plane_frame_info

    !> The bounding sphere: the origin, radius the plane's radius.
    subroutine plane_bound(self, center, radius)
        class(RadialPlane), intent(in) :: self
        real(wp), intent(out) :: center(3), radius          !! the bounding sphere, frame coordinates (scene units)
        center = 0.0_wp
        radius = self%radius_
    end subroutine plane_bound

    !> Draw the plane under `group` (the contents); with a window's `viewer`, later
    !> drawings are compiled.
    subroutine plane_build_contents(self, viewer, group)
        class(RadialPlane), intent(inout) :: self
        type(c_ptr), intent(in) :: viewer                   !! a window's viewer (null: headless)
        type(c_ptr), intent(in) :: group                    !! the frame's contents group
        self%nodes%contents = group
        call make_plane(self)
        self%nodes%shown = c_associated(viewer)
    end subroutine plane_build_contents

    !> Make the disk, its circles and its radial lines (1 degree steps), in place of
    !> the old.
    subroutine make_plane(self)
        type(RadialPlane), intent(inout) :: self
        integer, parameter :: n = 360                          ! points around a circle (1 degree)
        real(wp) :: disk(3, n + 1), circle(3, n), a, r
        integer :: idx(3, n), i, k
        real(wp), allocatable :: lines(:,:), main(:,:)
        type(c_ptr) :: group
        if (.not. c_associated(self%nodes%contents)) return
        ! the disk: a fan around the center
        disk(:, n + 1) = 0.0_wp
        do i = 1, n
            a = 2.0_wp * pi * (i - 1) / n
            disk(:, i) = self%radius_ * [cos(a), sin(a), 0.0_wp]
            idx(:, i) = [n + 1, i, modulo(i, n) + 1]
        end do
        allocate(lines(3, 0), main(3, 0))
        ! circles: the edge, and every radial_spacing from the center
        if (self%radial_spacing > 0.0_wp) then
            call append(lines, closed_segments(disk(:, 1:n)))
            k = 1
            do while (k * self%radial_spacing < self%radius_ * (1.0_wp - 1.0e-6_wp))
                r = k * self%radial_spacing
                circle = r / self%radius_ * disk(:, 1:n)
                call append(lines, closed_segments(circle))
                k = k + 1
            end do
        end if
        ! radial lines: along +x (thicker), then every lon_spacing
        if (self%lon_spacing > 0.0_wp) then
            k = 0
            do while (k * self%lon_spacing < 2.0_wp * pi - 1.0e-6_wp)
                a = k * self%lon_spacing
                if (k == 0) then
                    call append(main, reshape([0.0_wp, 0.0_wp, 0.0_wp, self%radius_ * [cos(a), sin(a), 0.0_wp]], [3, 2]))
                else
                    call append(lines, reshape([0.0_wp, 0.0_wp, 0.0_wp, self%radius_ * [cos(a), sin(a), 0.0_wp]], [3, 2]))
                end if
                k = k + 1
            end do
        end if
        group = vkf_switch_create()
        call add(group, surface_node(disk, idx, self%plane_color_))
        call add(group, line_node(lines, self%line_color_, 1.0_wp, .true.))
        call add(group, line_node(main, self%line_color_, 3.0_wp, .true.))
        call self%nodes%replace(group)
    end subroutine make_plane

    ! ==================================================================
    ! PolyhedralCone

    !> A polyhedral cone named `name`, of color `color` (its sides and lines):
    !> 3 vertices at clock angles 0, 120 and 240 degrees, cone angles 45
    !> degrees, 1 long (OpenFrames' defaults).
    function new_polyhedral_cone(name, color) result(c)
        character(len=*), intent(in) :: name                !! the frame's name
        real(wp), intent(in), optional :: color(:)          !! [r, g, b] or [r, g, b, a] (0 to 1)
        type(PolyhedralCone) :: c
        c%ReferenceFrame = ReferenceFrame(name, color)
        c%cone_color_ = c%color()
        c%line_color_ = c%color()
        c%clock = [0.0_wp, 120.0_wp, 240.0_wp] * pi / 180.0_wp
        c%cone = [45.0_wp, 45.0_wp, 45.0_wp] * pi / 180.0_wp
    end function new_polyhedral_cone

    !> The vertices: clock angles (from +x toward +y) and cone angles (from the
    !> boresight, -z; below 90 degrees), radians, the same number of each (>= 3).
    subroutine set_vertex_angles(self, clock, cone)
        class(PolyhedralCone), intent(inout) :: self
        real(wp), intent(in) :: clock(:), cone(:)           !! radians
        if (size(clock) /= size(cone) .or. size(clock) < 3) return
        self%clock = clock
        self%cone = cone
        call make_cone(self)
    end subroutine set_vertex_angles

    !> The vertices' clock angles (radians).
    pure function clock_angles(self) result(a)
        class(PolyhedralCone), intent(in) :: self
        real(wp), allocatable :: a(:)
        a = self%clock
    end function clock_angles

    !> The vertices' cone angles (radians).
    pure function cone_angles(self) result(a)
        class(PolyhedralCone), intent(in) :: self
        real(wp), allocatable :: a(:)
        a = self%cone
    end function cone_angles

    !> The cone's length along the boresight (> 0; default 1).
    subroutine set_cone_length(self, length)
        class(PolyhedralCone), intent(inout) :: self
        real(wp), intent(in) :: length                      !! scene units (<= 0: ignored)
        if (length <= 0.0_wp) return
        self%length_ = length
        call make_cone(self)
    end subroutine set_cone_length

    !> The cone's length along the boresight (scene units).
    pure real(wp) function cone_length(self)
        class(PolyhedralCone), intent(in) :: self
        cone_length = self%length_
    end function cone_length

    !> The sides' color ([r, g, b] or [r, g, b, a]; their alpha makes them
    !> transparent).
    subroutine set_cone_color(self, color)
        class(PolyhedralCone), intent(inout) :: self
        real(wp), intent(in) :: color(:)                    !! [r, g, b] or [r, g, b, a] (0 to 1)
        call set_rgba(self%cone_color_, color)
        call make_cone(self)
    end subroutine set_cone_color

    !> The sides' color [r, g, b, a].
    pure function cone_color(self) result(c)
        class(PolyhedralCone), intent(in) :: self
        real(wp) :: c(4)
        c = self%cone_color_
    end function cone_color

    !> The edges' and base outline's color.
    subroutine cone_set_line_color(self, color)
        class(PolyhedralCone), intent(inout) :: self
        real(wp), intent(in) :: color(:)                    !! [r, g, b] or [r, g, b, a] (0 to 1)
        call set_rgba(self%line_color_, color)
        call make_cone(self)
    end subroutine cone_set_line_color

    !> The edges' and base outline's color [r, g, b, a].
    pure function cone_line_color(self) result(c)
        class(PolyhedralCone), intent(in) :: self
        real(wp) :: c(4)
        c = self%line_color_
    end function cone_line_color

    !> What to draw: an ior of CONE_SIDES, CONE_EDGES, CONE_BASE_OUTLINE.
    subroutine set_draw_mode(self, mode)
        class(PolyhedralCone), intent(inout) :: self
        integer, intent(in) :: mode                         !! an ior of CONE_SIDES, CONE_EDGES, CONE_BASE_OUTLINE
        self%mode = mode
        call make_cone(self)
    end subroutine set_draw_mode

    !> What is drawn: an ior of CONE_SIDES, CONE_EDGES, CONE_BASE_OUTLINE.
    pure integer function draw_mode(self)
        class(PolyhedralCone), intent(in) :: self
        draw_mode = self%mode
    end function draw_mode

    !> Place the cone's apex at `apex` (in the parent frame) with its boresight
    !> along `dir` and its +y axis toward `up` (OpenFrames' makeConeLookAt).
    subroutine make_cone_look_at(self, apex, dir, up)
        class(PolyhedralCone), intent(inout) :: self
        real(wp), intent(in) :: apex(3)                     !! parent frame coordinates (scene units)
        real(wp), intent(in) :: dir(3)                      !! boresight direction (any length > 0)
        real(wp), intent(in) :: up(3)                       !! toward the cone's +y axis (any length)
        real(wp) :: f(3), s(3), u(3), m(3,3)
        if (norm2(dir) <= 0.0_wp) return
        f = dir / norm2(dir)
        s = cross(f, up)
        if (norm2(s) <= 1.0e-12_wp) s = cross(f, [1.0_wp, 0.0_wp, 0.0_wp])
        if (norm2(s) <= 1.0e-12_wp) s = cross(f, [0.0_wp, 1.0_wp, 0.0_wp])
        s = s / norm2(s)
        u = cross(s, f)
        m(:, 1) = s                                             ! the cone's axes in the parent frame
        m(:, 2) = u
        m(:, 3) = -f
        call self%set_position(apex)
        call self%set_attitude(quat_from_matrix(m))
    end subroutine make_cone_look_at

    !> Whether `point` (in the cone's frame) is inside the cone, between
    !> distances `min_distance` and `max_distance` along the boresight
    !> (defaults 0 and unlimited): inside the polygon of the vertices' directions
    !> (OpenFrames' isVisible; its PolyhedralCone always answers false).
    logical function cone_is_visible(self, point, min_distance, max_distance) result(visible)
        class(PolyhedralCone), intent(in) :: self
        real(wp), intent(in) :: point(3)                    !! in the cone's frame (scene units)
        real(wp), intent(in), optional :: min_distance, max_distance  !! along the boresight (scene units)
        real(wp) :: q(2), a(2), b(2)
        integer :: i, n, winding
        visible = in_range(point, min_distance, max_distance)
        if (.not. visible .or. norm2(point) == 0.0_wp) return
        q = point(1:2) / (-point(3))                            ! on the plane z = -1
        n = size(self%clock)
        winding = 0                                             ! winding number of the polygon around q
        do i = 1, n
            a = tan(self%cone(i)) * [cos(self%clock(i)), sin(self%clock(i))]
            b = tan(self%cone(modulo(i, n) + 1)) * [cos(self%clock(modulo(i, n) + 1)), sin(self%clock(modulo(i, n) + 1))]
            if (a(2) <= q(2)) then
                if (b(2) > q(2) .and. side(a, b, q) > 0.0_wp) winding = winding + 1
            else
                if (b(2) <= q(2) .and. side(a, b, q) < 0.0_wp) winding = winding - 1
            end if
        end do
        visible = winding /= 0
    end function cone_is_visible

    !> Along the boresight (-z) between the distances (a point ahead of the apex).
    pure logical function in_range(point, min_distance, max_distance)
        real(wp), intent(in) :: point(3)                    !! in the cone's frame (scene units)
        real(wp), intent(in), optional :: min_distance, max_distance  !! along the boresight (scene units)
        real(wp) :: d
        d = -point(3)
        in_range = d >= 0.0_wp
        if (present(min_distance)) in_range = in_range .and. d >= min_distance
        if (present(max_distance)) in_range = in_range .and. d <= max_distance
    end function in_range

    !> Which side of the line from a to b q is on: > 0 left, < 0 right, 0 on it
    !> (twice the signed area of the triangle a, b, q).
    pure real(wp) function side(a, b, q)
        real(wp), intent(in) :: a(2), b(2), q(2)            !! the line's two points, and the point tested
        side = (b(1) - a(1)) * (q(2) - a(2)) - (q(1) - a(1)) * (b(2) - a(2))
    end function side

    !> What kind of frame this is: 'PolyhedralCone'.
    function cone_frame_info(self) result(info)
        class(PolyhedralCone), intent(in) :: self
        character(len=:), allocatable :: info
        info = 'PolyhedralCone'
    end function cone_frame_info

    !> The bound of the apex and the base.
    subroutine cone_bound(self, center, radius)
        class(PolyhedralCone), intent(in) :: self
        real(wp), intent(out) :: center(3), radius          !! the bounding sphere, frame coordinates (scene units)
        real(wp) :: r
        r = self%length_ * maxval(tan(self%cone))
        center = [0.0_wp, 0.0_wp, -0.5_wp * self%length_]
        radius = sqrt((0.5_wp * self%length_)**2 + r**2)
    end subroutine cone_bound

    !> Draw the cone under `group` (the contents); with a window's `viewer`, later
    !> drawings are compiled.
    subroutine cone_build_contents(self, viewer, group)
        class(PolyhedralCone), intent(inout) :: self
        type(c_ptr), intent(in) :: viewer                   !! a window's viewer (null: headless)
        type(c_ptr), intent(in) :: group                    !! the frame's contents group
        self%nodes%contents = group
        call make_cone(self)
        self%nodes%shown = c_associated(viewer)
    end subroutine cone_build_contents

    !> Make the sides, edges and base outline (those the draw mode says), in place
    !> of the old.
    subroutine make_cone(self)
        class(PolyhedralCone), intent(inout) :: self
        real(wp), allocatable :: v(:,:), edges(:,:)
        integer, allocatable :: idx(:,:)
        type(c_ptr) :: group
        integer :: i, n
        if (.not. c_associated(self%nodes%contents)) return
        n = size(self%clock)
        allocate(v(3, n + 1), idx(3, n), edges(3, 2 * n))
        v(:, n + 1) = 0.0_wp                                    ! the apex
        do i = 1, n
            v(:, i) = self%length_ * [tan(self%cone(i)) * cos(self%clock(i)), tan(self%cone(i)) * sin(self%clock(i)), -1.0_wp]
            idx(:, i) = [n + 1, i, modulo(i, n) + 1]
            edges(:, 2*i-1) = 0.0_wp
            edges(:, 2*i) = v(:, i)
        end do
        group = vkf_switch_create()
        if (iand(self%mode, CONE_SIDES) /= 0) call add(group, surface_node(v, idx, self%cone_color_))
        if (iand(self%mode, CONE_EDGES) /= 0) call add(group, line_node(edges, self%line_color_, self%width_, .true.))
        if (iand(self%mode, CONE_BASE_OUTLINE) /= 0) &
            call add(group, line_node(closed_segments(v(:, 1:n)), self%line_color_, self%width_, .true.))
        call self%nodes%replace(group)
    end subroutine make_cone

    ! ==================================================================
    ! EllipticCone and RectangularCone

    !> An elliptic cone named `name`, of color `color` (sides and base outline drawn).
    function new_elliptic_cone(name, color) result(c)
        character(len=*), intent(in) :: name                !! the frame's name
        real(wp), intent(in), optional :: color(:)          !! [r, g, b] or [r, g, b, a] (0 to 1)
        type(EllipticCone) :: c
        c%PolyhedralCone = PolyhedralCone(name, color)
        c%mode = ior(CONE_SIDES, CONE_BASE_OUTLINE)
        call elliptic_vertices(c)
    end function new_elliptic_cone

    !> The half-angles along x (a) and y (b), radians (0 to 90 degrees).
    subroutine elliptic_set_primary_angles(self, a, b)
        class(EllipticCone), intent(inout) :: self
        real(wp), intent(in) :: a, b                        !! radians
        self%a = a
        self%b = b
        call elliptic_vertices(self)
        call make_cone(self)
    end subroutine elliptic_set_primary_angles

    !> The half-angles along x (a) and y (b), radians.
    pure subroutine elliptic_primary_angles(self, a, b)
        class(EllipticCone), intent(in) :: self
        real(wp), intent(out) :: a, b                       !! radians
        a = self%a
        b = self%b
    end subroutine elliptic_primary_angles

    !> 120 vertices (every 3 degrees), cone angle a b / sqrt(b^2 cos^2 + a^2 sin^2)
    !> (OpenFrames' createEllipticCone).
    subroutine elliptic_vertices(self)
        type(EllipticCone), intent(inout) :: self
        integer, parameter :: n = 120
        integer :: i
        real(wp) :: t
        if (allocated(self%clock)) deallocate(self%clock, self%cone)
        allocate(self%clock(n), self%cone(n))
        do i = 1, n
            t = 2.0_wp * pi * (i - 1) / n
            self%clock(i) = t
            self%cone(i) = self%a * self%b / sqrt(self%b**2 * cos(t)**2 + self%a**2 * sin(t)**2)
        end do
    end subroutine elliptic_vertices

    !> Inside the elliptic cone (exactly; OpenFrames' EllipticCone::isVisible).
    logical function elliptic_is_visible(self, point, min_distance, max_distance) result(visible)
        class(EllipticCone), intent(in) :: self
        real(wp), intent(in) :: point(3)                    !! in the cone's frame (scene units)
        real(wp), intent(in), optional :: min_distance, max_distance  !! along the boresight (scene units)
        real(wp) :: p(3), cone, clock
        visible = in_range(point, min_distance, max_distance)
        if (.not. visible .or. norm2(point) == 0.0_wp) return
        p = point / (-point(3))
        cone = atan(norm2(p(1:2)))
        clock = atan2(p(2), p(1))
        visible = cone <= self%a * self%b / sqrt(self%b**2 * cos(clock)**2 + self%a**2 * sin(clock)**2)
    end function elliptic_is_visible

    !> What kind of frame this is: 'EllipticCone'.
    function elliptic_frame_info(self) result(info)
        class(EllipticCone), intent(in) :: self
        character(len=:), allocatable :: info
        info = 'EllipticCone'
    end function elliptic_frame_info

    !> A rectangular cone named `name`, of color `color`.
    function new_rectangular_cone(name, color) result(c)
        character(len=*), intent(in) :: name                !! the frame's name
        real(wp), intent(in), optional :: color(:)          !! [r, g, b] or [r, g, b, a] (0 to 1)
        type(RectangularCone) :: c
        c%PolyhedralCone = PolyhedralCone(name, color)
        call rect_vertices(c)
    end function new_rectangular_cone

    !> The half-angles along x and y, radians (0 to 90 degrees).
    subroutine rect_set_primary_angles(self, x, y)
        class(RectangularCone), intent(inout) :: self
        real(wp), intent(in) :: x, y                        !! radians
        self%x = x
        self%y = y
        call rect_vertices(self)
        call make_cone(self)
    end subroutine rect_set_primary_angles

    !> The half-angles along x and y, radians.
    pure subroutine rect_primary_angles(self, x, y)
        class(RectangularCone), intent(in) :: self
        real(wp), intent(out) :: x, y                       !! radians
        x = self%x
        y = self%y
    end subroutine rect_primary_angles

    !> The 4 corners (OpenFrames' createRectangularCone).
    subroutine rect_vertices(self)
        type(RectangularCone), intent(inout) :: self
        real(wp) :: xl, yl, cone, clock
        xl = tan(self%x)
        yl = tan(self%y)
        cone = atan(sqrt(xl**2 + yl**2))
        clock = atan2(yl, xl)
        self%clock = [clock, pi - clock, pi + clock, -clock]
        self%cone = [cone, cone, cone, cone]
    end subroutine rect_vertices

    !> What kind of frame this is: 'RectangularCone'.
    function rect_frame_info(self) result(info)
        class(RectangularCone), intent(in) :: self
        character(len=:), allocatable :: info
        info = 'RectangularCone'
    end function rect_frame_info

    ! ==================================================================
    ! Polyhedron

    !> A polyhedron named `name`, of color `color` (its faces, edges and
    !> vertices), with no faces yet (set_faces, set_triangles or set_grid).
    function new_polyhedron(name, color) result(p)
        character(len=*), intent(in) :: name                !! the frame's name
        real(wp), intent(in), optional :: color(:)          !! [r, g, b] or [r, g, b, a] (0 to 1)
        type(Polyhedron) :: p
        p%ReferenceFrame = ReferenceFrame(name, color)
        p%face_color_ = p%color()
        p%line_color_ = p%color()
        p%point_color_ = p%color()
        allocate(p%v(3, 0), p%tri(3, 0), p%edge(2, 0))
    end function new_polyhedron

    !> Faces given as polygons: `faces` the vertex numbers (1-based) of each
    !> polygon in turn, `sizes` how many each has (>= 3; convex, or at least
    !> star-shaped from its first vertex: drawn as a fan of triangles from
    !> it); the edges are the polygons' sides. E.g. a cube: 8 vertices, 6
    !> faces of 4. A face's front: its vertices counterclockwise.
    subroutine set_faces(self, vertices, faces, sizes)
        class(Polyhedron), intent(inout) :: self
        real(wp), intent(in) :: vertices(:,:)               !! (3, n), frame coordinates
        integer, intent(in) :: faces(:)                     !! (sum(sizes)) vertex numbers
        integer, intent(in) :: sizes(:)                     !! each face's number of vertices
        integer :: f, k, j, nt, ne, first
        integer, allocatable :: sides(:,:)
        self%v = vertices
        nt = sum(max(sizes - 2, 0))
        if (allocated(self%tri)) deallocate(self%tri)
        allocate(self%tri(3, nt), sides(2, sum(sizes)))
        nt = 0
        ne = 0
        first = 0
        do f = 1, size(sizes)
            do k = 2, sizes(f) - 1                          ! a fan from the first vertex
                nt = nt + 1
                self%tri(:, nt) = [faces(first + 1), faces(first + k), faces(first + k + 1)]
            end do
            do k = 1, sizes(f)
                j = modulo(k, sizes(f)) + 1
                ne = ne + 1
                sides(:, ne) = [faces(first + k), faces(first + j)]
            end do
            first = first + sizes(f)
        end do
        self%edge = unique_edges(sides(:, 1:ne), size(vertices, 2))
        call drop_colors(self)
        call make_polyhedron(self)
    end subroutine set_faces

    !> Faces given as triangles (3, m) of vertex numbers (1-based); the edges
    !> are their sides (each once).
    subroutine set_triangles(self, vertices, triangles)
        class(Polyhedron), intent(inout) :: self
        real(wp), intent(in) :: vertices(:,:)               !! (3, n), frame coordinates
        integer, intent(in) :: triangles(:,:)               !! (3, m)
        self%v = vertices
        self%tri = triangles
        self%edge = unique_edges(reshape([triangles(1:2, :), triangles(2:3, :), triangles(3:1:-2, :)], &
                                         [2, 3 * size(triangles, 2)]), size(vertices, 2))
        call drop_colors(self)
        call make_polyhedron(self)
    end subroutine set_triangles

    !> Faces given as a grid of points (3, nu, nv): quadrilaterals between
    !> neighbours (two triangles each), joined around if `closed_u` (the last
    !> point i to the first, for each j; e.g. a tube's circumference) or
    !> `closed_v`. The edges: the grid's lines, those of constant i (along j)
    !> every lines(1) and of constant j (along i) every lines(2) (default:
    !> every one; 0 for none).
    subroutine set_grid(self, points, closed_u, closed_v, lines)
        class(Polyhedron), intent(inout) :: self
        real(wp), intent(in) :: points(:,:,:)               !! (3, nu, nv), frame coordinates
        logical, intent(in), optional :: closed_u, closed_v !! joined around in i, in j (default false)
        integer, intent(in), optional :: lines(2)           !! edges every lines(1) in i, lines(2) in j (default 1; 0: none)
        integer :: nu, nv, mu, mv, i, j, i2, j2, nt, ne, every(2)
        logical :: cu, cv
        nu = size(points, 2)
        nv = size(points, 3)
        cu = .false.
        cv = .false.
        if (present(closed_u)) cu = closed_u .and. nu > 2
        if (present(closed_v)) cv = closed_v .and. nv > 2
        every = 1
        if (present(lines)) every = max(lines, 0)
        self%v = reshape(points, [3, nu * nv])
        mu = merge(nu, nu - 1, cu)                          ! quadrilaterals along i, j
        mv = merge(nv, nv - 1, cv)
        if (allocated(self%tri)) deallocate(self%tri)
        allocate(self%tri(3, 2 * max(mu, 0) * max(mv, 0)))
        nt = 0
        do j = 1, mv
            j2 = modulo(j, nv) + 1
            do i = 1, mu
                i2 = modulo(i, nu) + 1
                self%tri(:, nt + 1) = [at(i, j), at(i2, j), at(i2, j2)]
                self%tri(:, nt + 2) = [at(i, j), at(i2, j2), at(i, j2)]
                nt = nt + 2
            end do
        end do
        if (allocated(self%edge)) deallocate(self%edge)
        allocate(self%edge(2, nu * nv * 2))
        ne = 0
        if (every(1) > 0) then                              ! lines of constant i, along j
            do i = 1, nu
                if (.not. shown_line(i, nu, every(1), cu)) cycle
                do j = 1, mv
                    ne = ne + 1
                    self%edge(:, ne) = [at(i, j), at(i, modulo(j, nv) + 1)]
                end do
            end do
        end if
        if (every(2) > 0) then                              ! lines of constant j, along i
            do j = 1, nv
                if (.not. shown_line(j, nv, every(2), cv)) cycle
                do i = 1, mu
                    ne = ne + 1
                    self%edge(:, ne) = [at(i, j), at(modulo(i, nu) + 1, j)]
                end do
            end do
        end if
        self%edge = self%edge(:, 1:ne)
        call drop_colors(self)
        call make_polyhedron(self)

    contains

        pure integer function at(i_, j_)
            integer, intent(in) :: i_, j_                       !! grid indices
            at = i_ + (j_ - 1) * nu
        end function at

        !> Line k of n drawn, every `step` (and the last, if the grid is open).
        pure logical function shown_line(k, n, step, closed)
            integer, intent(in) :: k, n, step                   !! the line, number of lines, and spacing
            logical, intent(in) :: closed                       !! the grid is joined around (else the last line drawn too)
            shown_line = modulo(k - 1, step) == 0 .or. (k == n .and. .not. closed)
        end function shown_line

    end subroutine set_grid

    !> The edges drawn (2, k) (vertex numbers), in place of the faces' (after
    !> set_faces, set_triangles or set_grid).
    subroutine set_edges(self, edges)
        class(Polyhedron), intent(inout) :: self
        integer, intent(in) :: edges(:,:)                   !! (2, k) vertex numbers (1-based)
        self%edge = edges
        call make_polyhedron(self)
    end subroutine set_edges

    !> Move the vertices (3, n), keeping the faces and edges (as many as
    !> before; e.g. a surface changing in time).
    subroutine set_vertices(self, vertices)
        class(Polyhedron), intent(inout) :: self
        real(wp), intent(in) :: vertices(:,:)               !! (3, n), frame coordinates (as many as before)
        if (size(vertices, 2) /= size(self%v, 2)) then
            write(*, '(a)') 'vkframes: Polyhedron '//self%name()//': set_vertices needs as many vertices as before'
            return
        end if
        self%v = vertices
        call make_polyhedron(self)
    end subroutine set_vertices

    !> The number of vertices.
    pure integer function num_vertices(self)
        class(Polyhedron), intent(in) :: self
        num_vertices = size(self%v, 2)
    end function num_vertices

    !> The number of triangles drawn (polygons: as fans).
    pure integer function num_triangles(self)
        class(Polyhedron), intent(in) :: self
        num_triangles = size(self%tri, 2)
    end function num_triangles

    !> The number of edges.
    pure integer function num_edges(self)
        class(Polyhedron), intent(in) :: self
        num_edges = size(self%edge, 2)
    end function num_edges

    !> A color per vertex (4, n) or (3, n) ([r, g, b(, a)], 0 to 1), for the
    !> faces, blended across them (in place of the face color; the opacity
    !> still applies). Dropped when the faces are set again.
    subroutine set_vertex_colors(self, colors)
        class(Polyhedron), intent(inout) :: self
        real(wp), intent(in) :: colors(:,:)                 !! (4, n) or (3, n): [r, g, b(, a)] (0 to 1)
        integer :: i
        if (size(colors, 2) /= size(self%v, 2) .or. size(colors, 1) < 3) then
            write(*, '(a)') 'vkframes: Polyhedron '//self%name()//': set_vertex_colors needs a color per vertex'
            return
        end if
        if (allocated(self%vertex_colors_)) deallocate(self%vertex_colors_)
        allocate(self%vertex_colors_(4, size(colors, 2)))
        do i = 1, size(colors, 2)
            call set_rgba(self%vertex_colors_(:, i), colors(:, i))
        end do
        call make_polyhedron(self)
    end subroutine set_vertex_colors

    !> The face color again, for every vertex.
    subroutine remove_vertex_colors(self)
        class(Polyhedron), intent(inout) :: self
        call drop_colors(self)
        call make_polyhedron(self)
    end subroutine remove_vertex_colors

    !> Whether the faces have a color per vertex.
    pure logical function has_vertex_colors(self)
        class(Polyhedron), intent(in) :: self
        has_vertex_colors = allocated(self%vertex_colors_)
    end function has_vertex_colors

    !> The faces' color (the alpha times the opacity).
    subroutine set_face_color(self, color)
        class(Polyhedron), intent(inout) :: self
        real(wp), intent(in) :: color(:)                    !! [r, g, b] or [r, g, b, a] (0 to 1)
        call set_rgba(self%face_color_, color)
        call make_polyhedron(self)
    end subroutine set_face_color

    !> The faces' color [r, g, b, a].
    pure function face_color(self) result(c)
        class(Polyhedron), intent(in) :: self
        real(wp) :: c(4)
        c = self%face_color_
    end function face_color

    !> The edges' color.
    subroutine poly_set_line_color(self, color)
        class(Polyhedron), intent(inout) :: self
        real(wp), intent(in) :: color(:)                    !! [r, g, b] or [r, g, b, a] (0 to 1)
        call set_rgba(self%line_color_, color)
        call make_polyhedron(self)
    end subroutine poly_set_line_color

    !> The edges' color [r, g, b, a].
    pure function poly_line_color(self) result(c)
        class(Polyhedron), intent(in) :: self
        real(wp) :: c(4)
        c = self%line_color_
    end function poly_line_color

    !> The vertices' color (POLY_VERTICES).
    subroutine set_point_color(self, color)
        class(Polyhedron), intent(inout) :: self
        real(wp), intent(in) :: color(:)                    !! [r, g, b] or [r, g, b, a] (0 to 1)
        call set_rgba(self%point_color_, color)
        call make_polyhedron(self)
    end subroutine set_point_color

    !> The vertices' color [r, g, b, a].
    pure function point_color(self) result(c)
        class(Polyhedron), intent(in) :: self
        real(wp) :: c(4)
        c = self%point_color_
    end function point_color

    !> The faces' opacity (0 to 1; default 1, opaque): below 1 translucent.
    subroutine poly_set_opacity(self, opacity)
        class(Polyhedron), intent(inout) :: self
        real(wp), intent(in) :: opacity                     !! 0 to 1 (1: opaque)
        self%opacity_ = max(0.0_wp, min(1.0_wp, opacity))
        call make_polyhedron(self)
    end subroutine poly_set_opacity

    !> The faces' opacity.
    pure real(wp) function poly_opacity(self)
        class(Polyhedron), intent(in) :: self
        poly_opacity = self%opacity_
    end function poly_opacity

    !> The edges' width (window units; default 1).
    subroutine poly_set_line_width(self, width)
        class(Polyhedron), intent(inout) :: self
        real(wp), intent(in) :: width                       !! window units
        self%width_ = max(0.0_wp, width)
        call make_polyhedron(self)
    end subroutine poly_set_line_width

    !> The edges' width (window units).
    pure real(wp) function poly_line_width(self)
        class(Polyhedron), intent(in) :: self
        poly_line_width = self%width_
    end function poly_line_width

    !> The vertices' size (window units; default 6).
    subroutine set_point_size(self, size_)
        class(Polyhedron), intent(inout) :: self
        real(wp), intent(in) :: size_                       !! window units
        self%point_size_ = max(0.0_wp, size_)
        call make_polyhedron(self)
    end subroutine set_point_size

    !> The vertices' size (window units).
    pure real(wp) function point_size(self)
        class(Polyhedron), intent(in) :: self
        point_size = self%point_size_
    end function point_size

    !> What to draw: an ior of POLY_FACES, POLY_EDGES, POLY_VERTICES (default
    !> faces and edges).
    subroutine poly_set_draw_mode(self, mode)
        class(Polyhedron), intent(inout) :: self
        integer, intent(in) :: mode                         !! an ior of POLY_FACES, POLY_EDGES, POLY_VERTICES
        self%mode = mode
        call make_polyhedron(self)
    end subroutine poly_set_draw_mode

    !> What is drawn: an ior of POLY_FACES, POLY_EDGES, POLY_VERTICES.
    pure integer function poly_draw_mode(self)
        class(Polyhedron), intent(in) :: self
        poly_draw_mode = self%mode
    end function poly_draw_mode

    !> Faces lit by the scene's lights (the default), or unlit (their color
    !> as is, transparent as a cone's sides: writing no depth).
    subroutine set_lit(self, lit)
        class(Polyhedron), intent(inout) :: self
        logical, intent(in) :: lit                          !! true: lit; false: unlit and transparent
        self%lit_ = lit
        call make_polyhedron(self)
    end subroutine set_lit

    !> Whether the faces are lit.
    pure logical function lit(self)
        class(Polyhedron), intent(in) :: self
        lit = self%lit_
    end function lit

    !> Shading smooth across the vertices (normals averaged over the faces
    !> around each: a curved surface, e.g. a grid's), or flat faces (the
    !> default: a polyhedron's).
    subroutine set_smooth(self, smooth)
        class(Polyhedron), intent(inout) :: self
        logical, intent(in) :: smooth                       !! true: smooth shading; false: flat faces
        self%smooth_ = smooth
        call make_polyhedron(self)
    end subroutine set_smooth

    !> Whether the shading is smooth.
    pure logical function smooth(self)
        class(Polyhedron), intent(in) :: self
        smooth = self%smooth_
    end function smooth

    !> Both sides of the faces drawn and lit (the default: open surfaces), or
    !> only their fronts (a closed surface, its faces counterclockwise seen
    !> from outside: translucent, its far side isn't drawn through it).
    subroutine set_two_sided(self, two_sided)
        class(Polyhedron), intent(inout) :: self
        logical, intent(in) :: two_sided                    !! true: both sides; false: fronts only
        self%two_sided_ = two_sided
        call make_polyhedron(self)
    end subroutine set_two_sided

    !> Whether both sides are drawn.
    pure logical function two_sided(self)
        class(Polyhedron), intent(in) :: self
        two_sided = self%two_sided_
    end function two_sided