Code covered by the BSD License

# A MATLAB Script for Predicting Lunar Eclipses

### David Eagle (view profile)

06 Dec 2012 (Updated )

Predicting the local circumstances of lunar eclipses.

oevent1 (objfunc, prtfunc, ti, tf, dt, dtsml)
```function oevent1 (objfunc, prtfunc, ti, tf, dt, dtsml)

% predict minimization/root-finding orbital events

% time argument in days

% 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

% Orbital Mechanics with Matlab

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

global trr

% initialization

tolm = 1e-8;              % 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 minimization/root-finding orbital events

% input

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

% Orbital Mechanics with Matlab

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

global jdatei jdate trr

% define root-bracketing and root-finding control parameters

factor = 0.25;            % geometric acceleration factor

dxmax = 300 / 86400;      % rectification interval

rtol = 0.00000001;        % convergence tolerance

% compute and display event start conditions

t1in = topt;

t2in = t1in - (10 / 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 conditions at optimum

froot = feval(objfunc, t1in);

jdate = jdatei + t1in;

feval(prtfunc, 2, jdate)

% compute and display event end conditions

t2in = t1in + (10 / 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);

trr = troot;```