path = '/Users/iljagolikov/Documents/ssh/study_Ilja/baka/period/fsum_area_'

pikkused = [];
aastad = 1993:1:2020;
algus= [];
lopp=[];
andmed = [];

for aasta = string(1993:1:2020)

ncf = append(path,aasta,'.nc');

data = ncread(ncf, 'siconc');
data = data(:);

st = find(data,1,'first');
ed = find(data,1,'last');

if sum(data)~=0
    
    algus(end+1) = st;
    lopp(end+1) = ed;
    pikkused(end+1) = ed - st + 1 ;
    andmed = cat(2,andmed,data); 
    
else
    pikkused(end+1) = NaN;     %kui polnud jaad siis pikkus = 0
    data(1:243)=0;
    %data(1:48)= NaN;
    %data(end:-1:220)=NaN;
    algus(end+1)=NaN;         % algus = NaN
    lopp(end+1)=NaN;           % lopp = NaN
    andmed = cat(2,andmed,data);
end



k = plot(data)
hold on
k.LineWidth = 0.8
title('Kogu jää pindala 1993/94-2020/21')
xlabel('Päevade arv alates 1-st oktoobrist') 
ylabel('Jää pindala [km^2]') 
disp(aasta);
end

%%
pikkused = pikkused(:);
aastad = aastad(:);
algus=algus(:);
lopp=lopp(:);

tabel = table(aastad,algus,lopp,pikkused);
%% Trendid

lm1 = fitlm(aastad,algus)
figure(1)
subplot(2,2,1)
plot(lm1)
title('aastad/algus')

lm2 = fitlm(aastad,lopp)
subplot(2,2,2)
plot(lm2)
title('aastad/lopp')

lm3 = fitlm(aastad,pikkused)
subplot(2,2,3)
plot(lm3)
title('aastad/pikkused')

% Kogu jaapindala aastas
sum_pindala = nansum(andmed)
lm4 = fitlm(aastad,sum_pindala)
subplot(2,2,4)
plot(lm4)
title('aastad/pindala')


%% Statistika tabel
stats = ["min";"max";"mean";"std";"trend";"p-value"];

algus_stat = [min(algus),max(algus),nanmean(algus),nanstd(algus),table2array(lm1.Coefficients(2,1)),table2array(lm1.Coefficients(2,4))];
algus_stat=algus_stat(:);
lopp_stat = [min(lopp),max(lopp),nanmean(lopp),nanstd(lopp),table2array(lm2.Coefficients(2,1)),table2array(lm2.Coefficients(2,4))];
lopp_stat=lopp_stat(:);
pikkused_stat = [min(pikkused),max(pikkused),nanmean(pikkused),nanstd(pikkused),table2array(lm3.Coefficients(2,1)),table2array(lm3.Coefficients(2,4))];
pikkused_stat=pikkused_stat(:);

% mean paevade kaupa ja std keskmistest paevade vaartustest
ice_ex = [min(min(andmed)),max(max(andmed)),nanmean(nanmean(andmed,2)),nanstd(nanmean(andmed,2)),table2array(lm4.Coefficients(2,1)),table2array(lm4.Coefficients(2,4))]
ice_ex=ice_ex(:)

stat = table(stats,algus_stat,lopp_stat,pikkused_stat,ice_ex)


%% main plot
data = andmed;
data(1:48,:) = NaN;  % teeme uue andmestiku kust kustutame 0, et joonis oleks ilusam
data(end:-1:220,:) = NaN;

f = figure(1);
f.Position = [100 100 1200 500];

max_sum = max(transpose(data));
plot(max_sum,'red','LineWidth',1,'LineStyle','--')
hold on

%keskmine jaa pindala paevas ule koikide aastate ja std-d
mean_fsum = nanmean(data,2);
%std = nanstd(mean_fsum);
std = nanstd(data, 1, 2);
std_plus = mean_fsum + std;
std_minus = mean_fsum - std;

p = 0:250;
p = p(:);

plot(mean_fsum,'k','LineWidth',2)
set(gca,'xtick',p(1:15:end));
set(gca,'ytick',0:3000:15000);
hold on

x = 1:243;
y1 = transpose(std_plus);
y2 = transpose(std_minus);
fill([x(49:219) fliplr(x(49:219))], [y1(49:219) fliplr(y2(49:219))],...
    'b','FaceAlpha',0.3,'EdgeColor','k','LineWidth',1)
hold on

 %min algus/lopp ja mean algus/lopp +- std plot
plot(min(algus),0,'.','Color','k','MarkerSize', 30)
hold on
plot(max(lopp),0,'.','Color','k','MarkerSize', 30)

hold on
errorbar(nanmean(algus),0,nanstd(algus),'horizontal','Color','k','Marker','v','MarkerFaceColor','red','MarkerSize',12,'LineWidth',2.5,"CapSize",16)
hold on
errorbar(nanmean(lopp),0,nanstd(lopp),'horizontal','Color','k','Marker','^','MarkerFaceColor','red','MarkerSize',12,'LineWidth',2.5,"CapSize",16)
hold on


grid on;
lgd = legend('$S^{\ast}_{max_i}$','$\bar{S}^{\ast}_i$','$\bar{S}^{\ast}_i \pm \sigma_{S^{\ast}_i}$','$c_{a_{min}}$','$c_{l_{max}}$','$\bar{c}_a \pm \sigma_{c_a}$',...
    '$\bar{c}_l \pm \sigma_{c_l}$','Interpreter','latex')


lgd.FontSize = 20;
lgd.Location = 'northwest';

axis([0 250 0 15000])

xlabel('Päevade arv alates 1-st oktoobrist','fontweight','bold','FontSize', 20) 
ylabel('Jää pindala [km^2]','fontweight','bold','FontSize', 20)

ax = gca;
ax.XAxis.FontSize = 16;
ax.YAxis.FontSize = 16;
ax.LineWidth = 1;


%% koik pindalad 3D

%paevad = 1:243;
% miks ta ei joonesta 2019 ja 2020
%fake = zeros(243,1)

%data = cat(2,data,fake);  % Kui runnid esimest korda votta komeaat maha
%years = 1993:2021;
%years = years(:);

%pcolor(paevad,years,transpose(data));


%shading flat
%cb=colorbar;
%grid on

%cmap=colormap;
%cmap(1,:)=[1,1,1];colormap(cmap)
%set(gca,'ytick',aastad(1:2:end-1))
%saxis([49 220 1993 2021 0 15000])






