function icce
clear;close all

data = read_nc_data;

intervall = 3600; %seconds
delta_t = 1/24/60/60*intervall;
salti_nr = [1 4];

clear epoch int_par U_threshold

t = 1; go =1;
while go
    if (data.time_begin(1) + (t-1)*delta_t) >data.time_end(6)
        break
    end
    t
    time_index = data.time_num>data.time_begin(1) + (t-1)*delta_t & ...
        data.time_num<data.time_begin(1)+t*delta_t;
    
    
    wind_int = data.wind1(time_index,1);
    
    u_max(t) = max(wind_int);
    u_min(t) = min(wind_int);
    
    %     time_index = data.time_num>data.time_begin(5) & data.time_num<data.time_end(5);
    
    % determine intermittency
    epoch(t) = data.time_begin(1) + t*delta_t-0.5*delta_t;
    epoch(t) = data.time_begin(1)-0.5/24 + t*delta_t-0.5*delta_t;
    
    [U_mean(t),U_std(t),U_mean_conf(t,:),U_std_conf(t,:)] = normfit(wind_int);
    
    % saltiphone nr 1
    zero_ind = data.Salti(time_index,salti_nr(1))==0;
    int_par(1,t) = 1-sum(zero_ind)/length(zero_ind);
    U_threshold(1,t) = icdf('Normal',1-int_par(1,t),U_mean(t),U_std(t));
    sum_counts(1,t) = sum(data.Salti(time_index,salti_nr(1)));
    
    if length(wind_int(logical(1-zero_ind)))>5
        top5 = sort(wind_int(logical(1-zero_ind)));
        Ut_min(1,t) = mean(top5(1:5));
    end
    if length(wind_int(logical(zero_ind)))>5
        top5 = sort(wind_int(logical(zero_ind)));
        Ut_max(1,t) = mean(top5(end-4:end));
    end
    
    % saltiphone nr 2
    zero_ind = data.Salti(time_index,salti_nr(2))==0;
    int_par(2,t) = 1-sum(zero_ind)/length(zero_ind);
    U_threshold(2,t) = icdf('Normal',1-int_par(2,t),U_mean(t),U_std(t));
    sum_counts(2,t) = sum(data.Salti(time_index,salti_nr(2)));
    
    if length(wind_int(logical(1-zero_ind)))>5
        top5 = sort(wind_int(logical(1-zero_ind)));
        Ut_min(2,t) = mean(top5(1:5));
    end
    if length(wind_int(logical(zero_ind)))>5
        top5 = sort(wind_int(logical(zero_ind)));
        Ut_max(2,t) = mean(top5(end-4:end));
    end
    
      if t == 21 | t ==25;
        calc_lin_rel(data,time_index,1);
        
        title([datestr(epoch(t)-0.5*delta_t,'dd-mmm-yyyy HH:MM')...
            ' - ' datestr(epoch(t)+0.5*delta_t,15)])
      end
    
    t=t+1;
    
end
legend(['Data ' datestr(epoch(21),'dd-mmm-yyyy HH:MM')],'Linear Fit','Cubic Fit',...
    ['Data ' datestr(epoch(25),'dd-mmm-yyyy HH:MM')],'Linear Fit','Cubic Fit','location','NW')

%%
figure
subplot(3,1,1)

plot(epoch,U_mean,'r*')
title([datestr(epoch(1),'dd') '-' datestr(epoch(end),'dd mmm yyyy')])

hold all
datetick
ylabel('u [m/s]')
ylim([0 12])

subplot(3,1,2)
plot(epoch,sum_counts(1,:),'g*')
hold all
% plot(epoch,250*U_mean.^3,'r.')
datetick
ylabel('Total counts')
% ylim([0 3e5])

subplot(3,1,3)
tides
box on
set(gcf,'PaperPosition',[0.634517 6.34517 10 10])



function [rc,Vt,P_corr] = calc_lin_rel(data,time_index,plots)
S_nr = 1;



%%
wind_line = round((0:0.25:20)*10)/10;
max_count = ones(length(wind_line),5);
min_count = ones(length(wind_line),5);
mean_count = ones(length(wind_line),5);
for j=1:length(wind_line)
    wind_counts(j) = sum(data.wind1(time_index)==wind_line(j));
    count_line = data.Salti(((data.wind1 == wind_line(j)) & time_index) ,:);
    if ~isempty(count_line)
        %         max_count(i,:) = max(count_line);
        %         min_count(i,:) = min(count_line);
        mean_count(j,:) = mean(count_line);
    else
        %         max_count(i,:) = nan;
        %         min_count(i,:) = nan;
        mean_count(j,:) = nan;
        
    end
