path = '/Users/iljagolikov/Documents/ssh/study_Ilja/baka/period/fsum_icevol_'
pikkused = [];
aastad = 1993:1:2020;
algus= [];

lopp=[];

andmed = [];

for aasta = string(1993:1:2020)

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

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

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

if sum(data)~=0
    %data(1:st)= NaN;
    %data(end:-1:ed)=NaN; 
    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;           % andmed  = 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ää ruumala 1993/94-2020/21')
xlabel('Päevade arv alates 1-st oktoobrist') 
ylabel('Jää ruumala [km^3]') 
disp(aasta);
end


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

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

%% Trendid
f = figure(2);
f.Position = [100 100 1200 1000];

lm1 = fitlm(aastad,algus)
subplot(3,1,1)
p1 = plot(lm1)

p1(1).LineWidth = 1.5;
p1(2).LineWidth = 2;
p1(3).LineWidth = 1.5;
p1(4).LineWidth = 1.5;

p1(1).MarkerSize = 8;
p1(3).LineStyle = '--';
p1(4).LineStyle = '--';

text(0.01, 0.9, '(a)', 'Units', 'normalized', 'FontSize', 24)
xlabel('Aastad','fontweight','bold','FontSize', 20) 
ylabel('Päevad','fontweight','bold','FontSize', 20)
set(gca,'xtick',aastad(1:3:end));
set(gca,'ytick',0:25:150);
axis([1990 2022 0 150])
grid on
lgd = legend()
lgd.Location = 'southwest';
lgd.FontSize = 14;

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

lm2 = fitlm(aastad,lopp)
subplot(3,1,2)
p2 = plot(lm2)

p2(1).LineWidth = 1.5;
p2(2).LineWidth = 2;
p2(3).LineWidth = 1.5;
p2(4).LineWidth = 1.5;

p2(1).MarkerSize = 8;
p2(3).LineStyle = '--';
p2(4).LineStyle = '--';

xlabel('Aastad','fontweight','bold','FontSize', 20) 
ylabel('Päevad','fontweight','bold','FontSize', 20)
set(gca,'xtick',aastad(1:3:end));
set(gca,'ytick',100:25:250);
axis([1990 2022 100 250])
grid on
text(0.01, 0.9, '(b)', 'Units', 'normalized', 'FontSize', 24)
lgd = legend()
lgd.Location = 'southwest';
lgd.FontSize = 14;

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

lm3 = fitlm(aastad,pikkused)
subplot(3,1,3)
p3 = plot(lm3)

p3(1).LineWidth = 1.5;
p3(2).LineWidth = 2;
p3(3).LineWidth = 1.5;
p3(4).LineWidth = 1.5;

p3(1).MarkerSize = 8;
p3(3).LineStyle = '--';
p3(4).LineStyle = '--';

xlabel('Aastad','fontweight','bold','FontSize', 20) 
ylabel('Päevad','fontweight','bold','FontSize', 20)
set(gca,'xtick',aastad(1:3:end));
set(gca,'ytick',0:25:200);
axis([1990 2022 0 200])
grid on
text(0.01, 0.9, '(c)', 'Units', 'normalized', 'FontSize', 24)
lgd = legend()
lgd.Location = 'southwest';
lgd.FontSize = 14;

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

%% 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),mean(pikkused),std(pikkused),table2array(lm3.Coefficients(2,1)),table2array(lm3.Coefficients(2,4))];
pikkused_stat=pikkused_stat(:);

% mean paevade kaupa ja std keskmistest paevade vaartustest
ice_v = [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_v=ice_v(:)

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



%% main plot
data = andmed;
data(1:48,:) = NaN;
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_Vsum = nanmean(data,2);
%std = nanstd(mean_fsum);
std = nanstd(data, 1, 2);
std_plus = mean_Vsum + std;
std_minus = mean_Vsum - std;

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

plot(mean_Vsum,'k','LineWidth',2)
set(gca,'xtick',p(1:15:end));
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; hold on
lgd = legend('$V^{\ast}_{max_i}$','$\bar{V}^{\ast}_i$','$\bar{V}^{\ast}_i \pm \sigma_{V^{\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';

%title('Kogu jää pindala 1993/94-2020/21')
xlabel('Päevade arv alates 1-st oktoobrist','fontweight','bold','FontSize', 20) 
ylabel('Jää ruumala [km^3]','fontweight','bold','FontSize', 20)

axis([0 250 0 7])

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


