%% eshadow.m determines if the inputted points are in the Earth's shadow or not.
%% Assume this "extends" out 250km above the surface of the Earth to represent 
%% absorption by the atmosphere.
%%

function pois=eshadow(day,month,year,timl,latmat,dismat)

% day = day of interest
% month = month of interest
% year = year of interest
% timl = local geographic time (fractional hours)
% latmat = geographic latitude (rad)
% dismat = radial distance (km)
%
% pois = point of interest shadow ( 1 = in shadow, 0 = not in shadow)

re=6371;                           %Earth radii (km)
rb=re+250;                         %radial extension of shadow (km)

%Rotate the coordinates of the points of interest to include the solar declination.
%Use Euler angles, with the single rotation taken to be around the x axis
%(p147 Goldstein 2nd Ed Classical Mech.)
sdec=sundec(day,month,year);                           %declination of Sun (rads)
[xco,yco,zco]=sph2cart(timl*pi/12,latmat,dismat);
yc=yco;
xc=xco.*cos(-sdec)+zco.*sin(-sdec);
zc=-xco.*sin(-sdec)+zco.*cos(-sdec);

%Determine if the point of interest is in the Earth's shadow
pois=abs(yc)<rb & abs(zc)<rb & xc>0;         
