%clear

IDEAL=1;
Ndim=3;
Amps=[0,1,2,3,4,5,6,7,10];
% Exps={'WA1Jan35_EQ1M40_','WA1Jan35_EQ1M40lon60180_','WA1Jan35_EQ1M40_','WA1Jan35_EQ1M40lon60180_'};
Exps={'WA1Jan35_EQ1M40_','WA1Jan35_EQ1M40lon60180_','WA1Jan35_EQ1M40lon180300_','WA1Jan35_EQ1M40lon30060_','WA1Jan35_EQ1M40_','WA1Jan35_EQ1M40lon60180_'};
appendsss={'_2','_2','','','',''};
linetype={'k','k--+','k--^','k--d','r*','r+'};
legends={'0-360E','60-180E','180-300E','300-60E','0-360E','60-180E'};

% Exps={'WA1Jan35_EQ1M40_','WA1Jan35_EQ1M40lon60180_','WA1Jan35_EQ1M40lon180300_','WA1Jan35_EQ1M40lon30060_'};
% appendsss={'_2','_2','',''};
% linetype={'k','k--+','k--^','k--d'};
% legends={'0-360E','60-180E','180-300E','300-60E'};

% Amps=[0,1,2,3,4,5,6,7,10,30];
% Exps={'WA1Jan35_EQ1M40lon60180_','WA1Jan35_EQ1M40lon180300_','WA1Jan35_EQ1M40lon30060_'};
% appendsss={'_2','',''};
% linetype={'k--+','k--^','k--d'};
% legends={'60-180E','180-300E','300-60E'};

NExps=length(Exps);
showerror=false;
showUk1=false;

addpath('/Users/wanyingkang/Desktop/research/matlab/code/basic/')
addpath('/Users/wanyingkang/Desktop/matlab/code/basic/')
SSW_cis1=zeros(length(Amps),2,NExps);
Uk1s1=zeros(length(Amps),NExps);
figure
append={''};
clf;
hold on;
  SSW_cis1(:,:,1)=(SSW_cis1(:,:,1)+SSW_cis1(:,:,5))./2;
  SSW_cis1(:,:,2)=(SSW_cis1(:,:,2)+SSW_cis1(:,:,6))./2;
  SSW_cis1(1,:,2)=SSW_cis1(1,:,3);
  SSW_cis1(1,:,1)=SSW_cis1(1,:,3);
for iExp=1:6
    [Uk1s1(:,iExp),SSW_cis1(:,:,iExp)]=SSW_confidence(Exps(iExp),Ndim,IDEAL,Amps,appendsss(iExp));
    %if iExp==1
        if showerror 
            errorbar(Amps,mean(squeeze(SSW_cis1(:,:,iExp)),2),(SSW_cis1(:,2,iExp)-SSW_cis1(:,1,iExp))./2,char(linetype(iExp)),'LineWidth',3);
        else
            plot(Amps,mean(squeeze(SSW_cis1(:,:,iExp)),2),char(linetype(iExp)),'LineWidth',3);
        end
%     else
%         if showerror
%             errorbar(Amps,mean(squeeze(SSW_cis1(:,:,iExp)),2),(SSW_cis1(:,2,iExp)-SSW_cis1(:,1,iExp))./2,'k--','LineWidth',3);
%         else
%             plot(Amps,mean(squeeze(SSW_cis1(:,:,iExp)),2),'k--','LineWidth',3);
%         end    
%     end
end
    ax1=gca;
    ax1.YColor='k';
    %set(gcf,'Position',[114   249   564   320]);
    %ax1.Position=pos;
    ylim([0,20])
    xlabel('MJO Amplitude (K/day)')
    ylabel('#SSW per decade')
    set(gca,'FontSize',16)
if showUk1
    ax2=axes('Position',ax1.Position,'YAxisLocation','right','Color','none','YColor','r');
    hold on
for iExp=1:NExps
    if iExp==1
        plot(Amps,Uk1s1(:,iExp),'r','LineWidth',2,'Parent',ax2);
        hold on
    else
        plot(Amps,Uk1s1(:,iExp),'r--','LineWidth',2,'Parent',ax2);
    end
end
     ylabel('zonal asym U k=1 amplitude')
