%% gfs analysis 21 days before

% m_map
%dir='/home/chris/Documentos/Toolbox_WRF_matlab';% path to m_map
%addpath([dir,'/funciones/m_map']);% m_map

dir = pwd;
addpath([pwd, '/../output'])

% We verify that the GFS operational file is complete
file = 'forecast_16_days_after_oper.nc';
time = ncread(file, 'time');% forecast last 65 times (6 hourly timestep)
if length(time)<65
    error('The file has less than 65 records, try the next execution!, this error message was made on purpose...')
end

file = 'gfs_data.nc';% main file GFS analysis
time = ncread(file, 'time');% 6 hourly timestep

% info=ncinfo(file);
% inicio=info.Variables(1).Attributes(3).Value;
% 
% a = datenum(datetime(str2num(inicio(12:15)),str2num(inicio(17:18))...
%     ,str2num(inicio(20:21)),str2num(inicio(23:24)),str2num(inicio(26:27)),str2num(inicio(29:30))));%'19-Dec-2009 23:00:00' Fecha de referencia

a = datenum(2019, 8, 13, 0, 0, 0);% reference date gfs_data.nc file

% clear hrs out out2
% for t = 1:length(time)
%     hrs = time(t);% hours from the reference date
%     out{t} = datevec(a + (hrs/24)); % vector [year month day hour]
%     out2(t) = a + (hrs/24);
% end
out2 =  datenum(a + hours(time));

% spatial coordinates
lon = ncread(file, 'lon');
lon = double(lon);
lat = ncread(file, 'lat');
lat = double(lat);

rm = 40;% days to average

% west center
lat1 = find(lat == -52);
lon1 = find(lon == 14);
startLocW = [lon1-16 lat1-8 1 1];% the resolution is 0.5 degrees, if this is not correct change the positions to extract! 
countW = [2*16+1 2*8+1 1 Inf]; 
strideW = [1 1 1 1];% the file is [lon lat lev time] 

gh = ncread(file, 'gh', startLocW, countW, strideW);% read geopotential height
gh = double(squeeze(gh));
aux = NaN(size(gh));
gh_west = gh;
gh_backup_west = gh;

for t = 1 + size(gh(:,:,1:end-64),3)-61*4:size(gh(:,:,1:end-64),3)% 40 days running mean 61 days before
    for i = 1:size(gh,1)
        for j = 1:size(gh,2)
            aux(i,j,t) = (gh(i,j,t)-nanmean(gh(i,j,t-rm*4:t-1),3))./nanstd(gh(i,j,t-rm*4:t-1),[],3);
        end
    end
end
ghwest = nanmean(nanmean(aux,1),2);
ghwest = squeeze(ghwest);

% east center
lat2 = find(lat == -58);
lon2 = find(lon == 66);
startLocE = [lon2-16 lat2-8 1 1]; 
countE = [2*16+1 2*8+1 1 Inf]; 
strideE = [1 1 1 1]; 

gh = ncread(file, 'gh', startLocE, countE, strideE);
gh = double(squeeze(gh));
aux = NaN(size(gh));
gh_east = gh;
gh_backup_east = gh;

for t=1+size(gh(:,:,1:end-64),3)-61*4:size(gh(:,:,1:end-64),3)% 40 days running mean 61 days before
    for i=1:size(gh,1)
        for j=1:size(gh,2)
            aux(i,j,t)=(gh(i,j,t)-nanmean(gh(i,j,t-rm*4:t-1),3))./nanstd(gh(i,j,t-rm*4:t-1),[],3);
        end
    end
end
gheast = nanmean(nanmean(aux,1),2);
gheast = squeeze(gheast);

% calc of ETI
eti = gheast-ghwest;% east-west

seti = NaN(size(eti));
for t=size(gh(:,:,1:end-64),3)-21*4:size(gh(:,:,1:end-64),3)% 40 dias running mean 21 days before
    seti(t)=(eti(t)-nanmean(eti(t-rm*4:t-1),1))./nanstd(eti(t-rm*4:t-1),[],1);%cambiar a 40
end

