Code covered by the BSD License  

Highlights from
Rise and Set of the Sun, Moon and Planets

Rise and Set of the Sun, Moon and Planets

by

 

05 Dec 2012 (Updated )

Topocentric rise and set of the Sun, Moon and planets. Source ephemeris is DE421 with NOVAS routines

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

% predict minimization/root-finding 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-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 celestial events

% input

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

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

global jdatei jdate rootflg trr

% define root-bracketing and root-finding control parameters

factor = 0.25;            % geometric acceleration factor

dxmax = 0.1;              % rectification interval

rtol = 0.00000001;        % convergence tolerance

% compute and display event start conditions

rootflg = 1;

t1in = topt;

t2in = t1in - 0.05;

[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

rootflg = 0;

froot = feval(objfunc, t1in);

jdate = jdatei + t1in;

feval(prtfunc, 2, jdate)

% compute and display event end conditions

rootflg = 1;

t2in = t1in + 0.05;

[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

trr = troot;

jdate = jdatei + troot;

feval(prtfunc, 3, jdate);

rootflg = 0;

Contact us