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.

occevent (objfunc, prtfunc, ti, tf, dt, dtsml)
function occevent (objfunc, prtfunc, ti, tf, dt, dtsml)

% predict occultation events 

% input

%  objfunc = objective function
%  prtfunc = display results function
%  ti      = initial simulation time
%  tf      = final simulation time
%  dt      = step size used for bounding minima
%  dtsml   = small step size used to determine whether
%            the function is increasing or decreasing

% Celestial Computing with MATLAB

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

global trr

% initialization

tolm = 1e-4;              % convergence tolerance

iend = 0;
   
% check the initial time for a minimum

fmin1 = feval(objfunc, ti);
   
tmin1 = ti;
   
ftemp = feval(objfunc, ti + dtsml);

df = ftemp - fmin1;

dflft = df;

el = ti;

er = el;

t = ti;
   
% if the slope is positive and the minimum is
% negative, calculate event conditions at the initial time

if (df > 0 & fmin1 < 0)
   events1(objfunc, prtfunc, ti, tf, tmin1);
end

while (1)
    % find where function first starts decreasing

    while (1)
       if (df <= 0)
          break;
       end
      
       t = t + dt;

       if (t >= tf)

          % check final time for a minimum
             
          if (iend == 1)
             y = 1;
             return;
          end

          fmin1 = feval(objfunc, tf);

          ftemp = feval(objfunc, tf - dtsml);

          df = fmin1 - ftemp;

          % set minimum time to final simulation time

          tmin1 = tf;

          if (df < 0)
             % if both the slope and minimum are negative,
             % calculate the event conditions at the final
             % simulation time

             if (fmin1 < 0)
                events1(objfunc, prtfunc, ti, tf, tmin1);
             end

             % otherwise, we're finished
             y = 1;
             return;
          end

          if (dflft > 0)
             y = 1;
             return;
          end

          er = tf;

          iend = 1;
       end

       ft = feval(objfunc, t);

       ftemp = feval(objfunc, t - dtsml);

       df = ft - ftemp;
    end

    % function decreasing - find where function
    % first starts increasing

    while (1)
        el = t;
        dflft = df;
   
        t = t + dt;
   
        if (t >= tf)
      
           % check final time for a minimum

           if (iend == 1)
              y = 1;
              return;
           end

           fmin1 = feval(objfunc, tf);

           ftemp = feval(objfunc, tf - dtsml);

           df = fmin1 - ftemp;

           % set minimum time to final simulation time

           tmin1 = tf;

           if (df < 0)
              % if both the slope and minimum are negative,
              % calculate the event conditions at the final
              % simulation time

              if (fmin1 < 0)
                 events1(objfunc, prtfunc, ti, tf, tmin1);
              end

              % otherwise, we're finished
              y = 1;
              return;
           end

           if (dflft > 0)
              y = 1;
              return;
           end

           er = tf;

           iend = 1;
        end

        ft = feval(objfunc, t);

        ftemp = feval(objfunc, t - dtsml);

        df = ft - ftemp;

        if (df > 0)
           break;
        end

    end

    er = t;

    % calculate minimum using Brent's method
    
    [tmin1, fmin1] = minima(objfunc, el, er, tolm);
    
    el = er;
   
    dflft = df;

    % if the minimum is negative,
    % calculate event conditions for this minimum

    if (fmin1 < 0)
       events1(objfunc, prtfunc, ti, tf, tmin1);
       t = trr;
    end

end

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

function events1 (objfunc, prtfunc, ti, tf, topt)

% compute and display celestial events

% input

%  objfunc = objective function
%  prtfunc = display results function
%  ti      = initial simulation time
%  tf      = final simulation time
%  topt    = extrema time

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

global jdatei jdate trr

% define root-bracketing and root-finding control parameters

factor = 0.25;            % geometric acceleration factor
dxmax = 300 / 86400;      % rectification interval
rtol = 1e-4;              % convergence tolerance

% compute and display event start conditions

t1in = topt;

t2in = t1in - 30 / 86400;

[t1out, t2out] = broot(objfunc, t1in, t2in, factor, dxmax);

[troot, froot] = brent(objfunc, t1out, t2out, rtol);

% set to initial time if before ti

if (troot < ti)
   troot = ti;
   froot = feval(objfunc, ti);
end

jdate = jdatei + troot;

feval(prtfunc, 1, jdate);

% compute and display event end conditions

t2in = t1in + 30 / 86400;

[t1out, t2out] = broot(objfunc, t1in, t2in, factor, dxmax);

[troot, froot] = brent(objfunc, t1out, t2out, rtol);

% set to final time if after tf

if (troot > tf)
   troot = tf;
   froot = feval(objfunc, tf);
end

jdate = jdatei + troot;

feval(prtfunc, 2, jdate);

trr = troot;

Contact us