% calculates the Julian day given year, month and day
function [yday] = md2yd(year, month, day)
day_tab=[31,28,31,30,31,30,31,31,30,31,30,31; 31,29,31,30,31,30,31,31,30,31,30,31];

%leap = (mod(year,4) == 0 && mod(year,100) ~= 0 || mod(year,400) == 0) + 1;
leap = (mod(year,4) == 0 & mod(year,100) ~= 0 | mod(year,400) == 0) + 1;
yday = day;
for i = 1:month-1
    yday = yday + day_tab(leap,i);
end
