function [diuCyc,diuStd]=diuCycle(inputData,startDate,finalDate)

% inputData can be an array containing a time series of Data with n time
% steps (rows) and m different variables (coloumns)

% startDate and final have to be in the following format:
% dd-mmm-yyyy HH:MM:SS (e.g. 04-May-1990 10:30:00)

% time units (in seconds)
m = 60;          % one minute in seconds
h = 60*m;   % one hour in seconds
d = 24*h;     % one day in seconds
y = 365*d;     % one year in seconds


dtTotal  = size(inputData,1);
nw       = size(inputData,2);        % number of datasets
nm       = 12;                       % number of months per year


startNum   = datenum(startDate);
finalNum   = datenum(finalDate);
timeData   = linspace(startNum,finalNum,dtTotal)';

dayNum = finalNum - startNum;
dt     = round(dayNum * d)/(dtTotal-1);
dtd    = d/dt;                              % number of time steps per day


% sort data by month
monthCount = month(timeData);
monthData  = cell(nm,nw);




for iw = 1:nw
    for im = 1:nm
        monthData{im,iw} = inputData(monthCount == im,iw);
    end
end   


% number of days of the respective month
nday = cellfun( @(x) length(nanmean(reshape(x,dtd,[]))'), monthData, 'UniformOutput',false );
nday = cell2mat(nday);

% reshape into diurnal cycles
monthData = cellfun( @(x) reshape(x,dtd,[]), monthData, 'UniformOutput',false);

 
ns = 4;                     % number of seasons per year
nh = 2;                     % number of halfs per year

nmps = nm/ns;               % number of months per season
nmph = nm/nh;               % number of months per half year

monSeas = [12 1:11];        % month order of seasons
monHalf = [10:12 1:9];      % month order of half years

seasNames = {'DJF' 'MAM' 'JJA' 'SON'};
halfNames = {'winter','summer'};

seasData = cell(ns,nw);
halfData = cell(nh,nw);


for iw = 1:nw 
    for im = 1:nm
        halfData{ceil(im/nmph),iw} = [halfData{ceil(im/nmph),iw}, monthData{monHalf(im),iw}];
        seasData{ceil(im/nmps),iw} = [seasData{ceil(im/nmps),iw}, monthData{monSeas(im),iw}];
    end   
end

% average by diurnal cycles
seasCyc = cellfun( @(x) nanmean(x')', seasData, 'UniformOutput',false );
halfCyc = cellfun( @(x) nanmean(x')', halfData, 'UniformOutput',false );
seasStd = cellfun( @(x) nanstd(x')',  seasData, 'UniformOutput',false );
halfStd = cellfun( @(x) nanstd(x')',  halfData, 'UniformOutput',false );



for i = 1:length(seasNames)
    diuCyc.(seasNames{i}) = cell2mat(seasCyc(i,:));
    diuStd.(seasNames{i}) = cell2mat(seasStd(i,:));
end

for i = 1:length(halfNames)
    diuCyc.(halfNames{i}) = cell2mat(halfCyc(i,:));
    diuStd.(halfNames{i}) = cell2mat(halfStd(i,:));
end
    




end