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.

occfunc(x)
function fx = occfunc(x)

% occultation 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 dmoon obslat glon glat ht stlocl

global ram decm pmra pmdec parlax radvel

global sdia azimuth elevation

% current utc julian date

jdutc = jdatei + x;

% tdt julian date

jdtdt = utc2tdt(jdutc);

% calculate topocentric unit position vector of the moon

[ramoon, decmoon, rmmoon] = applan1(jdtdt, jdutc, 11, 3, 1, glon, glat, ht);

umoon(1) = cos(htr * ramoon) * cos(dtr * decmoon);
umoon(2) = sin(htr * ramoon) * cos(dtr * decmoon);
umoon(3) = sin(dtr * decmoon);

% compute semidiameter of the moon

sdmoon = asin(dmoon / rmmoon);

% calculate topocentric coordinates of the selected body

if (itarget == 10)
    
   % star
   
   [rastar, decstar] = apstar1 (jdtdt, jdutc, 3, 1, glon, glat, ht, ...
                                ram, decm, pmra, pmdec, parlax, radvel);
                            
else
    
   % planet
   
   [raplanet, decplanet, rmplanet] = applan1(jdtdt, jdutc, itarget, ...
                                             3, 1, glon, glat, ht);
                                         
end

% calculate topocentric unit pointing vector to target planet or star

if (itarget == 10)
    
   % star
   
   utarget(1) = cos(htr * rastar) * cos(dtr * decstar);
   utarget(2) = sin(htr * rastar) * cos(dtr * decstar);
   utarget(3) = sin(dtr * decstar);

   sdbody = 0;
   
else
    
   % 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;
   
end

% apparent hour angle

hangle = stlocl - (htr * ramoon);

% topocentric elevation of the moon (degrees)

elevation = rtd * asin(sin(obslat) * sin(dtr * decmoon) ...
            + cos(obslat) * cos(dtr * decmoon) * cos(hangle));
            
% topocentric azimuth of the moon (degrees)

tmp1 = sin(hangle);

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

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

% compute moon/object separation angle

cpsi = dot(umoon, utarget);

% compute objective function

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

Contact us