The input day format is somewhat non-user-friendly. The function header describes the input date as the "day of the year, counted starting with the day of the December solstice in the first year of a Great Year." To maintain Travis' functionality and get this function to work with Matlab's datenum format, add this to the beginning of the function:
if nargin>0 && any(Day>366)
[y,~,~] = datevec(Day);
Day = 11+datenum(Day)-datenum(y,1,1);
The above assumes that if any input date exceeds 366, it must be in datenum format. Then the function *approximates* day length assuming a static December solstice. Results will not be exact on the daily timescale, but should be sufficient to put multi-year data in datenum format in context with solar cycles.