anl_interval = size(gh(:,:,1:end-64),3)-21*4:size(gh(:,:,1:end-64),3);% interval 21 days before

h = figure('Position', [100, 100, 960, 620]);
ax(1) = subplot(2,1,1);
plot(ax(1),out2(anl_interval),seti(anl_interval),'linewidth',1.5,'Color',[0 0.3 0.9])
hold on

series.eti_anl = eti(anl_interval);% backup serie!
series.seti_anl = seti(anl_interval);% backup serie!
series.time_anl = out2(anl_interval);% backup serie!
series.time_forecast = out2(end-64:end);% backup serie!

% series.eti_anl = eti;% backup serie!
% series.seti_anl = seti;% backup serie!
% series.time_anl = out2;% backup serie!
% series.time_forecast = out2(end-64:end);% backup serie!

clear i j aux a file lat lon lat1 lat2 lon1 lon2 t time hrs anl_interval gheast ghwest eti seti out

%% loop ensemble gfs forecast 16 days after

for em=1:20    
    file = ['forecast_16_days_after_e',num2str(em,'%02d'),'.nc'];
    
    eval(['series.eti_e',num2str(em,'%02d'),' = [];']);
    eval(['series.seti_e',num2str(em,'%02d'),' = [];']);
    if isfile(['../output/',file])
        disp(['Ensemble member ',num2str(em,'%02d'),' exists.'])
        
        % west center
        gh_e = ncread(file, 'gh', startLocW, countW, strideW);
        gh_e = double(squeeze(gh_e));
        
        if size(gh_e,3) == 65
            aux = NaN(size(gh_west));
            gh_west(:,:,end-64:end) = gh_e(:,:,end-64:end);
            
            for t=1+size(gh_west(:,:,1:end-64),3)-40*4:size(gh_west,3)
                for i=1:size(gh_west,1)
                    for j=1:size(gh_west,2)
                        aux(i,j,t) = (gh_west(i,j,t)-nanmean(gh_west(i,j,t-rm*4:t-1),3))./nanstd(gh_west(i,j,t-rm*4:t-1),[],3);
                    end
                end
            end
            ghwest = nanmean(nanmean(aux,1),2);
            ghwest = squeeze(ghwest);
            
            % east center
            gh_e = ncread(file,'gh',startLocE,countE,strideE);
            gh_e = double(squeeze(gh_e));
            
            aux = NaN(size(gh_east));
            gh_east(:,:,end-64:end)=gh_e(:,:,end-64:end);
            
            for t=1+size(gh_east(:,:,1:end-64),3)-40*4:size(gh_east,3)
                for i=1:size(gh_east,1)
                    for j=1:size(gh_east,2)
                        aux(i,j,t) = (gh_east(i,j,t)-nanmean(gh_east(i,j,t-rm*4:t-1),3))./nanstd(gh_east(i,j,t-rm*4:t-1),[],3);
                    end
                end
            end
            gheast = nanmean(nanmean(aux,1),2);
            gheast = squeeze(gheast);
            
            eti = gheast-ghwest;
            
            seti = NaN(size(eti));
            for t=size(gh(:,:,1:end-64),3):size(gh,3)
                seti(t) = (eti(t)-nanmean(eti(t-rm*4:t-1),1))./nanstd(eti(t-rm*4:t-1),[],1);%cambiar a 40
            end
            
            plot(ax(1),out2(end-64:end),seti(end-64:end),'linewidth',1.5,'Color',[0.9 0.6 0.5])
            
            eval(['series.eti_e',num2str(em,'%02d'),' = eti(end-64:end);']);
            eval(['series.seti_e',num2str(em,'%02d'),' = seti(end-64:end);']);            
        else
            disp(['Ensemble member ',num2str(em,'%02d'),' has less than 65 records... is not included.'])            
        end
    else
        disp(['Ensemble member ',num2str(em,'%02d'),' does not exist.'])        
    end
end

clear i j t seti eti file gheast ghwest em

%% gfs forecast 16 days after (ensemble average)
   
file=['forecast_16_days_after_eavg.nc'];

