%gravity asteroid
function dx=gra(t,x)

global R density

G=6.674*10^(-11);      %gravity constant
M=4/3*pi()*R^3*density*1000;         % mass of asteriod, unit: kg

dx=zeros(6,1);
l=sqrt(x(1)^2+x(3)^2+x(5)^2);
sin_theta=x(1)/l;
cos_theta=sqrt(x(3)^2+x(5)^2)/l;
sin_psi=x(3)/sqrt(x(3)^2+x(5)^2);
cos_psi=x(5)/sqrt(x(3)^2+x(5)^2);
 



dx(1)=x(2);             % vertical velocity
dx(2)=-G*M/l^2*sin_theta; % vertical accleration 
dx(3)=x(4);             % y velocity
dx(4)=-G*M/l^2*cos_theta*sin_psi; % y accleration
dx(5)=x(6);              % x velocity
dx(6)=-G*M/l^2*cos_theta*cos_psi; % x accleration


end