Code covered by the BSD License  

Highlights from
A MATLAB Script for Predicting Lunar Occultations

A MATLAB Script for Predicting Lunar Occultations

by

 

07 Dec 2012 (Updated )

Predict the local circumstances of lunar occultations of a planet or star.

applan1 (tjd, ujd, l, n, topo, glon, glat, ht)
function [ra, dec, dis] = applan1 (tjd, ujd, l, n, topo, glon, glat, ht)

% this function computes the apparent geocentric or topocentric place
% of a planet or other solar system body. rectangular coordinates of
% solar system bodies are obtained from function solsys.

% jpleph is the source ephemeris for this routine

% input

%  tjd  = tdt julian date for apparent geocentric place
%  ujd  = ut1 julian date for apparent topocentric place
%  l    = body identification number for desired planet
%  n    = body identification number for the earth
%  topo = type of apparent place calculation
%       = 0 ==> geocentric
%       = 1 ==> topocentric
%  glon = geodetic longitude of observer (east +, degrees)
%  glat = geodetic latitude of observer (north +, degrees)
%  ht   = height of observer (meters)

% output

%  ra  = apparent geocentric or topocentric right ascension,
%        referred to true equator and equinox of date (hours)
%  dec = apparent geocentric or topocentric declination,
%        referred to true equator and equinox of date (degrees)
%  dis = true distance from earth to planet (astronomical units)

% Celestial Computing with MATLAB

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

c = 173.1446333;

t0 = 2451545;

% compute t1, the tdb julian date corresponding to tjd

[x, secdif] = tdtimes (tjd);

t1 = tjd + secdif / 86400;

% get position and velocity of the earth wrt barycenter
% of solar system and wrt center of sun

[peb, veb, ierr] = solsys (t1, n, 0);

if (ierr ~= 0)
   ra = 0;
   dec = 0;
   dis = 0;
   return;
end

[pes, ves, ierr] = solsys (t1, n, 1);

if (ierr ~= 0)
   ra = 0;
   dec = 0;
   dis = 0;
   return;
end

for j = 1:1:3
    pb(j) = peb(j);
    vb(j) = veb(j);
    ps(j) = pes(j);
    vs(j) = ves(j);
end

lplan = l;

if (topo == 1)
   % get position and velocity of observer wrt center of earth

   st = gast2 (ujd, 0, 0);

   [x1, x2, eqeq, x3, x4] = etilt1 (t1);

   gast = st + eqeq / 3600;

   [pos1, vel1] = terra (glon, glat, ht, gast);

   pos2 = nutate1 (-t1, pos1);

   pog = precess (t1, pos2, t0);

   vel2 = nutate1 (-t1, vel1);

   vog = precess (t1, vel2, t0);

   % compute position and velocity of observer wrt barycenter of
   % solar system and wrt center of sun

   for j = 1:1:3
       pb(j) = peb(j) + pog(j);
       vb(j) = veb(j) + vog(j);
       ps(j) = pes(j) + pog(j);
       vs(j) = ves(j) + vog(j);
    end
end
 
% get position of planet wrt barycenter of solar system

[pos1, vel1, ierr] = solsys (t1, lplan, 0);

if (ierr ~= 0)
   ra = 0;
   dec = 0;
   dis = 0;
   return;
end

[pos2, tlight] = geocen (pos1, pb);

s = tlight * c;

t2 = t1 - tlight;

while(1)
  [pos1, vel1, ierr] = solsys (t2, lplan, 0);

  if (ierr ~= 0)
     ra = 0;
     dec = 0;
     dis = 0;
     return;
  end

  [pos2, tlight] = geocen (pos1, pb);

  t3 = t1 - tlight;

  if (abs(t3 - t2) < 1e-8)
     pos3 = sunfld (pos2, ps);

     pos4 = aberat (pos3, vb, tlight);

     pos5 = precess (t0, pos4, t1);

     pos6 = nutate1 (t1, pos5);

     [r, d] = angles (pos6);

     ra = r;

     dec = d;

     dis = s;

     break;
  end  

  t2 = t3;
end



Contact us