XML_makerを完成させて少し時間ができたので，太陽系の動きを再現してみました。面倒だったので閏年は考慮していませんが，見た目上は問題なく動きます。

追記：再現とかいってますけど、公転周期しか再現してません。

こんな感じ

参考にしたサイト

軌道長半径などの定数はWikipediaの数値を使用しました。

import numpy as np
import matplotlib.pyplot as plt

def solar_system(year):
    day = 365 * year
    year = 1
    day_num = 0

    # 軌道長半径（au表記）1[au] = 149597870700[m]
    Mercury_orbit  = 0.38709927
    Venus_orbit    = 0.72333566
    Earth_orbit    = 1.00000261
    Mars_orbit     = 1.52371034
    Jupiter_orbit  = 5.20288700
    Saturn_orbit   = 9.53667594
    Uranus_orbit   = 19.18916464
    Neputune_orbit = 30.06992276

    for day in range(day):
        day += 1
        day_num += 1
        print('\r','Year'+':'+str(year) ,'Day'+':'+str(day_num),end='    ')
        if day % 365 == 0:
            year += 1
            day_num = 0

        n = 1000
        plot_range = 100  # -100~100で表示

        x = np.linspace(-plot_range, plot_range, n)
        y = np.linspace(-plot_range, plot_range, n)

        X, Y = np.meshgrid(x, y)
        Z = np.sqrt(X**2 + Y**2)

        plt.title('Solar System')
        plt.plot(0, 0, marker='o')  # 太陽

        # 出力半径
        Mercury_radius  = (Mercury_orbit  / Neputune_orbit) * plot_range
        Venus_radius    = (Venus_orbit    / Neputune_orbit) * plot_range
        Earth_radius    = (Earth_orbit    / Neputune_orbit) * plot_range
        Mars_radius     = (Mars_orbit     / Neputune_orbit) * plot_range
        Jupiter_radius  = (Jupiter_orbit  / Neputune_orbit) * plot_range
        Saturn_radius   = (Saturn_orbit   / Neputune_orbit) * plot_range
        Uranus_radius   = (Uranus_orbit   / Neputune_orbit) * plot_range
        Neputune_radius = (Neputune_orbit / Neputune_orbit) * plot_range

        # 軌道線をプロット
        plt.contour(X, Y, Z, levels=[Mercury_radius,
                                     Venus_radius,
                                     Earth_radius,
                                     Mars_radius,
                                     Jupiter_radius,
                                     Saturn_radius,
                                     Uranus_radius,
                                     Neputune_radius])

        pi = np.pi
        # 公転周期
        Mercury_radian  = (day/87.965)   * 2*pi
        Venus_radian    = (day/224.275)  * 2*pi
        Earth_radian    = (day/365)      * 2*pi
        Mars_radian     = (day/686.565)  * 2*pi
        Jupiter_radian  = (day/4328.9)   * 2*pi
        Saturn_radian   = (day/10752.9)  * 2*pi
        Uranus_radian   = (day/30663.65) * 2*pi
        Neputune_radian = (day/60148.35) * 2*pi

        # --------各天体の座標---------- #
        Mercury_x = round(Mercury_radius, 5) * np.cos(round(Mercury_radian, 5))
        Mercury_y = round(Mercury_radius, 5) * np.sin(round(Mercury_radian, 5))
        plt.plot(Mercury_x, Mercury_y, marker='o')

        Venus_x = round(Venus_radius, 5) * np.cos(round(Venus_radian, 5))
        Venus_y = round(Venus_radius, 5) * np.sin(round(Venus_radian, 5))
        plt.plot(Venus_x, Venus_y, marker='o')

        Earth_x = round(Earth_radius, 5) * np.cos(round(Earth_radian, 5))
        Earth_y = round(Earth_radius, 5) * np.sin(round(Earth_radian, 5))
        plt.plot(Earth_x, Earth_y, marker='o')

        Mars_x = round(Mars_radius, 5) * np.cos(round(Mars_radian, 5))
        Mars_y = round(Mars_radius, 5) * np.sin(round(Mars_radian, 5))
        plt.plot(Mars_x, Mars_y, marker='o')

        Jupiter_x = round(Jupiter_radius, 5) * np.cos(round(Jupiter_radian, 5))
        Jupiter_y = round(Jupiter_radius, 5) * np.sin(round(Jupiter_radian, 5))
        plt.plot(Jupiter_x, Jupiter_y, marker='o')

        Saturn_x = round(Saturn_radius, 5) * np.cos(round(Saturn_radian, 5))
        Saturn_y = round(Saturn_radius, 5) * np.sin(round(Saturn_radian, 5))
        plt.plot(Saturn_x, Saturn_y, marker='o')

        Uranus_x = round(Uranus_radius, 5) * np.cos(round(Uranus_radian, 5))
        Uranus_y = round(Uranus_radius, 5) * np.sin(round(Uranus_radian, 5))
        plt.plot(Uranus_x, Uranus_y, marker='o')

        Neputune_x = round(Neputune_radius, 5) * np.cos(round(Neputune_radian, 5))
        Neputune_y = round(Neputune_radius, 5) * np.sin(round(Neputune_radian, 5))
        plt.plot(Neputune_x, Neputune_y, marker='o')

        plt.gca().set_aspect('equal')
        plt.pause(0.01)
        plt.draw()
        plt.cla()

# exsample of use
year = int(input('何年分を見る？'))
solar_system(year)

グラフの大きさは19行目のplot_rangeを弄れば変更できます。