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

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