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.

etilt1 (tjd)
```function [oblm, oblt, eqeq, dpsi, deps] = etilt1 (tjd)

% this function computes quantities related to the orientation
% of the earth's rotation axis at julian date tjd

% nutation parameters obtained from jpleph

% input

%  tjd = tdb julian date for orientation parameters

% output

%  oblm = mean obliquity of the ecliptic at date tjd (degrees)
%  oblt = true obliquity of the ecliptic at date tjd (degrees)
%  eqeq = equation of the equinoxes at date tjd (arc seconds)
%  dpsi = nutation in longitude at date tjd (arc seconds)
%  deps = nutation in obliquity at date tjd (arc seconds)

% Celestial Computing with MATLAB

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

t0 = 2451545;

seccon = 206264.8062470964;

t = (tjd - t0) / 36525;

t2 = t * t;

t3 = t2 * t;

% obtain nutation parameters in arcseconds

[x, x, x, x, omega] = funarg(t);

rrd = jplephem(tjd, 14, 0);

psi = rrd(1) * seccon;

eps = rrd(2) * seccon;

% compute mean obliquity of the ecliptic in arcseconds

obm = 84381.4480 - 46.8150 * t - 0.00059 * t2 + 0.001813 * t3;

% compute true obliquity of the ecliptic in arcseconds

obt = obm + eps;

% compute equation of the equinoxes in arcseconds, time seconds
% (iau 1994 and iers 1996 definition)

ee = psi * cos(obm / seccon) + 0.00264 * sin(omega) ...
+ 0.000063 * sin(2 * omega);

ee = ee / 15;

% convert obliquity values to degrees

obm = obm / 3600;

obt = obt / 3600;

oblm = obm;
oblt = obt;
eqeq = ee;
dpsi = psi;
deps = eps;

```