Code covered by the BSD License  

Highlights from
A MATLAB Script for Predicting Transits of Mercury and Venus

A MATLAB Script for Predicting Transits of Mercury and Venus

by

 

07 Dec 2012 (Updated )

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

tranfunc(x)
function fx = tranfunc(x)

% transit objective function

% input

%  x = elapsed simulation time (days)

% output

%  fx = objective function at x

% Celestial Computing with MATLAB

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

global rtd dtr htr jdatei itarget

global dsun obslat glon glat ht stlocl

global sdia azimuth elevation rmsun rmplanet

% current utc julian date

jdutc = jdatei + x;

% tdt julian date

jdtdt = utc2tdt(jdutc);

% calculate topocentric unit position vector of the sun

[rasun, decsun, rmsun] = applan1(jdtdt, jdutc, 10, 3, 1, glon, glat, ht);

usun(1) = cos(htr * rasun) * cos(dtr * decsun);
usun(2) = sin(htr * rasun) * cos(dtr * decsun);
usun(3) = sin(dtr * decsun);

% compute semidiameter of the sun

sdsun = asin(dsun / rmsun);

% calculate topocentric coordinates of the planet

[raplanet, decplanet, rmplanet] = applan1(jdtdt, jdutc, itarget, ...
    3, 1, glon, glat, ht);

% calculate topocentric unit pointing vector to target planet

utarget(1) = cos(htr * raplanet) * cos(dtr * decplanet);
utarget(2) = sin(htr * raplanet) * cos(dtr * decplanet);
utarget(3) = sin(dtr * decplanet);

sdbody = sdia(itarget) / rmplanet;

% apparent hour angle

hangle = stlocl - (htr * rasun);

% topocentric elevation of the sun (degrees)

elevation = rtd * asin(sin(obslat) * sin(dtr * decsun) ...
    + cos(obslat) * cos(dtr * decsun) * cos(hangle));

% topocentric azimuth of the sun (degrees)

tmp1 = sin(hangle);

tmp2 = cos(hangle) * sin(obslat) - tan(dtr * decsun) * cos(obslat);

azimuth = mod(180 + rtd * atan3(tmp1, tmp2), 360);

% compute sun/object separation angle

cpsi = dot(usun, utarget);

% compute objective function

fx = acos(cpsi) - (sdsun + sdbody);

Contact us