如何在Python中将3D天文坐标转换为2D极坐标用于太阳系模拟

正在开发一个2D太阳系模拟程序，需要为行星设置初始数据，现有数据格式如下（仅关注距离(AU)、速度(m/s)、**方向(弧度)**列）：

计划通过Astropy获取1750年1月1日的真实天体数据，但coordinates.get_body_barycentric_posvel()返回的是3D位置和速度，直接忽略Z轴降维效果不佳，需要将其转换为符合现有格式的2D极坐标数据，且不打算改为3D模型。现有代码中transform_data函数为临时错误实现，寻求正确的转换方法或替代方案。

正确的做法是将3D坐标投影到太阳系黄道平面（行星公转的基准平面），而非直接丢弃Z轴，步骤如下：

• 投影到黄道平面：利用Astropy将天体坐标转换到黄道坐标系，提取XY分量，确保所有天体的位置和速度都在同一2D平面上
• 计算2D极坐标参数：基于投影后的XY分量计算距离、速度模长，以及极角方向

修改后的核心转换函数及逻辑如下：

def transform_data(pos, vel):
    # 将3D日心坐标转换为黄道坐标系
    sky_pos = coord.SkyCoord(pos, frame='icrs')
    sky_pos_ecl = sky_pos.transform_to('ecliptic')
    # 提取黄道平面的XY分量
    pos_ecl_xy = coord.CartesianRepresentation(sky_pos_ecl.cartesian.x, sky_pos_ecl.cartesian.y, 0*u.au)
    
    # 对速度做同样的黄道投影
    sky_vel = coord.SkyCoord(vel, frame='icrs')
    sky_vel_ecl = sky_vel.transform_to('ecliptic')
    vel_ecl_xy = coord.CartesianDifferential(sky_vel_ecl.cartesian.x, sky_vel_ecl.cartesian.y, 0*u.au/u.s)
    
    # 计算距离、速度和极角
    dist = pos_ecl_xy.norm().to(u.AU).value
    vel = vel_ecl_xy.norm().to(u.m/u.s).value
    angle = np.arctan2(sky_pos_ecl.cartesian.y.value, sky_pos_ecl.cartesian.x.value)
    
    return dist, vel, angle

另外，关于1750年数据的有效性：Astropy的JPL星历表支持回溯到公元前2000年左右，直接设置time = Time('1750-01-01')即可正常获取数据。

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


def transform_data(pos, vel):
    # 将3D日心坐标转换为黄道坐标系
    sky_pos = coord.SkyCoord(pos, frame='icrs')
    sky_pos_ecl = sky_pos.transform_to('ecliptic')
    # 提取黄道平面的XY分量
    pos_ecl_xy = coord.CartesianRepresentation(sky_pos_ecl.cartesian.x, sky_pos_ecl.cartesian.y, 0*u.au)
    
    # 对速度做黄道投影
    sky_vel = coord.SkyCoord(vel, frame='icrs')
    sky_vel_ecl = sky_vel.transform_to('ecliptic')
    vel_ecl_xy = coord.CartesianDifferential(sky_vel_ecl.cartesian.x, sky_vel_ecl.cartesian.y, 0*u.au/u.s)
    
    # 计算距离、速度和极角
    dist = pos_ecl_xy.norm().to(u.AU).value
    vel = vel_ecl_xy.norm().to(u.m/u.s).value
    angle = np.arctan2(sky_pos_ecl.cartesian.y.value, sky_pos_ecl.cartesian.x.value)
    
    return dist, vel, angle


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, angle = transform_data(pos, vel)

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


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

    time = Time('1750-01-01')

    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()

内容的提问来源于stack exchange，提问作者Tmpecho