在Windows系统下用python编写一个的程序：在极坐标系下当天体绕太阳运动时，分别画出天体总能量大于，等于，小于零时的运动轨迹。

以下是一个实现的示例程序：

import numpy as np
import matplotlib.pyplot as plt

# 定义常量
G = 6.67408e-11    # 万有引力常数
M_sun = 1.9891e30   # 太阳质量
AU = 149.6e9        # 天文单位

# 定义运动函数
def motion(r, v):
    # 计算距离和加速度
    d = np.sqrt(np.sum(r**2))
    a = -G * M_sun / d**2 * r / d
    # 计算速度变化
    v = v + a * dt
    # 计算位置变化
    r = r + v * dt
    return r, v

# 定义初始条件
r0 = np.array([AU, 0])   # 初始位置
v0 = np.array([0, 30e3])   # 初始速度
E0 = 1/2 * np.sum(v0**2) - G * M_sun / np.sqrt(np.sum(r0**2))   # 初始总能量

# 定义时间步长和时间间隔
dt = 3600    # 时间步长
T = 365 * 24 * 3600   # 时间间隔

# 初始化轨迹和能量数组
r_traj = [r0]
v_traj = [v0]
E_traj = [E0]

# 计算轨迹和能量
t = 0
while t < T:
    r, v = motion(r_traj[-1], v_traj[-1])
    r_traj.append(r)
    v_traj.append(v)
    E = 1/2 * np.sum(v**2) - G * M_sun / np.sqrt(np.sum(r**2))
    E_traj.append(E)
    t += dt

# 转换为极坐标系
theta = np.arctan2([r[1] for r in r_traj], [r[0] for r in r_traj])
r = [np.sqrt(np.sum(r**2)) for r in r_traj]

# 绘制图形
fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
ax.plot(theta, r)
ax.set_rmax(3*AU)
ax.grid(True)

# 绘制能量曲线
fig, ax = plt.subplots()
ax.plot(np.arange(0, T, dt), E_traj)
ax.axhline(y=0, color='r')
ax.set_xlabel('Time (s)')
ax.set_ylabel('Total Energy (J)')
ax.set_title('Motion of a planet around the Sun')
plt.show()

这个程序可以画出天体在极坐标系下的运动轨迹，并可以同时绘制天体总能量随时间的变化曲线。你可以根据需要修改初始条件、时间间隔和时间步长等参数，来得到不同的轨迹和能量变化曲线。