function gst = gast_md2000(md2000)
% gst = gast_JD2000(md2000);
%
% Greenwich apparent sidereal time
% input
%  md2000 = Modified Julian Day 2000
% output
%  gst = Greenwich apparent sidereal time (radians)
%        (0 <= gst <= 2 pi)
%
% Celestial Computing with MATLAB

% conversion factors
dtr = pi/180;
atr = dtr/3600;

% time arguments
t = (md2000 - 0.5) / 36525;
t2 = t.* t;
t3 = t.* t2;

% fundamental trig arguments

l = modulo(dtr * (280.4665 + 36000.7698 * t));

lp = modulo(dtr * (218.3165 + 481267.8813 * t));

lraan = modulo(dtr * (125.04452 - 1934.136261 * t));

% nutations in longitude and obliquity

dpsi = atr * (-17.2 * sin(lraan) - 1.32 * sin(2 * l) ...
       - 0.23 * sin(2 * lp) + 0.21 * sin(2 * lraan));

deps = atr * (9.2 * cos(lraan) + 0.57 * cos(2 * l) ...
       + 0.1 * cos(2 * lp) - 0.09 * cos(2 * lraan));

% mean and apparent obliquity of the ecliptic

eps0 = modulo(dtr * (23 + 26 / 60 + 21.448 / 3600) ... 
       + atr * (-46.815 * t - 0.00059 * t2 + 0.001813 * t3)); 

obliq = eps0 + deps;

% greenwich mean and apparent sidereal time

%gstm = modulo(dtr * (280.46061837 + 360.98564736629 * (jdate - 2451545) ...
%       + 0.000387933 * t2 - t3 / 38710000));
gstm = modulo(dtr * (280.46061837 + 360.98564736629 * (md2000 - 0.5) ...
       + 0.000387933 * t2 - t3 / 38710000));
    
gst = modulo(gstm + dpsi.*cos(obliq));

function y = modulo (x)
% modulo 2 pi function
% input
%  x = argument (radians)
% output
%  y = x modulo 2 pi
% Celestial Computing with MATLAB

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

y = mod(x, 2*pi);


