在Windows系统下用python编写一个程序：以太阳为参考，自行给定初始条件使天体绕太阳的轨道分别为圆形，椭圆，双曲线和抛物线，用四个子图画出各轨道随时间演化的动态图。不用FuncAnimation。

好的，以下是代码实现：

import numpy as np 
import matplotlib.pyplot as plt 

# Constants 
G = 6.674e-11 # gravitational constant 
m_sun = 1.989e30 # mass of the sun 

# Initial conditions 
r0 = [147e9, 147e9, 147e9, 147e9] # distance from the sun [m] 
v0 = [0, 30.29e3, 42.1e3, 42.1e3] # velocity with respect to the sun [m/s] 
theta0 = np.pi/2 # initial angle [rad] 
t0 = 0 # initial time [s] 
tf = 3.154e7 # final time [s] (1 year) 
dt = 3600 # time step [s] (1 hour) 

# Function to calculate acceleration 
def accel(r, m): 
    r_mag = np.linalg.norm(r) 
    r_hat = r/r_mag 
    a_mag = G*m/r_mag**2 
    a = -a_mag*r_hat 
    return a 

# Function to calculate orbit 
def orbit(r0, v0, t0, tf, dt, m): 
    t = [t0] 
    r = np.array(r0) 
    v = np.array(v0) 
    theta = [theta0] 
    
    for i in range(int((tf-t0)/dt)): 
        a = accel(r, m) 
        v += a*dt 
        r += v*dt 
        t.append(t[-1]+dt) 
        theta.append(np.arctan(r[1]/r[0])) 
    
    return np.array(theta), np.linalg.norm(r, axis=0), t 

# Circular orbit 
theta_circ, r_circ, t_circ = orbit(r0, v0[0], t0, tf, dt, m_sun) 

# Elliptical orbit 
theta_ellip, r_ellip, t_ellip = orbit(r0, v0[1], t0, tf, dt, m_sun) 

# Hyperbolic orbit 
theta_hyp, r_hyp, t_hyp = orbit(r0, v0[2], t0, tf, dt, m_sun) 

# Parabolic orbit 
theta_para, r_para, t_para = orbit(r0, v0[3], t0, tf, dt, m_sun) 

# Plotting 
fig, axs = plt.subplots(nrows=2, ncols=2) 

axs[0,0].set_title('Circular Orbit') 
axs[0,0].set_xlim([-2, 2]) 
axs[0,0].set_ylim([-2, 2]) 
axs[0,0].set_aspect('equal') 
axs[0,0].set_xlabel('x [km]') 
axs[0,0].set_ylabel('y [km]') 

axs[0,1].set_title('Elliptical Orbit') 
axs[0,1].set_xlim([-2e8, 2e8]) 
axs[0,1].set_ylim([-2e8, 2e8]) 
axs[0,1].set_aspect('equal') 
axs[0,1].set_xlabel('x [km]') 
axs[0,1].set_ylabel('y [km]') 

axs[1,0].set_title('Hyperbolic Orbit') 
axs[1,0].set_xlim([-7e8, 7e8]) 
axs[1,0].set_ylim([-7e8, 7e8]) 
axs[1,0].set_aspect('equal') 
axs[1,0].set_xlabel('x [km]') 
axs[1,0].set_ylabel('y [km]') 

axs[1,1].set_title('Parabolic Orbit') 
axs[1,1].set_xlim([-2e9, 2e9]) 
axs[1,1].set_ylim([-2e9, 2e9]) 
axs[1,1].set_aspect('equal') 
axs[1,1].set_xlabel('x [km]') 
axs[1,1].set_ylabel('y [km]') 

for ax in axs.flat: 
    ax.grid() 

for i in range(len(t_circ)): 
    circ_x = r_circ*np.cos(theta_circ) 
    circ_y = r_circ*np.sin(theta_circ) 
    axs[0,0].cla() 
    axs[0,0].plot(circ_x/1e9, circ_y/1e9, 'b') 
    axs[0,0].plot([0], [0], 'yo') 

for i in range(len(t_ellip)): 
    ellip_x = r_ellip*np.cos(theta_ellip) 
    ellip_y = r_ellip*np.sin(theta_ellip) 
    axs[0,1].cla() 
    axs[0,1].plot(ellip_x/1e9, ellip_y/1e9, 'b') 
    axs[0,1].plot([0], [0], 'yo') 

for i in range(len(t_hyp)): 
    hyp_x = r_hyp*np.cos(theta_hyp) 
    hyp_y = r_hyp*np.sin(theta_hyp) 
    axs[1,0].cla() 
    axs[1,0].plot(hyp_x/1e9, hyp_y/1e9, 'b') 
    axs[1,0].plot([0], [0], 'yo') 

for i in range(len(t_para)): 
    para_x = r_para*np.cos(theta_para) 
    para_y = r_para*np.sin(theta_para) 
    axs[1,1].cla() 
    axs[1,1].plot(para_x/1e9, para_y/1e9, 'b') 
    axs[1,1].plot([0], [0], 'yo') 

plt.show() 

解释一下主要的代码：

1. 定义了常量“G”和“m_sun”，分别表示万有引力常数和太阳的质量。
2. 定义了四个不同轨道的初始条件：距离太阳的距离、相对于太阳的速度、初始角度、初始时间、最终时间和时间步长。
3. 定义了“accel”函数来计算加速度，其中“r_mag”表示距离的大小，“r_hat”表示距离的单位向量，“a_mag”表示加速度的大小，“a”表示加速度的向量。
4. 定义了“orbit”函数来计算轨道，其中“t”表示时间，“r”表示距离向量，“v”表示速度向量，“theta”表示角度。该函数使用欧拉方法来计算轨道，其中加速度由“accel”函数计算得出。
5. 调用“orbit”函数来计算四个不同轨道的轨迹。
6. 创建四个子图，每个子图分别表示不同的轨道。设置子图的标题、x轴和y轴标签、x轴和y轴的范围、网格和纵横比等。
7. 使用循环绘制四个轨道的动态图像。在每个时间步长中，计算每个轨道的位置，并在相应的子图中绘制它。同时，在每个子图中绘制太阳（用黄点表示）。

注意：这个程序只是一个简单的示例，使用欧拉方法计算轨道不是最准确的方法，而且每个轨道的时间步长都是固定的。如果您需要更准确和更高效的方法，请考虑使用其他轨道计算方法，如Runge-Kutta方法等。