function [sun_i,r_sun]=solar(yr,mth,day,hr,min,sec,dt,tf);
%function [sun_i,ra_dec]=solar(yr,mth,day,hr,min,sec,dt,tf);
%
% This function computes the solar ephemeris given date and time
%
%  The inputs are:      
%         yr = year, e.g. 1995
%        mth = month, e.g. Jan=1, Feb=2, etc.
%        day = day, e.g. 1-31
%         hr = hour, e.g. 0-23
%        min = minutes, e.g. 0-59
%        sec = seconds, e.g. 0-59
%         dt = sampling interval (sec)
%         tf = run time (sec)
%  The outputs are:
%      sun_i = x,y,z components of ECI unit vector
%      r_sun = magnitude of Sun position (km) 
%

% John L. Crassidis 11/19/12
% Algorithm derived from Vallado's third edition book

t=[0:dt:tf]';
sec=sec+t;
d2r = pi/180;
au=149597870.0;

% Get Julian Date and Centuries
jd=julian(yr,mth,day,hr,min,sec);
jd_cent=(jd-2451545)/36525;

% Mean Longtitude and Other Paramters
lam=280.460+36000.771*jd_cent;
m_sun=357.5277233+35999.050*jd_cent;
lam_ecl=lam+1.914666471*sin(m_sun*d2r)+0.019994643*sin(2*m_sun*d2r);
eps=23.439291-0.0130042*jd_cent;

% Sun Location and Magnitude
sun_i=[cos(lam_ecl*d2r) cos(eps*d2r).*sin(lam_ecl*d2r) sin(eps*d2r).*sin(lam_ecl*d2r)];
r_sun=(1.000140612-0.016708617*cos(m_sun*d2r)-0.000139589*cos(2*m_sun*d2r))*au;
