Code covered by the BSD License

# A MATLAB Script for Predicting Transits of Mercury and Venus

### David Eagle (view profile)

07 Dec 2012 (Updated )

Local circumstances of solar transits of the planets Mercury and Venus.

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);
```