平面限制性三体问题仿真

平面限制性三体问题指第三个天体质量远小于其他两个天体的问题，即第三个天体的质量可以忽略。比如说地月系统和人造卫星就是典型的例子。本文旨在仿真卫星在拉格朗日点受到扰动后的行为。

需要说明的是地月形成的坐标系不是惯性系，卫星会受到引力、惯性力和科氏力。

clc;clear;clear;
%%%%%%%%%%%%%%——第二版——%%%%%%%%%%%%
%上一版是比较粗略的估计算法，以地球为地月系统中心点，实际上中心点偏离地心在地球内部
% 图设置
figure('name','地月系统拉格朗日点位置');
axis equal;
grid on;
hold on
xlabel('x/m');
ylabel('y/m');
hold on
axis([-5e8,5e8,-5e8,5e8]);
hold on
m1=5.976e24;%地球质量
m2=7.35e22;%月球质量
G=6.67*10^(-11);%万有引力常数
R=3.844e8;%地月距离
a1=m2*R/(m1+m2);
a2=m1*R/(m1+m2);
omega=sqrt(G*(m1+m2)/power(R,3));%牵连角速度
L1_x=R*0.85-a1;L1_y=0;%拉格朗日点坐标
L2_x=R*1.17-a1;L2_y=0;
L3_x=-R*0.99-a1;L3_y=0;
L4_x=R/2-a1;L4_y=sqrt(3)/2*R;
L5_x=R/2-a1;L5_y=-sqrt(3)/2*R;
x=[L1_x L2_x L3_x L4_x L5_x ];
y=[L1_y L2_y L3_y L4_y L5_y ];
plot(x,y,'ro');%绘图
hold on
line([-a1,L4_x],[0,L4_y],'linestyle',':');
line([-a1,L5_x],[0,L5_y],'linestyle',':');
line([L3_x,L2_x],[L3_y,L2_y],'linestyle',':');
hold on
theta=0:pi/30:2*pi;
x1=-a1+R*cos(theta);y1=R*sin(theta);
plot(x1,y1,'b--');
title('地月系统拉格朗日点','fontsize',18);
hold on
text(L1_x,L1_y+0.4e8,'L1','fontsize',12);
text(L2_x,L2_y+0.4e8,'L2','fontsize',12);
text(L3_x+0.5e8,L3_y,'L3','fontsize',12);
text(L4_x+0.5e8,L4_y,'L4','fontsize',12);
text(L5_x+0.5e8,L5_y,'L5','fontsize',12);
hold on
%%
dt=0.005;
pos=[L4_x,L4_y];
v=[40,50];
plant=plot(pos(1),pos(2));
h=animatedline('color','b');
frame=0;
while true
   rho1=sqrt((pos(1)+a1)^2+pos(2)^2); 
   rho2=sqrt((pos(1)-a2)^2+pos(2)^2); 
   r=sqrt(pos(1)^2+pos(2)^2);
   a_x=-G*m1*(pos(1)+a1)/power(rho1,3)-G*m2*(pos(1)-a2)/power(rho2,3)+omega^2*pos(1)+2*omega*v(2);
   a_y=-G*m1*pos(2)/power(rho1,3)-G*m2*pos(2)/power(rho2,3)+omega^2*pos(2)-2*omega*v(1);
   
   v(1)=v(1)+a_x*dt;
   v(2)=v(2)+a_y*dt;
   
   pos(1)=pos(1)+v(1)*dt+1/2*a_x*dt^2;
   pos(2)=pos(2)+v(2)*dt+1/2*a_y*dt^2;
   frame=frame+1;
   if frame==1000000
      frame=0;
      addpoints(h,pos(1),pos(2));
      drawnow
   end
end

下图是在L4点添加速度扰动的结果：