Equation of Motion

* Türkçesi için;

In CR3BP, M1 and M2 are assumed very massive objects compared to the m3. So, gravitational effect of  to the third body can be neglected. The first two masses make circular Keplerian Motion together about their mutual centre of mass, so mean motion (ω) is constant and ω=2π/T (T is period of binary system; M1, M2). Therefore, the masses M1 and M2  affect the motion of m3, without in return being affect by themselves.

To simplify the problem, the equation of motion is arranged to non-dimensional form and equations motion can be written for one parameter, which is μ called mass ratio. Also becuase M1 and M2 make circular Keplerian motion in constant ω, then these bodies seem motionless in rotated axis, which rotates in “ω” speed relative to intertia axis. So non-dimensional equation motion can be practically written in rorated axis, see ;

Where μ1 and μ2 represent both non-dimensional masses of M1 and M2, and distance of M1 and M2 from O (their mutual centre of mass). And r1 and r2 are distance m3 from M1 and M3, respectivly.

Also, Ux, Uy, Uz represent derivation of Potantial Function (U), according to x,y,z; Such as Ux is derivation of U according to x, and so on. And the Potantial Function (U);

Also, Jacobi Integral, which is supplied from conversation of energy, provides an additional equation;

V²=2U-C

Where V is the velocity (V²=V²x+V²y+V²z), and C is Jacobi Constant, which represent energy constant. Simple Matlab code of Equation of Motion is below. Just run it and have the trajectory plot. You may also play the initial conditions and other parameter. Enjoy !

function CR3BP
global mue
% --- All in non-dimesional form ------

mue=0.01215; % mass ratio of Earth-Moon
xo=0.93; yo=0.00; zo=0.04; % positions 
T=50; % Time 
C=3.16; % Jacobi Constant

V=Velo(C,xo,yo,zo); % total velocity
Vx=0; Vy=-V; Vz=0; % velocities
Xo=[xo yo zo Vx Vy Vz]; % all initial conditions
options = odeset('RelTol',1e-11,'AbsTol',1e-11); % set timestep
[t,x]=ode45(@cr3bp_3d, [0 T], Xo ,options); 

 plot3(x(:,1),x(:,2),x(:,3))
axis equal
 grid on

%----- Velocity -------- 
function Vo=Velo(C,xi,yi,zi)
global mue ;
r10=sqrt((xi+mue).^2+yi.^2+zi.^2);r20=sqrt((xi-1+mue).^2+yi.^2+zi.^2);
OM0=0.5*(xi.^2+yi.^2)+(1-mue)/r10+mue/r20+0.5*mue*(1-mue);
Vo=sqrt(2*OM0-C);
%----- Equation Motion --------
function dx = cr3bp_3d(t, x)
global mue ;
dx=zeros(6,1);
r1=sqrt((x(1)+mue)^2+x(2)^2+x(3)^2);
r2=sqrt((x(1)-1+mue)^2+x(2)^2+x(3)^2);
dx(1)=x(4);
dx(2)=x(5);
dx(3)=x(6);
dx(4)=2*x(5)+x(1)-(1-mue)*(x(1)+mue)/r1^3-mue*(x(1)-1+mue)/r2^3;
dx(5)=-2*x(4)+x(2)*(1-(1-mue)/r1^3-mue/r2^3);
dx(6)=-x(3)*((1-mue)/r1^3+mue/r2^3);

You may also transfer the varibles from non-dimensional to dimensional form;

Where x,y,z and dot x,y,z are trajectory data in non-dimensional form, which can be obtained from the plot in above code. r12 is distance between M1 and M2 and ω=2π/T (T is period of binary system; M1, M2). Rxyz and Vxyz are are trajectory data in dimensional form.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s