function N = hoursdarkness( J, latitude )
%HOURSDARKNESS Calculating the <hours of darkness> for a given [latitude(deg)] and [day of the year]
% 
% J may be a vector.
% 
% Source: Ecological Modeling, volume 80 (1995) pp. 87-95, "A Model Comparison
% for Daylength as a Function of Latitude and Day of the Year." According to
% the article, the model's error is less than one minute within 40 degrees of
% the equator, and less than seven minutes within 60 degrees and is commonly
% within two minutes for these latitudes.
% 
% The model:
% 
%    N = hours of darkness
%    L = latitude in radians
%    J = day of the year (0 - 365)
%    P = asin[.39795*cos(.2163108 + 2*atan{.9671396*tan[.00860(J-186)]})]
% 
%                        sin(0.8333*pi/180) + sin(L)*sin(P)
%    N = (24/pi)*acos *  -----------------------------------------
%                                 cos(L)*cos(P)
% 
%   note: cos(..), acos(..), etc's arguments must all be in radians
%

degToRad = pi / 180;
L = latitude * degToRad;
P = asin(.39795*cos(.2163108 + 2*atan(.9671396*tan(.00860*(J-186.)))));
numerator = sin(degToRad*0.8333) + sin(L)*sin(P);
denominator = cos(L)*cos(P);
N = (24.0/pi)*acos(numerator./denominator);
end

