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.

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;

```