end

    if Ndim==3
        title('#SSW (U pass 0)')
    else
        title('#SSW (U_t exceed 5m/s for 5 days)')
    end
    legend(legends(1:NExps),'FontSize',16)
    set(gca,'FontSize',16)
    if length(Exps)>3
    nSSW=squeeze(mean(SSW_cis1(:,:,2:4),2));
    plot(Amps,sum(nSSW-nSSW(1,:),2)+nSSW(1,3),'r','LineWidth',2)
    end


% SSW_cis2=zeros(length(Amps),2,NExps);
% Uk1s2=zeros(length(Amps),NExps);
% figure
% append={'_2'};
% clf;
% hold on;
% for iExp=1:NExps
%     [Uk1s2(:,iExp),SSW_cis2(:,:,iExp)]=SSW_confidence(Exps(iExp),Ndim,IDEAL,Amps,append);
%     if iExp==1
%         if showerror 
%             errorbar(Amps,mean(squeeze(SSW_cis2(:,:,iExp)),2),(SSW_cis2(:,2,iExp)-SSW_cis2(:,1,iExp))./2,'k','LineWidth',3);
%         else
%             plot(Amps,mean(squeeze(SSW_cis2(:,:,iExp)),2),'k','LineWidth',3);
%         end
%     else
%         if showerror
%             errorbar(Amps,mean(squeeze(SSW_cis2(:,:,iExp)),2),(SSW_cis2(:,2,iExp)-SSW_cis2(:,1,iExp))./2,'k--','LineWidth',3);
%         else
%             plot(Amps,mean(squeeze(SSW_cis2(:,:,iExp)),2),'k--','LineWidth',3);
%         end    
%     end
% end
%     ax1=gca;
%     ax1.YColor='k';
%     %set(gcf,'Position',[114   249   564   320]);
%     %ax1.Position=pos;
%     ylim([0,20])
%     xlabel('MJO Amplitude (K/day)')
%     ylabel('#SSW per decade')
%     set(gca,'FontSize',16)
% if showUk1
%     ax2=axes('Position',ax1.Position,'YAxisLocation','right','Color','none','YColor','r');
%     hold on
% for iExp=1:NExps
%     if iExp==1
%         plot(Amps,Uk1s2(:,iExp),'r','LineWidth',2,'Parent',ax2);
%         hold on
%     else
%         plot(Amps,Uk1s2(:,iExp),'r--','LineWidth',2,'Parent',ax2);
%     end
% end
%      ylabel('zonal asym U k=1 amplitude')
% end
%      
% if Ndim==3
%     title('Group 2 #SSW (U pass 0)')
% else
%     title('Group 2 #SSW (U_t exceed 5m/s for 5 days)')
% end
% set(gca,'FontSize',16)
% 
% 
% figure
% SSW_cismean=zeros(length(Amps),2,NExps);
% Uk1smean=zeros(length(Amps),NExps);
% append={'','_2'};
% hold on;
% for iExp=1:NExps
%     [Uk1smean(:,iExp),SSW_cismean(:,:,iExp)]=SSW_confidence(Exps(iExp),Ndim,IDEAL,Amps,append);
%     if iExp==1
%         if showerror 
%             errorbar(Amps,mean(squeeze(SSW_cismean(:,:,iExp)),2),(SSW_cismean(:,2,iExp)-SSW_cismean(:,1,iExp))./2,'k','LineWidth',3);
%         else
%             plot(Amps,mean(squeeze(SSW_cismean(:,:,iExp)),2),'k','LineWidth',3);
%         end
%     else
%         if showerror
%             errorbar(Amps,mean(squeeze(SSW_cismean(:,:,iExp)),2),(SSW_cismean(:,2,iExp)-SSW_cismean(:,1,iExp))./2,'k--','LineWidth',3);
%         else
%             plot(Amps,mean(squeeze(SSW_cismean(:,:,iExp)),2),'k--','LineWidth',3);
%         end    
%     end
% end
%     ax1=gca;
%     ax1.YColor='k';
%     %set(gcf,'Position',[114   249   564   320]);
%     %ax1.Position=pos;
%     ylim([0,20])
%     xlabel('MJO Amplitude (K/day)')
%     ylabel('#SSW per decade')
%     set(gca,'FontSize',16)
% if showUk1
%     ax2=axes('Position',ax1.Position,'YAxisLocation','right','Color','none','YColor','r');
%     hold on
% for iExp=1:NExps
%     if iExp==1
%         plot(Amps,Uk1smean(:,iExp),'r','LineWidth',2,'Parent',ax2);
%         hold on
%     else
%         plot(Amps,Uk1smean(:,iExp),'r--','LineWidth',2,'Parent',ax2);
%     end
% end
%      ylabel('zonal asym U k=1 amplitude')
% end
% 
% if Ndim==3
%     title('Group 1+2 #SSW (U pass 0)')
% else
%     title('Group 1+2 #SSW (U_t exceed 5m/s for 5 days)')
% end
% set(gca,'FontSize',16)