if isfile(['../output/',file])
    disp(['Ensemble average exists.'])
    
    % west center
    gh_e = ncread(file,'gh',startLocW,countW,strideW);
    gh_e = double(squeeze(gh_e));
    if size(gh_e,3) == 65
        
        aux = NaN(size(gh_west));
        gh_west(:,:,end-64:end) = gh_e(:,:,end-64:end);
        
        for t=1+size(gh_west(:,:,1:end-64),3)-40*4:size(gh_west,3)
            for i=1:size(gh_west,1)
                for j=1:size(gh_west,2)
                    aux(i,j,t) = (gh_west(i,j,t)-nanmean(gh_west(i,j,t-rm*4:t-1),3))./nanstd(gh_west(i,j,t-rm*4:t-1),[],3);
                end
            end
        end
        ghwest = nanmean(nanmean(aux,1),2);
        ghwest = squeeze(ghwest);
        
        % east center
        gh_e = ncread(file,'gh',startLocE,countE,strideE);
        gh_e = double(squeeze(gh_e));
        
        aux = NaN(size(gh_east));
        gh_east(:,:,end-64:end) = gh_e(:,:,end-64:end);
        
        for t=1+size(gh_east(:,:,1:end-64),3)-40*4:size(gh_east,3)
            for i=1:size(gh_east,1)
                for j=1:size(gh_east,2)
                    aux(i,j,t) = (gh_east(i,j,t)-nanmean(gh_east(i,j,t-rm*4:t-1),3))./nanstd(gh_east(i,j,t-rm*4:t-1),[],3);
                end
            end
        end
        gheast = nanmean(nanmean(aux,1),2);
        gheast = squeeze(gheast);
        
        eti = gheast-ghwest;
        
        seti=NaN(size(eti));
        for t=size(gh(:,:,1:end-64),3):size(gh,3)
            seti(t) = (eti(t)-nanmean(eti(t-rm*4:t-1),1))./nanstd(eti(t-rm*4:t-1),[],1);%cambiar a 40
        end
        pl1 = plot(ax(1),out2(end-64:end),seti(end-64:end),'linewidth',1.5,'Color',[1 0 0]);
        
        series.eti_eavg = eti(end-64:end);
        series.seti_eavg = seti(end-64:end);
    else
        disp(['Ensemble average has less than 65 records... is not included.'])
        series.eti_eavg = [];
        series.seti_eavg = [];
        pl1=[];
    end
else
    disp(['Ensemble average does not exist.'])
    series.eti_eavg = [];
    series.seti_eavg = [];
    pl1=[];
end

clear i j t seti eti file gheast ghwest

%% gfs forecast 16 days after (gfs operational)

file = 'forecast_16_days_after_oper.nc';

if isfile(['../output/',file])
    disp(['GFS operationnal exists.'])
    
    % west center
    gh_e = ncread(file,'gh',startLocW,countW,strideW);
    gh_e = double(squeeze(gh_e));
    if size(gh_e,3)==65
        
        aux = NaN(size(gh_west));
        gh_west(:,:,end-64:end)=gh_e(:,:,end-64:end);
        
        for t=1+size(gh_west(:,:,1:end-64),3)-40*4:size(gh_west,3)
            for i=1:size(gh_west,1)
                for j=1:size(gh_west,2)
                    aux(i,j,t) = (gh_west(i,j,t)-nanmean(gh_west(i,j,t-rm*4:t-1),3))./nanstd(gh_west(i,j,t-rm*4:t-1),[],3);
                end
            end
        end
        ghwest = nanmean(nanmean(aux,1),2);
        ghwest = squeeze(ghwest);
        
        % east center
        gh_e = ncread(file,'gh',startLocE,countE,strideE);
        gh_e = double(squeeze(gh_e));
        
        aux = NaN(size(gh_east));
        gh_east(:,:,end-64:end) = gh_e(:,:,end-64:end);
        
        for t=1+size(gh_east(:,:,1:end-64),3)-40*4:size(gh_east,3)
            for i=1:size(gh_east,1)
                for j=1:size(gh_east,2)
                    aux(i,j,t) = (gh_east(i,j,t)-nanmean(gh_east(i,j,t-rm*4:t-1),3))./nanstd(gh_east(i,j,t-rm*4:t-1),[],3);
                end
            end
        end
        gheast = nanmean(nanmean(aux,1),2);
        gheast = squeeze(gheast);
        
        eti = gheast-ghwest;
        
        seti = NaN(size(eti));
        for t=size(gh(:,:,1:end-64),3):size(gh,3)
            seti(t) = (eti(t)-nanmean(eti(t-rm*4:t-1),1))./nanstd(eti(t-rm*4:t-1),[],1);%cambiar a 40
        end
        pl2 = plot(ax(1),out2(end-64:end),seti(end-64:end),'linewidth',1.5,'Color',[0 0.5 0.2]);
        
        series.eti_oper = eti(end-64:end);
        series.seti_oper = seti(end-64:end);
    else
        disp(['GFS operationnal has less than 65 records... is not included.'])
        pl2=[];
    end
