作业九Luna(Vpython, matplotlib)

2013301020084 许晗
作业 双星 Euler-Cromer方法 Vpython

摘要　

月亮绕地球公转同时绕自身轴自转，并且一直保持同一面朝向地球．在月球初形成时，它的自转速度可能比公转速度大很多，但是随时间流逝，月球的质量分布可能会变得不均匀－根据受到地球引力影响的不同。可能到了某种分布状态就可以造成目前看不到月球背面。

正文　

程序目的是完成chapter4.20

计算原理

• 在质心系下，月亮一受力 $$\vec{F}_1=-\frac{GM_{earth}m_1}{r_1^3}\vec{r}_1 \tag1$$
  受到力矩 $$\vec{\tau}_1=(\vec{r}_1-\vec{r}_c) \times \vec{F}_1 \tag2$$
  月亮二受力也是类似．
  月亮一和二构成的系统在地球坐标系下受的力矩就是 $$\tau = \vec{\tau_1}+\vec{\tau_2}+\vec{\tau_c}=\vec{\tau_1}+\vec{\tau_2} \tag3$$
  因为质心受力指向地球与位矢方向一致，故 $\vec{\tau}_c=\vec{F}_c \times \vec{r}_c=0$。
• 然后，可以得到角动量的变化
  $$\frac{dJ}{dt}=\vec{\tau}_1+\vec{\tau}_2 \tag4$$
  $$\vec{J} = I\vec{\omega}_{total} = I(\vec{\omega}+\vec{\omega}_c) \tag5$$
  转动惯量 $I=m_1|r_1|^2+m_2|r_2|^2$，又因为$\frac{d\omega_c}{dt}=0$
  $$\frac{d\omega}{dt}=\frac{\vec{\tau}_1+\vec{\tau}_2}{I} \tag6$$
• 设月亮一、月亮二到质心的距离分别为$l1$、$l2$, 由质心的定义
  $$\vec{r}_c=\frac{m_1\vec{r_1}+m_2\vec{r_2}}{m_1+m_2} \tag7$$
  得到 $l_1/l_2=m_2/m_1$，利用这个条件以及$r_1\approx r_c,r_2\approx r_c$化简(6)式得
  $$\frac{d\omega}{dt}\approx -\frac{3GM_{earth}}{r_c^5}(x_csin\theta-y_ccos\theta)(x_ccos\theta+y_csin\theta) \tag 8$$

数值计算单位和方法

• 若取做圆周运动的半径长度(即$r_c$)为单位长度(HU)，再令公式中$GM_{earth}=4\pi^2$，则做一整圈圆周运动的时间为$T=2\pi\sqrt{\frac{r^3}{GM_{earth}}}=1$，故取做整圈圆周的时间为单位时间(yr)。
• Euler-Cromer方法计算质心的运动，以及$\theta$的变化，并通过式(7)给出月亮一和二的位置。

程序使用

目前是计算两次，计算的初值需要在程序里修改:

omega1 = [0] #here
theta1 = [0.00] #here
t1 = [0]
reset()
run(moon1, moon2, omega1, theta1, t1, dt=0.0001) #keyboard'q'
t2 = [0]
omega2 = [0] #here
theta2 = [0.01] #here
reset()
run(moon1, moon2, omega2, theta2, t2, dt=0.0001)

程序每运行1s(1个yr时间长度)，需要使用键盘：按'q'确认画图或进行下一个初值，其他键是继续这次计算。

结果展示与分析

计算取0.001(yr)为间隔，只计算了椭圆轨道。
1. - 按书上的要求展示$\Delta \theta=|\theta_1-\theta_2|$的变化：

- Vpython中的轨迹

2.
按题目4.20将程序中对角度的调整注释掉：

#if (theta <= - math.pi):    theta = theta + 2*math.pi
#if (theta > math.pi):    theta = theta - 2*math.pi `

• 得到连续的$\Delta \theta$变化
  
  可以看出在比较小的值时有与修改前一致的突降，因为这时两个角度都比较小，不用对角度调整。越到后来，两个角度越大，且两次计算的差值也更容易变大。需要说明在第二幅图中曲线超过10的部分，虽然变化幅度看起来很小，但是由于这是指数坐标，其实该角度的变化是非常大的．第二幅图像中$\Delta \theta$曲线的下跌时由与系统本身的性质决定。
• Vpython中的轨迹，当然与未修改前的一致
  
• 演示，计算了2(yr)时间中的运动