function [] = CoherencyWRFplot(dataMatrix,lon,lat)

% lon: longitude in UTM [m]
% lat: longitude in UTM [m]
% dataMatrix: each column represents power output data at one location

%% Colormap:
colorcodedistance=0:0.1:100;
xxx = 150; % I'm cropping the colormap
colorcode = hsv(length(colorcodedistance)+xxx);
colorcode = colorcode(1:end-xxx,:);

range_f{1} = [1/(3*3600),1/(2*3600)]; % [1/3h, 1/2h]
range_f{2} = [1/(9*3600),1/(6*3600)]; % [1/9h, 1/6h]
range_f{3} = [1/(24*3600),1/(16*3600)]; % [1/24h, 1/16h]
range_f{4} = [1/(3*24*3600),1/(2*24*3600)]; % [1/3days, 1/2days]

A = 1:length(dataMatrix(1,:));
n = 2; % For each calculation, we calculate the coherence between two sites
ind = nchoosek(A,n);
for i = 1:length(ind(:,1));
    x = dataMatrix(:,ind(i,1));
    y = dataMatrix(:,ind(i,2));
    [coh,f]       = Coherence(x,y,3600,256);
    coh_dist(i,1) = sqrt((lon(ind(i,1))-lon(ind(i,2)))^2+(lat(ind(i,1))-lat(ind(i,2)))^2)/1000; % [m] -> [km]
    colorcodeind = find(abs(colorcodedistance-coh_dist(i))==min(abs(colorcodedistance-coh_dist(i))));
    if i == 1;figure();else; hold on;end
    semilogx(f,coh,'color',colorcode(colorcodeind,:))
    % ---------------------------------------------------
    %   Average coherence between selected frequencies
    % ---------------------------------------------------
    for j = 1:length(range_f);
        ind_frange = find(f>=range_f{j}(1) & f<=range_f{j}(2));
        coh_mean{j}(i) = mean(coh(ind_frange));
    end
end
%  Decorate coherence plot:
c=colorbar;
colormap(colorcode);
caxis([min(colorcodedistance),max(colorcodedistance)])
c.YLabel.String = 'Inter-site distance [km]';
xlabel('f [Hz]')
ylabel('\gamma^2')



%% ---------------------------------------------------
%  Average coherence between 1/3h and 1/2h (range_f)
% ---------------------------------------------------
figure()
for j = 1:length(range_f)
    hold on;plot(coh_dist,coh_mean{j},'.','MarkerSize',10);
    if j == 1 | j == 2;
        leg{j} = ['\,\,\,$(',num2str(1/(3600*range_f{j}(1)),'%2.0f'),'\,\textrm{h})^{-1} \leq$ f $\leq (',num2str(1/(3600*range_f{j}(2))),'\,\textrm{h})^{-1}$'];
    else
        leg{j} = ['$(',num2str(1/(3600*range_f{j}(1)),'%2.0f'),'\,\textrm{h})^{-1} \leq$ f $\leq (',num2str(1/(3600*range_f{j}(2))),'\,\textrm{h})^{-1}$'];
    end    
end
xlabel('Inter-site distance [km]')
ylabel('$\overline{\gamma^2}$','Interpreter','latex')
% Legend:
Leg=legend(leg,'Interpreter','latex');
h=gca;
Leg.Location='SouthOutside';
box on