else
    disp(['GFS operationnal does not exist.'])
    pl2=[];
end

clear i j t seti eti file gheast ghwest

%% figure settings (top panel)

% title
%t = title(['From GFS analysis and forecast (21+16 days) ',datestr(out2(end-64),'dd-mmm-yyyy HH:MM:SS')],'FontWeight','b','FontSize',13);
%t1 = get(t, 'position');
%set(t, 'position', [t1(1) t1(2)+2.3 t1(3)])

% ylabel
%ylabel({'Standardized';'Extra Tropical Index [SETI]'},'FontWeight','b','FontSize',12)
ylabel({'SETI [\sigma]'},'FontWeight','b','FontSize',12)

% settings ylabel
ylim([-5 5])

NumTicks = 11;
L = get(gca,'YLim');
set(gca,'YTick',linspace(L(1),L(2),NumTicks))

% xlabel
xlabel([],'FontWeight','b','FontSize',12)

% settings xlabel
xlim([out2(size(gh(:,:,1:end-64),3)-21*4) out2(end)])
set(ax(1),'XTick',out2(size(gh(:,:,1:end-64),3)-21*4:8:end))

datetick(ax(1),'x','dd','keeplimits','keepticks')

% general settigns
set(gca,'TickDir','out');
grid on

plot([1 1]*out2(end-64),ylim,'--','Color',[0 0 0.3],'linewidth',1.5)

% horizontal patch settings
limx = xlim;
f1 = patch([limx(1) limx(2) limx(2) limx(1)],[-1 -1 1 1],[0.5608    0.6980    0.5373],'LineStyle','none','FaceAlpha',0.4);
uistack(f1,'bottom')

% legend
if (~isempty(pl1) & ~isempty(pl2))
    legend([pl1 pl2],'GEFS avg','GFS','location','southeast')
elseif (~isempty(pl1) & isempty(pl2))
    legend([pl1],'GEFS avg','location','southeast')
elseif (isempty(pl1) & ~isempty(pl2))
    legend([pl2],'GFS','location','southeast')
else
end

dim = [.14 .61 .3 .3];

texto = {['From GFS analysis and forecast (21+16 days) '];...
    ['                    ',datestr(out2(end-64),'dd-mmm-yyyy HH:MM:SS')]};
ann = annotation('textbox',dim,'String',texto,...
    'backgroundcolor','w','FitBoxToText','on');
ann.FontWeight = 'bold';
ann.FontSize = 10;

%% figure settings (bottom panel)

%data = readtable([pwd,'/../backup/mjo_index'], 'HeaderLines', 2);
%data = readtable([pwd,'/../data/series/mjo_index'], 'HeaderLines', 2);
%if (~isempty(data))
%    data = table2array(data(:,1:7));
%    auxtime = datenum(data(:,1),data(:,2),data(:,3),0,0,0);
%    amp = data(:,7);
%    pha = data(:,6);
%end
%clear data

fname = [pwd,'/../data/series/amp.nc'];
time1 = ncread(fname, 'T');
time1 = datenum(datetime(time1, 'convertfrom', 'juliandate'));
amp1 = ncread(fname, 'amplitude');

