Python 中两颗行星围绕一颗恒星的轨道，速度为 verlet

Question

我正在做一个模拟，其中有两个行星移动和一个固定的恒星。假设两颗行星中的一颗行星 p1 围绕恒星运行，而另一颗行星 p2 围绕 p1 运行（因此也围绕恒星运行）。该脚本是用 python 编写的，我使用的算法是：velocity-verlet

到目前为止我写了以下脚本

    import numpy as np
    from matplotlib import pyplot as plt
    import time as tm

    start_time = tm.time()

    def v_verlet(x_p1, y_p1, vx_p1, vy_p1, x_p2, y_p2, vx_p2, vy_p2, r_TS, r_TP, r_PS, Nstep, tau):

        GM = G = 4*np.pi**2 # G*M in astronomical unit. We suppose M to be the mass of the sun and M = 1
        #if we want GM with M the mass of a generic star, we do M*4*np.pi**2
    
        #arrays to store position of the first planet p1
        xP1_ARR = np.zeros(Nstep+1)
        yP1_ARR = np.zeros(Nstep+1)
        vxP1_ARR = np.zeros(Nstep+1)
        vyP1_ARR = np.zeros(Nstep+1)
    
        #arrays to store position of the second planet p2
        xP2_ARR = np.zeros(Nstep+1)  
        yP2_ARR = np.zeros(Nstep+1)
        vxP2_ARR = np.zeros(Nstep+1)
        vyP2_ARR = np.zeros(Nstep+1)
    
        #saving initial conditions for the two planets
        xP1_ARR[0] = x_p1
        yP1_ARR[0] = y_p1
        vxP1_ARR[0] = vx_p1
        vyP1_ARR[0] = vy_p1
    
        xP2_ARR[0] = x_p2
        yP2_ARR[0] = y_p2
        vxP2_ARR[0] = vx_p2
        vyP2_ARR[0] = vy_p2
    
        #velocity-verlet
        for i in range(Nstep):
    
            #acceleration p1
            Fx_p1 = (-GM*x_p1/r_TS**3) + (-G*P*(x_p1-x_p2)/r_TP**3)
            Fy_p1 = (-GM*y_p1/r_TS**3) + (-G*P*(y_p1-y_p2)/r_TP**3)
    
            #acceleration p2
            Fx_p2 = (-GM*x_p2/r_PS**3) + (-G*T*(x_p2-x_p1)/r_TP**3)
            Fy_p2 = (-GM*y_p2/r_PS**3) + (-G*T*(y_p2-y_p1)/r_TP**3)
    
            #update position p1
            x_p1 += vx_p1*tau+0.5*Fx_p1*tau**2
            y_p1 += vy_p1*tau+0.5*Fy_p1*tau**2
    
            #update position p2
            x_p2 += vx_p2*tau+0.5*Fx_p2*tau**2
            y_p2 += vy_p2*tau+0.5*Fy_p2*tau**2
    
            #update distances
            r_P1S = np.sqrt(x_p1**2+y_p1**2)
            r_P1S = np.sqrt(x_p2**2+y_p2**2)
            r_P1P2 = np.sqrt((x_p1-x_p2)**2 + (y_p1-y_p2)**2)
    
            #update acceleration p1
            F1x_p1 = (-GM*x_p1/r_TS**3) + (-G*P*(x_p1-x_p2)/r_TP**3)
            F1y_p1 = (-GM*y_p1/r_TS**3) + (-G*P*(y_p1-y_p2)/r_TP**3)
    
            #update acceleration p2 i+1
            F1x_p2 = (-GM*x_p2/r_PS**3)+(-G*T*(x_p2-x_p1)/r_TP**3)
            F1y_p2 = (-GM*y_p2/r_PS**3)+(-G*T*(y_p2-y_p1)/r_TP**3)
    
            #update velocity p1
            vx_p1 += 0.5*tau*(Fx_p1+F1x_p1)
            vy_p1 += 0.5*tau*(Fy_p1+F1y_p1)
    
            #update velocity p2
            vx_p2 += 0.5*tau*(Fx_p2+F1x_p2)
            vy_p2 += 0.5*tau*(Fy_p2+F1y_p2)
    
            #saving all
            xP1_ARR[i+1] = x_p1
            yP1_ARR[i+1] = y_p1
            vxP1_ARR[i+1] = vx_p1
            vyP1_ARR[i+1] = vy_p1
    
            xP2_ARR[i+1] = x_p2
            yP2_ARR[i+1] = y_p2
            vxP2_ARR[i+1] = vx_p2
            vyP2_ARR[i+1] = vy_p2
    
        return xP1_ARR, yP1_ARR, xP2_ARR, yP2_ARR

    #Constants
    M = 2e30 # sun mass
    T = 5.9722E+24/M # earth mass
    P = 7.348e22/M #moon mass

    #initial position for the star
    x_s = 0
    y_s = 0

    #initial position and velocity in UA for the first planet p1
    x_p1 = 1
    y_p1 = 0.0
    r_P1S = np.sqrt(x_p1**2+y_p1**2) # distance first planet to the star
    vx_p1 = 0
    vy_p1 = 2\*np.pi #casual

    #initian position and velocity in UA for the second planet p2
    x_p2 = 1.1
    y_p2 = 0
    vx_p2 = 0
    vy_p2 = np.pi

    r_P1P2 = np.sqrt((x_p1-x_p2)**2 + (y_p1-y_p2)2) #distance from p1 to p2
    r_P2S = np.sqrt(x_p2**2 + y_p22) #distance second planet to the star

    #number of steps and timestep
    Nstep = 100000
    tau = 0.0001 # timestep

    x_p1_ARR, y_p1_ARR, X_p2_ARR, Y_p2_ARR = v_verlet(x_p1, y_p1, vx_p1, vy_p1, x_p2, y_p2, vx_p2, vy_p2, r_P1S, r_P1P2, r_P2S, Nstep, tau)

    #Plot trajectories
    plt.xlim(-3,3)
    plt.ylim(-3,3)
    plt.xlabel("x (AU)")
    plt.ylabel("y (AU)")
    #initial position
    plt.scatter(x_s, y_s, color = "orange")
    plt.scatter(x_p1_ARR\[0\], y_p1_ARR\[0\], color = "brown")
    plt.scatter(X_p2_ARR\[0\], Y_p2_ARR\[0\], color = "olive")
    #trajectories
    plt.plot(x_p1_ARR, y_p1_ARR, color = "maroon")
    plt.plot(X_p2_ARR, Y_p2_ARR, color = "darkgrey")

    plt.axis("equal")

    plt.grid(alpha = 0.3)

    end_time = tm.time()

    print("Time: ", end_time-start_time)

    plt.show()

但我得到这个结果：

但是，我想要像this

这样的东西

有人能帮我理解为什么p1星球不拉p2吗？ 如何让p2绕p1转？