function output = SmootingLimitations(NormalizedDataMatrix);
% Input:
%    NormalizedDataMatrix: Each column is a Normalized Power output times
%                          series at one location, with a wind farm power
%                          capacity of one.
% Output:
%    Spectrum + std, 5th and 95th percentiles of the hourly step change
%    function. Spectral plots are also visualized. 
% Assumption:
%    The wind farm capacity is equally distributed over all considered wind
%    farms.


%% ========================================================================
%                  Aggregated spectra (takes a while to run)
% =========================================================================
nts = length(NormalizedDataMatrix(1,:)); % Number of total sites
ndp = length(NormalizedDataMatrix(:,1)); % Number of data points per site
[s,f] = AverageSpectrum(NormalizedDataMatrix(:,1),3600,256);
f_ind = find(f>1/3/3600);
ind5  = round((ndp-1)*0.05);
ind95 = round((ndp-1)*0.95);

tic
output.info = '{n wind farms}(DifferentCombinations,SiteIndex)';
A = 1:nts;
for n = 1:length(A) % number of aggregated time series
    ind = nchoosek(A,n);
    for i = 1:length(ind(:,1));
        Ind = ind(i,:);
        lump    = mean(NormalizedDataMatrix(:,Ind),2);
        [s,f] = AverageSpectrum(lump,3600,256);
        output.spectrum{n}(i,:) = s;
        output.trapz{n}(i) = trapz(f(f_ind),s(f_ind));
        dp = sort(diff(lump));
        output.PCTL5{n}(i)  = dp(ind5);
        output.PCTL95{n}(i) = dp(ind95);
        output.std{n}(i) = std(dp);
    end
    [n,toc]
end
output.f = f;
tic

%% Visualization
Col_map     = hsv(nts+3);
Col_map     = Col_map(1:nts,:);
caxis_range = [1,nts];
for i = A;
    if i == 1;figure();else;hold on;end
    temp = output.spectrum{i};
    s=mean(output.spectrum{i},1);
    loglog(f(2:end),f(2:end).*s(2:end),'-','LineWidth',1,'color',Col_map(i,:));
end

colormap(Col_map);
c=colorbar;
c.Title.String='{   }n';
caxis([0,nts]);
c.Ticks = [1:nts]-0.5;
c.TickLabels = num2str([c.Ticks+0.5]');

xlim([2e-6,2e-4]);
ylim([9e-5,2e-2]);
xx = xlim;
yy = ylim;
tekstur = {'3 days   ','   1 day','2 hours{  }','  2 min'};
T = [1/(3*24*60*60),1/(1*24*60*60),1/(2*3600),1/(2*60)];
index = find(T>=xx(1) & T<=xx(2));
tekstur = tekstur(index);
T = T(index);
for i = 1:length(tekstur)
    hold on,plot(T([i,i]),yy,'k:');
    text(T(i),yy(2),tekstur{i},'Horizontalalignment','center','VerticalAlignment','bottom','FontSize',9)
end
xlabel('f [Hz]')
ylabel('f.S')

%% Plot all compinations of aggregated spectra for selected number of sites
n = [1,3,10,14];
tel = 0;
for i = n;
    tel = tel+1;
    if tel == 1;figure();else;hold on;end
    temp = output.spectrum{i}(:,2:end);
    loglog(f(2:end),repmat(f(2:end),length(temp(:,1)),1).*output.spectrum{i}(:,2:end),'-','LineWidth',1,'color',Col_map(i,:));
end

colormap(Col_map);
c=colorbar;
c.Title.String='{   }n';
caxis([0,nts]);
c.Ticks = [1:nts]-0.5;
c.TickLabels = num2str([c.Ticks+0.5]');

xticks(10.^[-100:100])
yticks(10.^[-100:100])
xlim([2e-6,2e-4]);
ylim([9e-5,2e-2]);
xx = xlim;
yy = ylim;
tekstur = {'3 days   ','   1 day','2 hours{  }','  2 min'};
T = [1/(3*24*60*60),1/(1*24*60*60),1/(2*3600),1/(2*60)];
index = find(T>=xx(1) & T<=xx(2));
tekstur = tekstur(index);
T = T(index);
for i = 1:length(tekstur)
    hold on,plot(T([i,i]),yy,'k:');
    text(T(i),yy(2),tekstur{i},'Horizontalalignment','center','VerticalAlignment','bottom','FontSize',9)
end
xlabel('f [Hz]')
ylabel('f.S')