fname = [pwd,'/../data/series/pha.nc'];
time2 = ncread(fname, 'T');
time2 = datenum(datetime(time2, 'convertfrom', 'juliandate'));
pha2 = ncread(fname, 'phase');


auxtime = [datenum(date)-100: datenum(date)]';
auxtime = datenum(auxtime + duration(12, 0, 0));

amp = NaN(size(auxtime));
pha = NaN(size(auxtime));

for i = 1:length(auxtime)
    idx1 = auxtime(i) == time1;
    idx2 = auxtime(i) == time2;
    
    if ~isempty(time1(idx1))
        amp(i) = amp1(idx1);
    end
    if ~isempty(time2(idx2))
        pha(i) = pha2(idx2);
    end
end

auxtime = auxtime(~isnan(pha));
amp = amp(~isnan(pha));
pha = pha(~isnan(pha));

auxtime = auxtime(~isnan(amp));
pha = pha(~isnan(amp));
amp = amp(~isnan(amp));

nphase = 8;

auxcmap = [0.1961    0.5333    0.7412;
           0.4000    0.7608    0.6471;
           0.6706    0.8667    0.6431;
           0.9020    0.9608    0.5961;
           0.9961    0.8784    0.5451;
           0.9569    0.4275    0.2627;
           0.7098    0.2314    0.1294;
           0.8118    0.5020    0.8745];

ax(2) = subplot(2,1,2);
yyaxis right
ax(2).YAxis(2).Color = 'k';
f1 = patch([limx(1) limx(2) limx(2) limx(1)],[0 0 1 1],[0.5608    0.6980    0.5373],'LineStyle','none','FaceAlpha',0.4);
hold on
if exist('auxtime','var')
    if (~isempty(auxtime))
        for t=1:size(auxtime,1)
            if pha(t)>8
            else
                scatter(auxtime(t),amp(t),'ok','SizeData',80,'MarkerFaceColor',[auxcmap(pha(t),:)],'MarkerEdgeColor','k')
                hold on
            end
        end
    end
end
ylabel(['MJO Amplitude'],'FontWeight','b','FontSize',12)
ylim([0 6])
NumTicks = 7;
L = get(gca,'YLim');
set(gca,'YTick',linspace(L(1),L(2),NumTicks))

yyaxis left
ax(2).YAxis(1).Color = 'k';

clear p

% codigo estacion general freire
cod = 340031;

% read Tx
data = readtable([pwd,'/../data/series/T_max'], 'HeaderLines', 0);
data = table2array(data);

cods = NaN(length(data(1,:)),1);
for i = 2:length(data(1,:))
    cods(i) = str2double(data(1,i));    
end

idx = find(cods == cod);

if (~isempty(idx))
    clear time
    clear value
    cont = 1;
    for i = 16:length(data(:,1))
        time(cont,1) = datenum(data{i,1});
        value(cont,1) = str2double(data(i,idx));
        cont = cont + 1;
    end    
    value(value==-9999)=NaN;
    p(1) = plot(time,value,'-','linewidth',1.5,'Color',[1 0.3 0.3],'DisplayName','T max');
    hold on
end

% read Tm
data = readtable([pwd,'/../data/series/T_mean'], 'HeaderLines', 0);
data = table2array(data);

cods = NaN(length(data(1,:)),1);
for i = 2:length(data(1,:))
    cods(i) = str2double(data(1,i));    
end

idx = find(cods == cod);

if (~isempty(idx))
    clear time
    clear value
    cont = 1;
    for i = 16:length(data(:,1))
        time(cont,1) = datenum((data{i,1}));
        value(cont,1) = str2double(data(i,idx));
        cont = cont + 1;
    end    
    value(value==-9999)=NaN;
    p(2) = plot(time,value,'-','linewidth',1.5,'Color',[0.2 0.5 0.5],'DisplayName','T mean');
    hold on
end

% read Tn
data = readtable([pwd,'/../data/series/T_min'], 'HeaderLines', 0);
data = table2array(data);

cods = NaN(length(data(1,:)),1);
for i = 2:length(data(1,:))
    cods(i) = str2double(data(1,i));    
