function xprime = traj_prime(t, x)
% function traj_prime.m -- specifies the first-order differential equations used in 
% integration. Takes into account gravity and rotation, and is based on equations of motion.

% Serina, last modified 21 June 2004
% potential_mat.m

rp = [x(1) x(3) x(5)];
global density_1 density_sur ell_x ell_y ell_z lay_depth omega
global density ellip_x ellip_y ellip_z laplacian

% for inner layer: solid, fractured body
density = density_1;
ellip_x = ell_x - lay_depth;
ellip_y = ell_y - lay_depth;
ellip_z = ell_z - lay_depth;
[gforce,jacobian,laplacian,U] = potential_mat(rp);
Ux = gforce(1);
Uy = gforce(2);
Uz = gforce(3);

% for surface layer: regolith of depth layer_depth
density = density_sur;
ellip_x = ell_x;
ellip_y = ell_y;
ellip_z = ell_z;
[gforce,jacobian,laplacian,U] = potential_mat(rp);
Ux = Ux + gforce(1);
Uy = Uy + gforce(2);
Uz = Uz + gforce(3);

xprime = [ x(2);
2*omega*x(4) + omega^2*x(1) + Ux;
x(4);
-2*omega*x(2) + omega^2*x(3) + Uy;
x(6);
Uz];