function [ftot,xmoon]=odefuncorbit(t,x);
%
% Objective: Computes the RHS of the ODE for the asteroid problem.
% That is, given the current asteroid state (position and velocity),
% compute the time-derivative of that state. Note that for simplicity,
% the position and velocity are each only two-dimensional.
%
% input variables:
%   t - time
%   x - 4x1 state. The state is comprised of a 2-dimensional position
%       appended with a 2-dimensional velocity, i.e., [x_1;x_2;v_1;v_2]
%       where [x_1;x_2] is position and [v_1;v_2] is velocity.
%
% output variables:
%   ftot - the derivative of the state wrt time at the current time.
%   xmoon - position of moon and the current time.
%
% functions called:
%   none
%
%
% Obtain the acceleration due to planet's gravity.
%
mu=200000; % planet Gm.
rad=norm(x(1:2,1),2);
c=-mu/(rad^3);
accp=c.*x(1:2,1);
%
% Obtain the current position of moon.
%
rmoon=45500; % moon dist from planet center
cmoon=mu/(rmoon^3);
ommoon=sqrt(cmoon); % moon angular rate
pmoon=2*pi/ommoon; % moon period (unused)
moonang=(-0.79*pi)+ommoon*t;    
xmoon=rmoon*[cos(moonang);sin(moonang)]; % current moon position
%
% Obtain the acceleration due to planet's moon.
%
mumoon=80000; % moon Gm
xdel=x(1:2,1)-xmoon; % position relative to the moon
radmoon2sc=norm(xdel,2);
cmoon2sc=-mumoon/(radmoon2sc^3);
accmoon=cmoon2sc.*xdel; % acceleration due to the moon
%
% Obtain total acceleration.
%
acctot=accp+accmoon;
%
% Generate the RHS of ODE.
%
ftot=[x(3:4,1);acctot];



