Code covered by the BSD License  

Highlights from
A MATLAB Script for Predicting Lunar Eclipses

A MATLAB Script for Predicting Lunar Eclipses

by

 

06 Dec 2012 (Updated )

Predicting the local circumstances of lunar eclipses.

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

% predict lunar eclipse 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

global psi pmoon psun sdsun sdmoon

% define root-bracketing and root-finding control parameters

factor = 0.25;            % geometric acceleration factor

dxmax = 600 / 86400;      % rectification interval

rtol = 1.0e-4;            % root-finding convergence tolerance

% display type of eclipse

if (psi < 1.02 * (0.99834 * pmoon + psun - sdsun) - sdmoon)
    
   fprintf('\n\ntotal lunar eclipse');
   fprintf('\n===================\n');
   
elseif (psi < 1.02 * (0.99834 * pmoon + psun - sdsun) + sdmoon)
    
   fprintf('\n\npartial lunar eclipse \n');
   fprintf('\n=====================\n');
   
elseif (psi < 1.02 * (0.99834 * pmoon + psun + sdsun) + sdmoon)
    
   fprintf('\n\npenumbral lunar eclipse \n');
   fprintf('\n=======================\n');
   
end

% 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 greatest eclipse conditions

jdate = jdatei + topt;

feval(prtfunc, 2, 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, 3, jdate);

% save current value of root

trr = troot;

Contact us