end

if plots
    figure
    
    subplot(4,1,1)
    plot(data.wind1(time_index),data.Salti(time_index,S_nr),'.')
    legend('Loc. G','location','NW');%legend('boxoff')
    
    ylabel(['counts'])
    xlim([5 14])
    
    subplot(4,1,2:3)
    xlim([5 14])
    box on
    % Add wind speed bins
    
    hold on
    plot(wind_line,mean_count(:,1:5),'*')
    legend('Loc. G','Loc. E','Loc. C','Loc. J','Loc. D','location','NW')
    %     legend boxoff
    ylim([0 1200])
    ylabel('Avg counts [-]')
    subplot(4,1,4)
    plot(wind_line,wind_counts/sum(wind_counts)*100,'r*')
    xlim([5 14])
    % hline(5)
    % hline(2.5)
    % vline(mean(data.wind1(time_index))-2*std(data.wind1(time_index)))
    ylabel({'Occurence';'[%]'})
    xlabel('Wind speed [m/s]')
    vline(mean(data.wind1(time_index)),'-k')
    
    set(gcf,'PaperPosition',[0.634517 6.34517 10 15])
    
    
end
%% fit linear trends on wind averages.

for S_nr = 1:5
    index = ~isnan(mean_count(:,S_nr)) & ...
        wind_line'>mean(data.wind1(time_index))-2*std(data.wind1(time_index)) & ...
        wind_line'<mean(data.wind1(time_index))+2*std(data.wind1(time_index));
    
    l_fit(S_nr,:) = polyfit(wind_line(index)',...
        mean_count((index),S_nr),1);
    P_corr(S_nr,:) = corr(wind_line(index)',...
        mean_count((index),S_nr));
    if plots
        subplot(4,1,2:3)
        plot([0 15],polyval(l_fit(S_nr,:),[0 15]));
        hold all
    end
end


rc = l_fit(:,1);
Vt = -l_fit(:,2)./l_fit(:,1);

if plots
    if 1
        %%
        figure(56)
        xlim([5 14])
        box on
        % Add wind speed bins
        
        hold all
        ho = plot(wind_line,mean_count(:,1),'*')
        legend('Data Loc. A','Data Loc. B','Data Loc. C','location','NW')
        
        for S_nr = [1]% 2 5]
            index = ~isnan(mean_count(:,S_nr)) & ...
                wind_line'>mean(data.wind1(time_index))-2*std(data.wind1(time_index)) & ...
                wind_line'<mean(data.wind1(time_index))+2*std(data.wind1(time_index));
            
            % apply linear fit
            l_fit(S_nr,:) = polyfit(wind_line(index)',...
                mean_count((index),S_nr),1);
            P_corr(S_nr,:) = corr(wind_line(index)',...
                mean_count((index),S_nr));
            % apply 3d order fit
            nl_fit(S_nr,:) = polyfit(wind_line(index).^3',...
                mean_count((index),S_nr),1);
            P_corr_nl(S_nr,:) = corr(wind_line(index).^3',...
                mean_count((index),S_nr));
            
            
            if plots
                %                 subplot(4,1,2:3)
                plot([0 13.5],polyval(l_fit(S_nr,:),[0 13.5]),'-','linewidth',2,'Color',get(ho,'Color'));
                
                hold all
            end
        end
        
        %         plot 3d power plot (outside loop because of colors)
        hold on
        for S_nr = [1]
            plot([0:0.5:13.5], polyval(nl_fit(S_nr,:),[0:0.5:13.5].^3),':','linewidth',2,'Color',get(ho,'Color'))
            hold all
        end
        
        ylim([0 1200])
        ylabel('Avg counts [-]')
        xlabel('Wind speed [m/s]')
        
        %         legend('Data Loc. A','Data Loc. B','Data Loc. C',...
        %             ['Linear Fit A (r = ' num2str(round(P_corr(1)*100)/100) ')'],...
        %             ['Linear Fit B (r = ' num2str(round(P_corr(2)*100)/100) ')'],...
        %             ['Linear Fit C (r = ' num2str(round(P_corr(5)*100)/100) ')'],...
        %             ['Cubic Fit A (r = ' num2str(round(P_corr_nl(1)*100)/100) ')'],...
        %             ['Cubic Fit B (r = ' num2str(round(P_corr_nl(2)*100)/100) ')'],...
        %             ['Cubic Fit C (r = ' num2str(round(P_corr_nl(5)*100)/100) ')'],...
        %             'location','NW')
        set(gcf,'PaperPosition',[0.634517 6.34517 12 6])
        
    end
    
    
end

