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

etilt1 (tjd)
function [oblm, oblt, eqeq, dpsi, deps] = etilt1 (tjd)

% this function computes quantities related to the orientation
% of the earth's rotation axis at julian date tjd

% nutation parameters obtained from jpleph

% input

%  tjd = tdb julian date for orientation parameters

% output

%  oblm = mean obliquity of the ecliptic at date tjd (degrees)
%  oblt = true obliquity of the ecliptic at date tjd (degrees)
%  eqeq = equation of the equinoxes at date tjd (arc seconds)
%  dpsi = nutation in longitude at date tjd (arc seconds)
%  deps = nutation in obliquity at date tjd (arc seconds)

% Celestial Computing with MATLAB

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

t0 = 2451545;

seccon = 206264.8062470964;

t = (tjd - t0) / 36525;

t2 = t * t;

t3 = t2 * t;

% obtain nutation parameters in arcseconds

[x, x, x, x, omega] = funarg(t);

rrd = jplephem(tjd, 14, 0);

psi = rrd(1) * seccon;

eps = rrd(2) * seccon;

% compute mean obliquity of the ecliptic in arcseconds

obm = 84381.4480 - 46.8150 * t - 0.00059 * t2 + 0.001813 * t3;

% compute true obliquity of the ecliptic in arcseconds

obt = obm + eps;

% compute equation of the equinoxes in arcseconds, time seconds
% (iau 1994 and iers 1996 definition)

ee = psi * cos(obm / seccon) + 0.00264 * sin(omega) ...
     + 0.000063 * sin(2 * omega);

ee = ee / 15;

% convert obliquity values to degrees

obm = obm / 3600;

obt = obt / 3600;

oblm = obm;
oblt = obt;
eqeq = ee;
dpsi = psi;
deps = eps;


Contact us