In CR3BP, the third body (m3) is actually moving in a field which is also called effective potential, and the potential function (-U) gives this effective potential surface, which is shown in below;
As it is seen, L1, L2 and L3 are on the saddle in this surface, which means the motion are unstable in there. But L4 and L5 are on the peaks, which is actually dip in possitive U plot (+U), which means the motion are stable around there. The Matlab code of effective potantial is in below;
function pot_fun figure mue=0.2; %---mass ratio [X,Y] = meshgrid(-2:.01:2); hold on %potential function Z=-(X.^2+Y.^2)-(2.*(1-mue)./sqrt((X+mue).^2+Y.^2)+... 2.*mue./sqrt((X-1+mue).^2+Y.^2))-0*mue*(1-mue); surf(X,Y,Z,'FaceColor',[1 0.99 1],'EdgeColor','none') %[1 0.99 1] camlight left; lighting phong zlim([-3.8 -2.83]) alpha(0.5)% --- for transpancy view([60 80]) grid on
A horizontal cut of the plot of the effective potential (2U-C) gives Hill’s region, where the third body can only move in there at a given energy level C (Jacobi constant). Hill’s region is also called as zero velocity curves (or surfece in 3D), which is also obtined from Jacobi integral function (V²=2U-C), making velocity (V) zero.
In this figure, the three realms; namely the inner region, the neck region and the outher region that are connected by the neck region. For other values of the energy (C), one or more of these realms may be prohibited due to conservation of energy; such as, the necks may close off. The Matlab code of Hill’s region (or can be said as zero velocity curves in 2D or or zero velocity surfaces in 3D) is below. Copy-paste the code in your matlab work folder, then write in the commant line as follow;
- for 2 dimensional plots; zrvs(‘2d’, μ, C)
- for 3 dimensional plots; zrvs(‘3d’, μ, C, transparancy)
Where μ is mass ratio, C is jacobi constant. Transparancy is set from 0 to 1; 0 for clear, 1 for opaque (so set as 0.5). Such as zrvs(‘3d’,0.01215,3.16, 0.5).
function zrvs(dim,mue,C,tra)
% plot the zero velocity surface
%dim : '2d' or '3d'
%mue : mass ratio
%C : jacobi cost.
%tra : transparancy
%
%exp: zrvs('3d',0.01215,3.16,0.6)
hold on
switch dim
case '3d'
%---------------------------------------------------------------
d=-2:0.05:2;
[X,Y,Z] = meshgrid(d,d,d);
r1=((X + mue).^2+Y.^2+Z.^2).^(1/2);
r2=((X + mue-1).^2+Y.^2+Z.^2).^(1/2);
V=2*[(1/2)*(X.^2 + Y.^2)+(1-mue)./r1+mue./r2];
% you may change the color in below
h = patch(isosurface(X,Y,Z,V,C),'FaceColor','white','EdgeColor','none');
camlight; lighting phong
alpha(tra)
%view(3)
axis off
axis equal
%-----------------------------------------------------------------
case '2d'
%-----------------------------------------------------------------
syms X Y
f=(X.*X+Y.*Y)+2.*(1-mue)./sqrt((X+mue).^2+Y.^2)+...
2.*mue./sqrt((X-1+mue).^2+Y.^2)+mue.*(1-mue)-C;
ezcontour(f,[-2,2],[-1.5,1.5],300);
axis equal
grid on;
%-----------------------------------------------------------------
end
******** ho ho ho *****
Renkli Şeyler 