caeqm (time, y)

function yddot = caeqm (time, y)
   
% heliocentric equations of motion
% with planetary perturbations

% required by cae2ho.m

% Celestial Computing with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global aunit xmu jdatei iephtype

% current julian date

jdate = jdatei + time;

% distance from sun to object

rs2sc = sqrt(y(1) * y(1) + y(2) * y(2) + y(3) * y(3));
   
rrs2sc = -xmu(1) / rs2sc^3;
   
% calculate planetary position vectors

% Mercury

if (iephtype == 1)
   result = jplephem(jdate, 2, 11);   
else
   [result, ierr] = slp96 (jdate, 1, 11, 1, 1, 1);
end

for i = 1:1:3
    rmercury(i) = result(i);
end

rpm(1) = sqrt(dot(rmercury, rmercury));

for i = 1:1:3
    rp(1, i) = rmercury(i);
end

% Venus

if (iephtype == 1)
   result = jplephem(jdate, 2, 11);
else
   [result, ierr] = slp96 (jdate, 2, 11, 1, 1, 1);
end

for i = 1:1:3
   rvenus(i) = result(i);
end

rpm(2) = sqrt(dot(rvenus, rvenus));

for i = 1:1:3
    rp(2, i) = rvenus(i);
end

% Earth

if (iephtype == 1)
   result = jplephem(jdate, 3, 11);
else
   [result, ierr] = slp96 (jdate, 12, 11, 1, 1, 1);
end

for i = 1:1:3
    rearth(i) = result(i);
end

rpm(3) = sqrt(dot(rearth, rearth));

for i = 1:1:3
    rp(3, i) = rearth(i);
end

% Mars

if (iephtype == 1)
   result = jplephem(jdate, 4, 11);
else
   [result, ierr] = slp96 (jdate, 4, 11, 1, 1, 1);
end

for i = 1:1:3
   rmars(i) = result(i);
end

rpm(4) = sqrt(dot(rmars, rmars));

for i = 1:1:3
    rp(4, i) = rmars(i);
end

% Jupiter

if (iephtype == 1)
   result = jplephem(jdate, 5, 11);
else
   [result, ierr] = slp96 (jdate, 5, 11, 1, 1, 1);   
end

for i = 1:1:3
   rjupiter(i) = result(i);
end

rpm(5) = sqrt(dot(rjupiter, rjupiter));

for i = 1:1:3
    rp(5, i) = rjupiter(i);
end

% Saturn

if (iephtype == 1)
   result = jplephem(jdate, 6, 11);
else
   [result, ierr] = slp96 (jdate, 6, 11, 1, 1, 1);   
end

for i = 1:1:3
   rsaturn(i) = result(i);
end

rpm(6) = sqrt(dot(rsaturn, rsaturn));

for i = 1:1:3
    rp(6, i) = rsaturn(i);
end

for i = 1:1:6
    rrpm(i) = xmu(i + 1) / (rpm(i)^3);
end

% compute planetocentric position vectors of spacecraft

for i = 1:1:6
    rp2sc(i, 1) = y(1) - rp(i, 1);
    rp2sc(i, 2) = y(2) - rp(i, 2);
    rp2sc(i, 3) = y(3) - rp(i, 3);
end

% compute planetocentric distances of spacecraft

for i = 1:1:6
    rp2scm(i) = sqrt(rp2sc(i, 1)^2 + rp2sc(i, 2)^2 + rp2sc(i, 3)^2);
    rrp2scm(i) = xmu(i + 1) / (rp2scm(i)^3);
end

% compute planetary perturbations

for j = 1:1:3
    accp(j) = 0;
    for i = 1:1:6
        accp(j) = accp(j) - rp2sc(i, j) * rrp2scm(i) - rp(i, j) * rrpm(i);
    end
end

% compute integration vector

yddot(1) = y(4);

yddot(2) = y(5);

yddot(3) = y(6);

yddot(4) = accp(1) + y(1) * rrs2sc;
    
yddot(5) = accp(2) + y(2) * rrs2sc;
    
yddot(6) = accp(3) + y(3) * rrs2sc;