end

idx = find(cods == cod);

if (~isempty(idx))
    clear time
    clear value
    cont = 1;
    for i = 16:length(data(:,1))
        time(cont,1) = datenum(data{i,1});
        value(cont,1) = str2double(data(i,idx));
        cont = cont + 1;
    end    
    value(value==-9999)=NaN;
    p(3) = plot(time,value,'-','linewidth',1.5,'Color',[0.1 0 0.8],'DisplayName','T min');
    hold on
end


altura = data{5,idx};
latitud = abs(str2double(data{6,idx}));
longitud = abs(str2double(data{7,idx}));
name_estacion = data{4,idx};

% t=title(['Curico DGA 34.9811Â°S, 71.2361Â°W, 195m'],'FontWeight','b','FontSize',10.3);
t = title([name_estacion,' ',num2str(latitud),'°S, ',num2str(longitud),'°W, ',altura,'m'],'FontWeight','b','FontSize',10.3);
set(t, 'horizontalAlignment', 'right')
set(t, 'units', 'normalized')
t1 = get(t, 'position');
set(t, 'position', [1 t1(2) t1(3)])
ylabel(['Temperature [°C]'],'FontWeight','b','FontSize',12)
ylim([-10 40])
xlim([out2(size(gh(:,:,1:end-64),3)-21*4) out2(end)])
set(ax(2),'XTick',out2(size(gh(:,:,1:end-64),3)-21*4:8:end))
datetick('x','dd','keeplimits','keepticks')
set(gca,'TickDir','out');
grid on
plot([1 1]*out2(end-64),ylim,'--','Color',[0 0 0.3],'linewidth',1.5)

if (~isempty(p))
    legend([p(1) p(2) p(3)],'T max','T mean','T min')
end

barra = colorbar('Southoutside','FontWeight','bold','Linewidth',1,'XTickLabel',0:8);
ylabel(barra,'MJO Phase','FontWeight','bold','FontSize',12);

colormap(auxcmap)
caxis([0 8])

saveas(h,[pwd,'/../figures/SETI forecast From GFS ',datestr(out2(end-64),'yyyy-mm-dd HH:MM:SS'),'.png'])
close(h)

save('../output/series.mat','series');
save('../upload/series.mat','series');

%%

time = [series.time_anl(1:end-1, 1); series.time_forecast];
Time = datetime(datevec(time));
Time.Format =  'uuuu-MM-dd HH:mm:ss';
aux = NaN(size(time));

anl = aux;
anl(1:length(series.seti_anl), 1) = series.seti_anl;

oper = aux;
if ~isempty(series.seti_oper)
oper(length(series.seti_anl):end, 1) = series.seti_oper;
end

eavg = aux;
if ~isempty(series.seti_eavg)
eavg(length(series.seti_anl):end, 1) = series.seti_eavg;
end

e1 = aux;
if ~isempty(series.seti_e01)
e1(length(series.seti_anl):end, 1) = series.seti_e01;
end

e2 = aux;
if ~isempty(series.seti_e02)
e2(length(series.seti_anl):end, 1) = series.seti_e02;
end

e3 = aux;
if ~isempty(series.seti_e03)
e3(length(series.seti_anl):end, 1) = series.seti_e03;
end

e4 = aux;
if ~isempty(series.seti_e04)
e4(length(series.seti_anl):end, 1) = series.seti_e04;
end

e5 = aux;
if ~isempty(series.seti_e05)
e5(length(series.seti_anl):end, 1) = series.seti_e05;
end

e6 = aux;
if ~isempty(series.seti_e06)
e6(length(series.seti_anl):end, 1) = series.seti_e06;
end

e7 = aux;
if ~isempty(series.seti_e07)
e7(length(series.seti_anl):end, 1) = series.seti_e07;
end

e8 = aux;
if ~isempty(series.seti_e08)
e8(length(series.seti_anl):end, 1) = series.seti_e08;
end

e9 = aux;
if ~isempty(series.seti_e09)
e9(length(series.seti_anl):end, 1) = series.seti_e09;
end

