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

rsfunc(x)
function fx = rsfunc(x)

% topocentric elevation angle objective function

% input

%  x = elapsed simulation time (days)

% output

%  fx = objective function at x

% Celestial Computing with MATLAB

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

global jdatei itarg icent rootflg

global obslat glon glat ht stlocl

global aunit elevation azimuth

rtd = 180 / pi;

dtr = pi / 180;

jdutc = jdatei + x;

% compute current tdt julian date

jdtdt = utc2tdt(jdutc);

% compute apparent topocentric coordinates of the celestial body

[ra, dec, dis] = applan1(jdtdt, jdutc, itarg, icent, 1, glon, glat, ht);

% apparent hour angle

hangle = stlocl - (2.0 * pi * ra / 24.0);

% topocentric elevation (degrees)

elevation = rtd * asin(sin(obslat) * sin(dtr * dec) ...
            + cos(obslat) * cos(dtr * dec) * cos(hangle));
            
% topocentric azimuth (degrees)

tmp1 = sin(hangle);

tmp2 = cos(hangle) * sin(obslat) - tan(dtr * dec) * cos(obslat);

azimuth = mod(180 + rtd * atan3(tmp1, tmp2), 360);

if (rootflg == 1)
    
   % correct for horizontal refraction
   
   dref = (34 / 60);

   if (itarg == 11)
       
      % correct for semidiameter of the moon
      
      tsdia = rtd * asin((1738 / aunit) / dis);
      
   elseif (itarg == 10)
       
      % correct for semidiameter of the sun
      
      tsdia = rtd * asin((0.5 * 696000 / aunit) / dis);
      
   else
       
      % no semidiameter correction for planets
      
      tsdia = 0;
      
   end
   
else
    
   dref = 0;
   
   tsdia = 0;
   
end

% evaluate objective function

fx = -(elevation + tsdia + dref);

Contact us