Code covered by the BSD License  

Highlights from
A MATLAB Script for Predicting Lunar Occultations

A MATLAB Script for Predicting Lunar Occultations

by

 

07 Dec 2012 (Updated )

Predict the local circumstances of lunar occultations of a planet or star.

loccult.m
% loccult.m            December 5, 2012

% lunar occultations of a planet or star

% source ephemeris is DE421

% Celestial Computing with MATLAB

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

clear all;

global rtd dtr htr sdia iephem ephname km

global itarget aunit

global dmoon jdatei obslat glon glat ht

global ram decm pmra pmdec parlax radvel

% 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 moon (au)

dmoon = 1738 / aunit;

% semidiameter of the planets (radians)

sdia(1) = atr * 3.36;
sdia(2) = atr * 8.41;
sdia(3) = atr * 4.68;
sdia(5) = atr * 98.44;
sdia(6) = atr * 82.73;
sdia(7) = atr * 35.02;
sdia(8) = atr * 33.50;
sdia(9) = atr * 2.07;

% define target body name vector

pdata = ['Mercury'; ...
    'Venus  '; ...
    'Mars   '; ...
    'Jupiter'; ...
    'Saturn '; ...
    'Uranus '; ...
    'Neptune'; ...
    'Pluto  '; ...
    'Star   '];

tname = cellstr(pdata);

% begin simulation

clc; home;

fprintf('\n\nlocal circumstances of lunar occultations');
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('\n    target body menu\n');
    
    fprintf('\n  <1> Mercury');
    fprintf('\n  <2> Venus');
    fprintf('\n  <3> Mars');
    fprintf('\n  <4> Jupiter');
    fprintf('\n  <5> Saturn');
    fprintf('\n  <6> Uranus');
    fprintf('\n  <7> Neptune');
    fprintf('\n  <8> Pluto');
    fprintf('\n  <9> Star');
    
    fprintf('\n\nplease select the target body\n');
    
    itarg = input('? ');
    
    if (itarg >= 1 && itarg <= 11)
        break;
    end
end

if (itarg == 9)
    
    [filename, pathname] = uigetfile('*.dat', 'please select an input data file');
    
    [fid, sname, ram, decm, pmra, pmdec, parlax, radvel] = readstar(filename);
    
end

if (itarg >= 3 && itarg <= 9)
    
    itarget = itarg + 1;
    
else
    
    itarget = itarg;
    
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 lunar occultation conditions

occevent('occfunc', 'occprint', ti, tf, dt, dtsml);

fprintf('\n\n');

Contact us