function [Julianday, Juneday, Everday, Everweek, Week]=daynumber(year, month, day, firstyear)
%Usage:
%[Julianday, Juneday, Everday, Everweek, Week]=daynumber(year, month, day, firstyear)
%or
%[Julianday, Juneday, Everday, Everweek, Week]=daynumber()
%or
%[Month, Week, Day]=daynumber(Julianday,NoDisplay)
%or
%[Month, Week, Day]=daynumber(Julianday)
%
%Input with 3 arguments
%year  [n]       year
%month [n]       month
%day   [n]       day
%firstyear [1]   [1950] year to start cumulating daynumbers
%
%Output
%Julianday      Julian day number
%Juneday        day number from 1 June
%Everday        day number since firstyear or 1 Jan 1950 (default)
%Everweek       week number since firstyear or 1 Jan 1950 (default)
%Week           week of the year
%
%Tom Schut, 2005

if nargin < 4 firstyear=1950; end;

md  = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31];
mds = [31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31];
jmd  = [md(:,6:12)  md(:,1:5) ];
jmds = [mds(:,6:12) mds(:,1:5)];

if nargin <3
    days=0;
    for i=1:length(md)
        for j=1:md(i)       
            days=days+1;
            Year(days,1)=1;
            month(days,1)=i;
            day(days,1)=j;
        end
    end
else
    Year=year;
end
n=length(Year);



mj=(month > 5) .* (month-5) + (month <= 5) .* (month + 7);

for i=1:length(month)
    smd(i)=sum(md(:,1:month(i)-1));
    smds(i)=sum(mds(:,1:month(i)-1));
    sjmd(i)=sum(jmd(:,1:mj(i)-1));
    sjmds(i)=sum(jmds(:,1:mj(i)-1));
end

%errors occur in weatherfiles in 1900
a=find(Year==1900);
if length(a == 365)
    sy=(mod(Year,4) ==0 & Year~=1900);
     y=(mod(Year,4) ~=0 | Year==1900);
else
    sy=(mod(Year,4) ==0);
     y=(mod(Year,4) ~=0);
end

if size(y)~=size(smd)
    smd=smd';
    smds=smds';
    sjmd=sjmd';
    sjmds=sjmds';
end

Julianday = y .* smd + sy .* smds;
Juneday = y .* sjmd + sy .* sjmds;

Julianday = Julianday + day;
Juneday = Juneday + day;

for i=1:n
    yrs=firstyear:Year(i)-1;
    sy=(mod(yrs,4) ==0) .* 366;
     y=(mod(yrs,4) ~=0) .* 365;
    Everday(i) = sum(sy) + sum(y) + Julianday(i);
end
if size(Year)~=size(Everday) Everday=Everday';end
    
Everweek=floor(Everday/7) + 1;
Week=floor(Julianday/7) + 1;

if nargin < 1
    disp('Year month day Julianday Juneday Week');
    disp([Year month day Julianday Juneday Week]);
    Julianday=[];
    Juneday=[];
elseif nargin < 2
    jday=Julianday;
    Julianday=[];Juneday=[];Everday=[];Everweek=[];
    disp(' Month  Week   Day');
    for i=1:length(year)
        a=find(year(i)==jday);
        Julianday(i)=month(a);
        Juneday(i)=Week(a);
        Everday(i)=day(a);
        disp([month(a) Week(a) day(a)]);
    end   
elseif nargin < 3
    jday=Julianday;
    Julianday=[];Juneday=[];Everday=[];Everweek=[];
    for i=1:length(year)
        a=find(year(i)==jday);
        Julianday(i)=month(a);
        Juneday(i)=Week(a);
        Everday(i)=day(a);
    end   
end
