Add new `angular` coordinate system, and use as a basis for `xy` and `km` coordinate systems

[ortk95]

The radec coordinate system has singularities at the celestial poles (declination = ±90), causing warping when e.g. plotted on a cartesian axis (see #323). This is inherent in the radec coordinate system, so probably something that isn't worth (or necessarily even possible) to change. However, the km and xy coordinates are currently both just affine transformations of the radec coordinates, which means that they both inherit the warping at high declination angles when they shouldn't. Therefore, we should change how these 'tangent frame' coordinates are calculated to fix this issue.

This is possibly just fixable by rotating the spherical coordinates of the target away from the pole before doing the tangent plane projection to the km/xycoordinates (though will need to properly check and test this).

Things that need doing:

• [x] Implement new tangent coordinate system based on tangent plane centred on target body
• [x] Settle on a good name for this new tangent coordinate system
• [x] Decide on what units to use for tangent coordinate system - maybe arcsec (then name system arcsec for consistency with km)?
• [x] Allow the origin of the tangent coordinates to be customised - default to target centre, but would be useful if all relevant functions could take appropriate origin_ra, origin_dec and rotation kwargs for full customisability (also need to make sure these kwargs don't conflict with anything else important)
• [x] Use tangent coords as basis for km coords
• [x] Use tangent coords as basis for xy coords
• [x] Update relevant code in km and xy wireframe plots, matplotlib transforms etc. - will probably require adding some function to maybe transform ra, dec coordinates before passing them to the matplotlib transforms
• [x] Add wireframe plot for new coordinates
• [x] Update any relevant documentation for general coordinate systems and specifics to km and xy coordinates
• [x] Check image backplanes (which are effectively all based on the xy coordinate system) work properly
• [x] Update matplotlib transforms
• [x] Update xy wireframe
• [x] Make transforms more efficient (e.g. remove any round-trips to radec)
• [x] Check names don't conflict with any existing naming system (both the name of the coordinate system, and any kwargs used to construct it)

[ortk95]

Looks like this should be doable by defining a new coordinate system using a tangent plane centred on the target body. This new tangent coordinate system can then be used as a basis for the xy and km coordinate systems with simple affine transformations (i.e. a drop-in replacement of radec).

import functools
from typing import cast

import matplotlib.pyplot as plt
import numpy as np
import spiceypy as spice

import planetmapper

mkpath = './scratchpad/JUICE/kernels/mk/juice_plan.tm'
planetmapper.base.prevent_kernel_loading()
spice.kclear()
spice.furnsh(mkpath)

times = [
    '2032-08-22T00:00:00',
    '2032-09-01T00:00:00',
    '2032-09-24T05:30:00',
    '2032-09-24T06:30:00',
]

bodies = [
    planetmapper.Body(
        'Jupiter',
        time,
        observer='JUICE_SPACECRAFT',
        observer_frame='JUICE_SPACECRAFT_PLAN',
    )
    for time in times
]
bodies.extend(
    [
        planetmapper.Body(
            'Jupiter',
            times[0],
            observer='JUICE_SPACECRAFT',
            observer_frame='J2000',
        ),
        planetmapper.Body(
            'Jupiter',
            times[0],
        ),
    ]
)

fig, axs = plt.subplots(nrows=len(bodies), ncols=2, figsize=(8, 12), squeeze=False)
axs = cast(list[list[plt.Axes]], axs)

def targvec2tangent(
    targvec: np.ndarray, body: planetmapper.Body
) -> tuple[float, float]:
    if any(np.isnan(targvec)):
        return np.nan, np.nan

    matrix = get_transform_matrix(body)
    vec = matrix @ body._targvec2obsvec(targvec)
    _, x, y = spice.recrad(vec)
    x = np.rad2deg(x)
    x *= -1
    x %= 360
    if x > 180:
        x -= 360
    y = np.rad2deg(y)

    return x, y

@functools.cache
def get_transform_matrix(body: planetmapper.Body) -> np.ndarray:
    origin_vec = body._targvec2obsvec(np.array([0.0, 0.0, 0.0]))

    _, ra_angle, _ = spice.recrad(origin_vec)
    ra_matrix = spice.rotate(ra_angle, 3)

    _, _, dec_angle = spice.recrad(ra_matrix @ origin_vec)
    dec_matrix = spice.rotate(-dec_angle, 2)

    return dec_matrix @ ra_matrix

for body_idx, body in enumerate(bodies):
    ax = axs[body_idx][0]
    body.plot_wireframe_radec(
        ax=ax,
        add_title=False,
        add_axis_labels=False,
    )
    target_targvec = body.lonlat2targvec(
        *body.radec2lonlat(body.target_ra, body.target_dec)
    )
    ax.scatter(
        *body.lonlat2radec(*body.targvec2lonlat(target_targvec)), color='b', marker='x'
    )

    ax = axs[body_idx][1]
    target_targvec = np.array([0.0, 0.0, 0.0])

    for targvec in body._limb_targvec(npts=30):
        x, y = targvec2tangent(targvec, body)
        ax.scatter(x, y, color='r', marker='.')

    for ra, dec in zip(*body.terminator_radec(npts=30)):
        x, y = targvec2tangent(body.lonlat2targvec(*body.radec2lonlat(ra, dec)), body)
        ax.scatter(x, y, color='g', marker='.')

    for lon, lat, s in body.get_poles_to_plot():
        x, y = targvec2tangent(body.lonlat2targvec(lon, lat), body)
        ax.annotate(s, (x, y), ha='center', va='center')

    target_tangent = targvec2tangent(target_targvec, body)
    ax.scatter(*target_tangent, color='b', marker='x')
    fmt = 'g'
    print(
        f'({body.target_ra:.3f}, {body.target_dec:.3f}) -> ({target_tangent[0]:+{fmt}}, {target_tangent[1]:+{fmt}})'
    )
    ax.set_aspect('equal')

    radius = body.target_diameter_arcsec / 2 / 60 / 60
    for r in (-radius, 0, radius):
        ax.axvline(r, color='c', zorder=0, alpha=0.25)
        ax.axhline(r, color='c', zorder=0, alpha=0.25)

plt.show()