Code covered by the BSD License  

Highlights from
A MATLAB Script for Predicting Transits of Mercury and Venus

A MATLAB Script for Predicting Transits of Mercury and Venus

by

 

07 Dec 2012 (Updated )

Local circumstances of solar transits of the planets Mercury and Venus.

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

% predict transit 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 rmsun rmplanet

if (rmplanet > rmsun)
    
   trr = topt;
   
   return
   
end

% define root-bracketing and root-finding control parameters

factor = 0.25;            % geometric acceleration factor
dxmax = 300/86400;        % rectification interval (days)
rtol = 1e-4;              % convergence tolerance

% compute and display ingress 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 least distance conditions

froot = feval(objfunc, t1in);

jdate = jdatei + t1in;

feval(prtfunc, 2, jdate)

% compute and display egress 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);

trr = troot;

Contact us