function [rasc, decl, r_sun] = sun_md2000(mjd2000)
% [rasc, decl] = sun_md2000(mjd2000);
% solar ephemeris
% input
%  mjd2000 = modified julian day 2000
% output
%  rasc = right ascension of the sun (radians)
%         (0 <= rasc <= 2 pi)
%  decl = declination of the sun (radians)
%         (-pi/2 <= decl <= pi/2)
%  rsun = eci position vector of the sun (km)
% note
%  coordinates are inertial, geocentric,
%  equatorial and true-of-date

% Celestial Computing with MATLAB
% Modified by Nils Olsen, DSRI
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

atr = pi/648000;

mjd2000 = mjd2000(:);

% time arguments
%djd = jdate - 2451545;
djd = mjd2000-0.5;

t = (djd / 36525) + 1;

% fundamental arguments (radians)
gs = r2r(0.993126 + 0.0027377785 * djd);
lm = r2r(0.606434 + 0.03660110129 * djd);
ls = r2r(0.779072 + 0.00273790931 * djd);
g2 = r2r(0.140023 + 0.00445036173 * djd);
g4 = r2r(0.053856 + 0.00145561327 * djd);
g5 = r2r(0.056531 + 0.00023080893 * djd);
rm = r2r(0.347343 - 0.00014709391 * djd);

% geocentric, ecliptic longitude of the sun (radians)

plon = 6910 * sin(gs) + 72 * sin(2 * gs) - 17 * t.*sin(gs);
plon = plon - 7 * cos(gs - g5) + 6 * sin(lm - ls) ... 
       + 5 * sin(4 * gs - 8 * g4 + 3 * g5);
plon = plon - 5 * cos(2 * (gs - g2)) - 4 * (sin(gs - g2) ... 
       - cos(4 * gs - 8 * g4 + 3 * g5));
plon = plon + 3 * (sin(2 * (gs - g2)) - sin(g5) - sin(2 * (gs - g5)));
plon = ls + atr * (plon - 17 * sin(rm));

%mod(plon,2*pi)*180/pi
% geocentric distance of the sun (kilometers)

rsm = 149597870.691 * (1.00014 - 0.01675 * cos(gs) - 0.00014 * cos(2 * gs));

% obliquity of the ecliptic (radians)

obliq = atr * (84428 - 47 * t + 9 * cos(rm));

% geocentric, equatorial right ascension and declination (radians)
a = sin(plon).*cos(obliq);
b = cos(plon);
   
rasc = atan2(a, b);
decl = asin(sin(obliq).*sin(plon));

% geocentric position vector of the sun (kilometers)

r_sun = [rsm.*cos(rasc).*cos(decl)  rsm.*sin(rasc).*cos(decl) rsm.*sin(decl)];


function y = r2r (x)
% revolutions to radians function
% input
%  x = argument (revolutions; 0 <= x <= 1)
% output
%  y = equivalent x (radians; 0 <= y <= 2 pi)

y = 2*pi*mod(x,1);