e10 = aux;
if ~isempty(series.seti_e10)
e10(length(series.seti_anl):end, 1) = series.seti_e10;
end

e11 = aux;
if ~isempty(series.seti_e11)
e11(length(series.seti_anl):end, 1) = series.seti_e11;
end

e12 = aux;
if ~isempty(series.seti_e12)
e12(length(series.seti_anl):end, 1) = series.seti_e12;
end

e13 = aux;
if ~isempty(series.seti_e13)
e13(length(series.seti_anl):end, 1) = series.seti_e13;
end

e14 = aux;
if ~isempty(series.seti_e14)
e14(length(series.seti_anl):end, 1) = series.seti_e14;
end

e15 = aux;
if ~isempty(series.seti_e15)
e15(length(series.seti_anl):end, 1) = series.seti_e15;
end

e16 = aux;
if ~isempty(series.seti_e16)
e16(length(series.seti_anl):end, 1) = series.seti_e16;
end

e17 = aux;
if ~isempty(series.seti_e17)
e17(length(series.seti_anl):end, 1) = series.seti_e17;
end

e18 = aux;
if ~isempty(series.seti_e18)
e18(length(series.seti_anl):end, 1) = series.seti_e18;
end

e19 = aux;
if ~isempty(series.seti_e19)
e19(length(series.seti_anl):end, 1) = series.seti_e19;
end

e20 = aux;
if ~isempty(series.seti_e20)
e20(length(series.seti_anl):end, 1) = series.seti_e20;
end

TT = timetable(Time, anl, oper, eavg, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,...
    e11, e12, e13, e14, e15, e16, e17, e18, e19, e20);

TT2 = TT;
TT3 = timetable2table(TT2); 
writetable(TT3,'../output/SETI_all.csv');
writetable(TT3,'../upload/SETI_all.csv');

%%

