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.