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.

applan1 (tjd, ujd, l, n, topo, glon, glat, ht)
```function [ra, dec, dis] = applan1 (tjd, ujd, l, n, topo, glon, glat, ht)

% this function computes the apparent geocentric or topocentric place
% of a planet or other solar system body. rectangular coordinates of
% solar system bodies are obtained from function solsys.

% jpleph is the source ephemeris for this routine

% input

%  tjd  = tdt julian date for apparent geocentric place
%  ujd  = ut1 julian date for apparent topocentric place
%  l    = body identification number for desired planet
%  n    = body identification number for the earth
%  topo = type of apparent place calculation
%       = 0 ==> geocentric
%       = 1 ==> topocentric
%  glon = geodetic longitude of observer (east +, degrees)
%  glat = geodetic latitude of observer (north +, degrees)
%  ht   = height of observer (meters)

% output

%  ra  = apparent geocentric or topocentric right ascension,
%        referred to true equator and equinox of date (hours)
%  dec = apparent geocentric or topocentric declination,
%        referred to true equator and equinox of date (degrees)
%  dis = true distance from earth to planet (astronomical units)

% Celestial Computing with MATLAB

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

c = 173.1446333;

t0 = 2451545;

% compute t1, the tdb julian date corresponding to tjd

[x, secdif] = tdtimes (tjd);

t1 = tjd + secdif / 86400;

% get position and velocity of the earth wrt barycenter
% of solar system and wrt center of sun

[peb, veb, ierr] = solsys (t1, n, 0);

if (ierr ~= 0)
ra = 0;
dec = 0;
dis = 0;
return;
end

[pes, ves, ierr] = solsys (t1, n, 1);

if (ierr ~= 0)
ra = 0;
dec = 0;
dis = 0;
return;
end

for j = 1:1:3
pb(j) = peb(j);
vb(j) = veb(j);
ps(j) = pes(j);
vs(j) = ves(j);
end

lplan = l;

if (topo == 1)
% get position and velocity of observer wrt center of earth

st = gast2 (ujd, 0, 0);

[x1, x2, eqeq, x3, x4] = etilt1 (t1);

gast = st + eqeq / 3600;

[pos1, vel1] = terra (glon, glat, ht, gast);

pos2 = nutate1 (-t1, pos1);

pog = precess (t1, pos2, t0);

vel2 = nutate1 (-t1, vel1);

vog = precess (t1, vel2, t0);

% compute position and velocity of observer wrt barycenter of
% solar system and wrt center of sun

for j = 1:1:3
pb(j) = peb(j) + pog(j);
vb(j) = veb(j) + vog(j);
ps(j) = pes(j) + pog(j);
vs(j) = ves(j) + vog(j);
end
end

% get position of planet wrt barycenter of solar system

[pos1, vel1, ierr] = solsys (t1, lplan, 0);

if (ierr ~= 0)
ra = 0;
dec = 0;
dis = 0;
return;
end

[pos2, tlight] = geocen (pos1, pb);

s = tlight * c;

t2 = t1 - tlight;

while(1)
[pos1, vel1, ierr] = solsys (t2, lplan, 0);

if (ierr ~= 0)
ra = 0;
dec = 0;
dis = 0;
return;
end

[pos2, tlight] = geocen (pos1, pb);

t3 = t1 - tlight;

if (abs(t3 - t2) < 1e-8)
pos3 = sunfld (pos2, ps);

pos4 = aberat (pos3, vb, tlight);

pos5 = precess (t0, pos4, t1);

pos6 = nutate1 (t1, pos5);

[r, d] = angles (pos6);

ra = r;

dec = d;

dis = s;

break;
end

t2 = t3;
end

```