% seti_all = NaN(length([series.seti_anl; NaN(size(series.time_forecast))]), 24);
% seti_all(:,1) = [series.time_anl; series.time_forecast];
% seti_all(:,2) = [series.seti_anl; NaN(size(series.time_forecast))];
% 
% if size(series.seti_oper,1)==65
% seti_all(:,3) = [NaN(size(series.time_anl)); series.seti_oper];
% end
% if size(series.seti_eavg,1)==65
% seti_all(:,4) = [NaN(size(series.time_anl)); series.seti_eavg];
% end
% 
% if size(series.seti_e01,1)==65
% seti_all(:,5) = [NaN(size(series.time_anl)); series.seti_e01];
% end
% if size(series.seti_e02,1)==65
% seti_all(:,6) = [NaN(size(series.time_anl)); series.seti_e02];
% end
% if size(series.seti_e03,1)==65
% seti_all(:,7) = [NaN(size(series.time_anl)); series.seti_e03];
% end
% if size(series.seti_e04,1)==65
% seti_all(:,8) = [NaN(size(series.time_anl)); series.seti_e04];
% end
% if size(series.seti_e05,1)==65
% seti_all(:,9) = [NaN(size(series.time_anl)); series.seti_e05];
% end
% if size(series.seti_e06,1)==65
% seti_all(:,10) = [NaN(size(series.time_anl)); series.seti_e06];
% end
% if size(series.seti_e07,1)==65
% seti_all(:,11) = [NaN(size(series.time_anl)); series.seti_e07];
% end
% if size(series.seti_e08,1)==65
% seti_all(:,12) = [NaN(size(series.time_anl)); series.seti_e08];
% end
% if size(series.seti_e09,1)==65
% seti_all(:,13) = [NaN(size(series.time_anl)); series.seti_e09];
% end
% if size(series.seti_e10,1)==65
% seti_all(:,14) = [NaN(size(series.time_anl)); series.seti_e10];
% end
% if size(series.seti_e11,1)==65
% seti_all(:,15) = [NaN(size(series.time_anl)); series.seti_e11];
% end
% if size(series.seti_e12,1)==65
% seti_all(:,16) = [NaN(size(series.time_anl)); series.seti_e12];
% end
% if size(series.seti_e13,1)==65
% seti_all(:,17) = [NaN(size(series.time_anl)); series.seti_e13];
% end
% if size(series.seti_e14,1)==65
% seti_all(:,18) = [NaN(size(series.time_anl)); series.seti_e14];
% end
% if size(series.seti_e15,1)==65
% seti_all(:,19) = [NaN(size(series.time_anl)); series.seti_e15];
% end
% if size(series.seti_e16,1)==65
% seti_all(:,20) = [NaN(size(series.time_anl)); series.seti_e16];
% end
% if size(series.seti_e17,1)==65
% seti_all(:,21) = [NaN(size(series.time_anl)); series.seti_e17];
% end
% if size(series.seti_e18,1)==65
% seti_all(:,22) = [NaN(size(series.time_anl)); series.seti_e18];
% end
% if size(series.seti_e19,1)==65
% seti_all(:,23) = [NaN(size(series.time_anl)); series.seti_e19];
% end
% if size(series.seti_e20,1)==65
% seti_all(:,24) = [NaN(size(series.time_anl)); series.seti_e20];
% end
% 
% % aux = length(series.time_anl);
% % for i = 3:size(seti_all,2)
% % seti_all (aux-3:aux,i) = series.time_anl(end-3:end);
% % end
% 
% D = datevec(seti_all(:,1));
% Time = datetime(D(:,1),D(:,2),D(:,3),D(:,4),D(:,5),D(:,6));
% SETI = [seti_all(:,2:end)];
% tt2 = timetable(Time, SETI(:,4), SETI(:,5), SETI(:,6), SETI(:,7),...
%      SETI(:,8), SETI(:,9), SETI(:,10), SETI(:,11), SETI(:,12),...
%       SETI(:,13), SETI(:,14), SETI(:,15), SETI(:,16), SETI(:,17), SETI(:,18),...
%       SETI(:,19), SETI(:,20), SETI(:,21), SETI(:,22), SETI(:,23),...
%       SETI(:,1), SETI(:,2), SETI(:,3));
% % tt = synchronize(tt1,tt2);
% TT2 = retime(tt2,'daily');
% TT3 = timetable2table(TT2); 
% TT3 = table(datestr(datenum(TT2.Time),'yyyy-mm-dd'), TT2.Var1, TT2.Var2, TT2.Var3...
%     , TT2.Var4, TT2.Var5, TT2.Var6, TT2.Var7, TT2.Var8, TT2.Var9, TT2.Var10, TT2.Var11...
%     , TT2.Var12, TT2.Var13, TT2.Var14, TT2.Var15, TT2.Var16, TT2.Var17, TT2.Var18, TT2.Var19...
%     , TT2.Var20, TT2.Var21, TT2.Var22, TT2.Var23);
% writetable(TT3,'../output/SETI_all.csv');
% 
% % D = datevec([series.time_anl; series.time_forecast]);
% % Time = datetime(D(:,1),D(:,2),D(:,3),D(:,4),D(:,5),D(:,6));
% % SETI = [series.seti_anl; NaN(size(series.time_forecast))];
% % tt2 = timetable(Time, SETI);
% % % tt = synchronize(tt1,tt2);
% % TT2 = retime(tt2,'daily');
% % TT3 = timetable2table(TT2); 
% % TT3 = table(datestr(datenum(TT2.Time),'yyyy-mm-dd'), TT2.SETI);
% % writetable(TT3,'../output/SETI_anl.csv');
% 
% % D = datevec([series.time_anl; series.time_forecast]);
% % Time = datetime(D(:,1),D(:,2),D(:,3),D(:,4),D(:,5),D(:,6));
% % SETI = [NaN(size(series.time_anl)); series.seti_oper];
% % tt2 = timetable(Time, SETI);
% % % tt = synchronize(tt1,tt2);
% % TT2 = retime(tt2,'daily');
% % TT3 = timetable2table(TT2); 
% % TT3 = table(datestr(datenum(TT2.Time),'yyyy-mm-dd'), TT2.SETI);
% % writetable(TT3,'../output/SETI_oper.csv');

return
