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.

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;

Contact us