Trouble Rotating Coordinate System

Astropy
1

First post here: hello, hi, howdy!

I have a script which is designed to give the coordinates for solar system objects, but it outputs based on a coordinate system which has the Z-axis running parallel to Earth’s tilt and I’d like to use the ecliptic coordinate system. Now I programmed in a function to rotate the coordinates, but I’m assuming there’s a functionality similar to this built into Astropy and I may have missed it. Here is my code:

# Import
from astropy.coordinates import solar_system_ephemeris, get_body_barycentric, ICRS, PrecessedGeocentric
from astropy.time import Time
import astropy.units as u
import numpy as np


# Set Variables
start_date = '2023-12-25 00:00:00'
target_body = 'earth'
target_step_unit = 'day'
target_duration = 72
target_step_size = 6
cap = True


# Get Positions Function
def get_positions(body, start, duration, step_size, step_unit, if_cap):

    # Reformat Start Time
    start_time = Time(start, scale='tdb')

    # Time Mapping
    unit_mapping = {
        'second': u.second,
        'minute': u.minute,
        'hour': u.hour,
        'day': u.day,
        'week': u.week,
        'year': u.year,
        'year_sidereal': 365.25636 * u.day
    }

    # Set Steps and Duration
    step_duration = step_size * unit_mapping[step_unit]
    total_duration = duration * unit_mapping[step_unit]
    num_steps = int(total_duration / step_duration) + (1 if if_cap else 0)

    # Initialize Array
    positions = []

    # Set the Ephemeris to Heliocentric
    with solar_system_ephemeris.set('builtin'):
        for step in range(num_steps):

            # Calculate Current Step Time
            current_time = start_time + step * step_duration

            # Get Body Position - Heliocentric Frame
            body_position = get_body_barycentric(body, current_time)

            # Extract the XYZ Coordinates
            xyz_km = np.array([body_position.x.to(u.km).value,
                               body_position.y.to(u.km).value,
                               body_position.z.to(u.km).value])

            # Rotate the Coordinate System (Earth Tilt)
            xyz_rotated = rotate_coordinates(xyz_km, -23.4393)
            x_rotated, y_rotated, z_rotated = [u.Quantity(coord, unit=u.km) for coord in xyz_rotated]

            # Append the rotated position as a tuple
            positions.append((current_time, (x_rotated, y_rotated, z_rotated)))

    return positions


# Rotate System Function
def rotate_coordinates(coords, angle_deg):

    # Convert Degrees to Radians
    angle_rad = np.radians(angle_deg)

    # Rotation Matrix Around the X-axis
    rotation_matrix = np.array([
        [1, 0, 0],
        [0, np.cos(angle_rad), -np.sin(angle_rad)],
        [0, np.sin(angle_rad), np.cos(angle_rad)]
    ])

    # Apply the Rotation
    return np.dot(rotation_matrix, coords)


# Run: Retrieve Positions
target_positions = get_positions(
    target_body,
    start_date,
    target_duration,
    target_step_size,
    target_step_unit,
    cap)


# Rounding and Display
for time, pos in target_positions:
    x = f"{pos[0].value:.2f}"
    y = f"{pos[1].value:.2f}"
    z = f"{pos[2].value:.2f}"

    print(f"Time: {time.iso}, X: {x} km, Y: {y} km, Z: {z} km")

Thank you!