Code covered by the BSD License  

Highlights from
Rise and Set of the Sun, Moon and Planets

Rise and Set of the Sun, Moon and Planets

by

 

05 Dec 2012 (Updated )

Topocentric rise and set of the Sun, Moon and planets. Source ephemeris is DE421 with NOVAS routines

precess (tjd1, pos1, tjd2)
function pos2 = precess (tjd1, pos1, tjd2)

% this function precesses equatorial rectangular coordinates from
% one epoch to another.  the coordinates are referred to the mean
% equator and equinox of the two respective epochs. see pages 30-34
% of the explanatory supplement to the ae, lieske, et al. (1977)
% astronomy and astrophysics 58, 1-16, and lieske (1979) astronomy
% and astrophysics 73, 282-284.

% input

%  tjd1 = tdb julian date of first epoch

%  pos1 = position vector, geocentric equatorial rectangular
%         coordinates, referred to mean equator and equinox of
%         first epoch
%  tjd2 = tdb julian date of second epoch

% output

%  pos2 = position vector, geocentric equatorial rectangular
%         coordinates, referred to mean equator and equinox of
%         second epoch

% Celestial Computing with MATLAB

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

seccon = 206264.8062470964;

% t0 and t below correspond to Lieske's big t and little t

t0 = (tjd1 - 2451545) / 36525;

t = (tjd2 - tjd1) / 36525;

t02 = t0 * t0;

t2 = t * t;

t3 = t2 * t;

% zeta0, zee, and theta below correspond to Lieske's zeta-sub-a,
% z-sub-a, and theta-sub-a

zeta0 = (2306.2181 + 1.39656 * t0 - 0.000139*t02) * t ...
        + (0.30188 - 0.000344 * t0) * t2 +  0.017998 * t3;

zee = (2306.2181 + 1.39656 * t0 - 0.000139 * t02) * t ...
      + (1.09468 + 0.000066 * t0) * t2 + 0.018203 * t3;

theta = (2004.3109 - 0.85330 * t0 - 0.000217 * t02) * t ...
        + (-0.42665 - 0.000217 * t0) * t2 - 0.041833 * t3;

zeta0 = zeta0 / seccon;

zee = zee / seccon;

theta = theta / seccon;

czeta0 = cos(zeta0);

szeta0 = sin(zeta0);

czee = cos(zee);

szee = sin(zee);

ctheta = cos(theta);

stheta = sin(theta);

% precession rotation matrix follows

xxp = czeta0 * ctheta * czee - szeta0 * szee;
yxp = -szeta0 * ctheta * czee - czeta0 * szee;
zxp = -stheta * czee;

xyp = czeta0 * ctheta * szee + szeta0 * czee;
yyp = -szeta0 * ctheta * szee + czeta0 * czee;
zyp = -stheta * szee;

xzp = czeta0 * stheta;
yzp = -szeta0 * stheta;
zzp = ctheta;

% perform rotation

pos2(1) = xxp * pos1(1) + yxp * pos1(2) + zxp * pos1(3);
pos2(2) = xyp * pos1(1) + yyp * pos1(2) + zyp * pos1(3);
pos2(3) = xzp * pos1(1) + yzp * pos1(2) + zzp * pos1(3);

Contact us