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.

transit.m
% transit.m            December 6, 2012     

% transits of Mercury and Venus

% source ephemeris is DE421

% Celestial Computing with MATLAB

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

global rtd dtr htr sdia iephem ephname

global km itarget aunit dsun

global jdatei obslat glon glat ht

% initialize jpl ephemeris

iephem = 1;

ephname = 'de421.bin';

% request output in au and au/day

km = 0;

% read leap seconds data file

readleap;

% conversion factor - astronomical unit to kilometers

aunit = 149597870.691;

% conversion factor - radians to degrees

rtd = 180 / pi;

% conversion factor - degrees to radians

dtr = pi / 180;

% conversion factor - hours to radians

htr = pi / 12;

% conversion factor - arc seconds to radians

atr = dtr / 3600;

% radius of the sun (au)

dsun = 696000 / aunit;

% semidiameter of the planets (radians)

sdia(1) = atr * 3.36;
sdia(2) = atr * 8.41;

% begin simulation

clc; home;

fprintf('\n\ntransits of Mercury and Venus');
fprintf('\n=============================\n');

% request calendar date

fprintf('\nplease input an initial UTC calendar date\n');

[month, day, year] = getdate;
   
jdutc = julian(month, day, year);

% initial UTC Julian date

jdatei = jdutc;

% request search duration

while (1)
    
   fprintf('\nplease input the search duration (days)\n');

   ndays = input('? ');
   
   if (ndays > 0)
      break;
   end  
   
end

% select target body
   
while(1)

   fprintf('\ntarget body menu\n');
       
   fprintf('\n  <1> Mercury\n');
   fprintf('\n  <2> Venus');
   
   fprintf('\n\nplease select the target body\n');
          
   itarget = input('? ');
            
   if (itarget == 1 || itarget == 2)
      break;
   end
   
end
   
% request observer coordinates

[obslat, obslong, obsalt] = getobs;

% convert observer east longitude to degrees

glon = obslong * rtd;

% convert observer geodetic latitude to degrees

glat = obslat * rtd;

% convert observer altitude to meters

ht = 1000 * obsalt;
   
% define search parameters

ti = 0;

tf = ndays;

dt = 0.25;

dtsml = 0.1;

% find transit conditions

traevent('tranfunc', 'tranprt', ti, tf, dt, dtsml);

fprintf('\n\n');

Contact us