How to flatten 3d astronomical coordinates into 2d polar form in Python

Question

I first asked this at Stack Overflow but got redirected here

I am making a 2d simulation of the solar system. For that you need a starting point for all planets. Right now I am storing that data in this form (ignore the color, radius, mass and Max_Trail_Length columns as they are cosmetics for the model):

Name	Color (RGB)	Radius (km)	Mass (kg)	Distance (AU)	Velocity (m/s)	Direction (radians)	Max_Trail_Length
Sun	(255, 255, 0)	6963400	1.98847e+30	0	0	0	100
Mercury	(255, 0, 0)	2440	3.3011e+23	0.393	47362	pi / 2	200
Venus	(255, 165, 0)	6052	4.8675e+24	0.723	35021.4	pi / 2	600
Earth	(0, 0, 255)	6371	5.9722e+24	1	29784.8	pi / 2	800
Mars	(255, 0, 0)	3390	6.4171e+23	1.52	24130.8	pi / 2	2000
Jupiter	(255, 222, 173)	69911	1.8982e+27	5.203	13070	pi / 2	10000
Saturn	(210, 180, 140)	58232	5.6834e+26	9.737	9690	pi / 2	20000
Uranus	(0, 0, 128)	25362	8.681e+25	19.61	6810	pi / 2	50000
Neptune	(0, 0, 255)	24622	1.02413e+26	29.897	5430	pi / 2	200000

As you can see, the direction is the same, so all planets are aligned, and the distance and velocity data is just the average.

For the record, I do not plan to make this into a 3d model, even though it would mitigate this problem entirely.

I want to collect real data from a date (January 1. 1750) using astropy, however when i collect the position and velocity using pos, vel = coordinates.get_body_barycentric_posvel(name, time) (coordinates from astropy), I get 3d data, as you can expect. But I want to convert it into data similar as I at the moment have.

I can't seem to do it right, as flattening the z axis doesn't work (I though it might be ignored at first as the solar system is relatively flat).

Here is the code I have right now for collecting this data.

import csv
import os

import numpy as np
from astropy import coordinates as coord
from astropy import units as u
from astropy.time import Time


# FIXME: Temporary wrong solution
def transform_data(pos, vel):
    pos = coord.CartesianRepresentation(pos.x, pos.y, 0 * u.au)
    vel = coord.CartesianDifferential(vel.x, vel.y, 0 * u.au / u.s)
    vel = vel.norm().to(u.m / u.s).value
    dist = pos.norm().to(u.AU).value

    return dist, vel


def get_path():
    parent_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
    data_dir = os.path.join(parent_dir, 'data')
    os.makedirs(data_dir, exist_ok=True)
    file_path = os.path.join(data_dir, 'solsystem_data_1750.csv')
    return file_path


def process_data(time, writer):
    for name in ['Sun', 'Mercury', 'Venus', 'Earth', 'Mars', 'Jupiter', 'Saturn', 'Uranus', 'Neptune']:
        pos, vel = coord.get_body_barycentric_posvel(name, time)

        dist, vel = transform_data(pos, vel)

        angle = np.arctan2(pos.y, pos.x).value

        writer.writerow([name, f'{dist:.5f}', f'{vel:.5f}', f'{angle:.5f}'])


def main():
    coord.solar_system_ephemeris.set('jpl')

    time = Time('1970-01-01')  # TODO: Change to 1750-01-01, but 1750 is a dubious year

    file_path = get_path()

    with open(file_path, 'w', newline='') as file:
        writer = csv.writer(file)
        writer.writerow(['Name', 'Distance (AU)', 'Velocity (m/s)', 'Direction (radians)'])

        process_data(time, writer)


if __name__ == '__main__':
    main()

I want to know if I can either modify my transform_data function, or collect the data in some sort of other way.

I couldn't find any 2d data from astropy, and I haven't found a way to convert the coordinates in a good way yet. I researched some alternatives to astropy, like just collecting the data from NASA's Horizons System, however I there run into the same problem of flattening the data.

Answer 1

The code in question outputs:

Name,Distance (AU),Velocity (m/s),Direction (radians)
Earth,0.90583,30207.28054,1.76660

Direction is reasonable for January 1 but Distance is too short. get_body_barycentric_posvel returns ICRS coordinates, which are equatorial. If I compute pos.z / pos.y for each planet, all are between 0.30 and 0.45, and Earth's is 0.4336 = tan 23.44°.

It sounds like you want ecliptic coordinates instead. BarycentricMeanEcliptic may help. Its default representation is spherical. If you just want the planets' speeds, you can leave their velocities in equatorial form, but include the z component. I suggest something like:

pos_eq, vel_eq = coord.get_body_barycentric_posvel(name, time)

pos_ecl = coord.BarycentricMeanEcliptic(pos_eq)
dist = pos_ecl.distance.to(u.AU).value
angle = pos_ecl.lon.to(u.radian).value

vel = vel_eq.norm().to(u.m / u.s).value

If you assume the planets' initial radial velocities are 0, that's true only at perihelion and aphelion, so their orbits won't be exactly right, but at least they won't look too strange.