Code covered by the BSD License  

Highlights from
Encke's Method

Encke's Method

by

 

Encke's method for solving the orbital initial value problem (IVP).

earth (jdate)
function rearth = earth (jdate)

% true-of-date heliocentric, ecliptic
% position vector of the Earth

% input
    
%  jdate = julian day

% output

%  rearth = position vector of the earth (kilometers)

% Orbital Mechanics with Matlab

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

atr = pi / 648000;

rearth = zeros(3, 1);

% fundamental time arguments

djd = jdate - 2451545;

t = djd / 36525 + 1;

gs = r2r(0.993126 + 0.0027377785 * djd);
lm = r2r(0.606434 + 0.03660110129 * djd);
ls = r2r(0.779072 + 0.00273790931 * djd);
g2 = r2r(0.140023 + 0.00445036173 * djd);
g4 = r2r(0.053856 + 0.00145561327 * djd);
g5 = r2r(0.056531 + 0.00023080893 * djd);

pl = 6910 * sin(gs) + 72 * sin(2 * gs) - 17 * t * sin(gs);
pl = pl - 7 * cos(gs - g5) + 6 * sin(lm - ls);
pl = pl + 5 * sin(4 * gs - 8 * g4 + 3 * g5);
pl = pl - 5 * cos(2 * (gs - g2));
pl = pl - 4 * sin(gs - g2) + 4 * cos(4 * gs - 8 * g4 + 3 * g5);
pl = pl + 3 * sin(2 * (gs - g2)) - 3 * sin(g5) - 3 * sin(2 * (gs - g5));
   
plon = ls + atr * pl;

plat = 0;

r = 149597870.691 * (1.00014 - 0.01675 * cos(gs) - 0.00014 * cos(2 * gs));
   
rearth(1) = -r * cos(plat) * cos(plon);
rearth(2) = -r * cos(plat) * sin(plon);
rearth(3) = -r * sin(plat);

Contact us