function [Uk1s,SSW_cis]=SSW_confidence(ExpName,Ndim,IDEAL,Amps,append)
N=length(Amps);
SSW_cis=zeros(N,2);
Uk1s=zeros(N,1)';
for icase=1:N
%% SSW    

% read data and match date
file1=strcat('./SSW_',ExpName,sprintf('%d',Amps(icase)),append,'.nc');
var='U';
offset=-22;
disp(file1)
SSW=[];
t=[];
for ifile=1:length(file1)
    SSW_temp=double(squeeze(ncread(char(file1(ifile)),var)));
    t_temp=ncread(char(file1(ifile)),'time');
    shift=offset;
    if length(t)~=0
        shift=max([offset,t(end)-t_temp(1)+1000]);
    end
    t_temp=t_temp+shift;
    t_temp=t_temp(1+11*365:60*365);
    SSW_temp=SSW_temp(1+11*365:60*365);
    SSW=[SSW;SSW_temp];
    t=[t;t_temp];
end
disp(size(SSW))
U_mean1(icase)=mean(SSW);


    
if IDEAL
    conseason=month(t)>0;
else
    conseason=(month(t)==1|month(t)==2|month(t)==12);
end

if Ndim==3
%find major SSW events
p_Urev=find(SSW<0&conseason);
nssw=0;
l_SSW=p_Urev.*0;
r_SSW=l_SSW;
if (length(p_Urev)>1)
    nssw=1;
    l_SSW(nssw)=p_Urev(1);
    r_SSW(nssw)=p_Urev(1);
    for i_urev=2:length(p_Urev)
        if (p_Urev(i_urev)-p_Urev(i_urev-1)<10)
            r_SSW(nssw)=p_Urev(i_urev);
        else
            nssw=nssw+1;
            l_SSW(nssw)=p_Urev(i_urev);
            r_SSW(nssw)=p_Urev(i_urev);
        end
    end
    l_SSW=l_SSW(1:nssw);
    r_SSW=r_SSW(1:nssw);
end
else % Ndim==2
    SSWsmt=smooth(SSW,5);
    SSWsmtdiff=diff(SSWsmt);
    [Slopemin,pSSW]=findpeaks(-SSWsmtdiff,'MinPeakDistance',40);
    l_SSW=pSSW(Slopemin>5.);
end
    SSW_series=SSW.*0;
    SSW_series(l_SSW)=1;
    
    %sampling
    yr=365;
    Nt=length(t);
    sample_size=16*yr;
    overlap=8*yr;
    Nsample=floor((Nt-sample_size)/(sample_size-overlap))+1;
    Nsample
    temp=zeros(Nsample,1);
    for isample=1:Nsample
        temp(isample)=mean(SSW_series(1+(isample-1)*(sample_size-overlap):sample_size+(isample-1)*(sample_size-overlap)))*365*10;
    end
    err=std(temp);
    SSW_cis(icase,:)=[mean(SSW_series)*365*10-err,mean(SSW_series)*365*10+err];
    disp(SSW_cis(icase,:))
    
    
%% stationary wave 
    file2=strcat('./All_',ExpName,sprintf('%d',Amps(icase)),append,'.nc');
    for ifile=1:length(file2)
    Utropo__=double(ncread(char(file2(ifile)),'U'));
    lat=double(squeeze(ncread(char(file2(ifile)),'lat')));
    lev=double(squeeze(ncread(char(file2(ifile)),'lev')));
    lon=double(squeeze(ncread(char(file2(ifile)),'lon')));
    lat0=45.;
    lev0=400.;
    [~,ilat]=min(abs(lat-lat0));
    [~,ilev]=min(abs(lev-lev0));
    r1=-5;
    r2=5;
    Utropolon=smoothn(squeeze(mean(Utropo__(:,ilat+r1:ilat+r2,ilev-1:ilev+1),3)),'robust');
    Uk1s(icase)=Uk1s(icase)+max(abs(sum(Utropolon.*repmat(exp(1i*lon/180*pi),[1,r2-r1+1]),1))/length(lon))/length(file2);
    end
    
    
end
end