%cd /Users/wanyingkang/Desktop/matlab/

IDEAL=true;
x_day_year=true;        %true: show year; false: show day
indx0=0;
%file_mjopc1='./WACCMSC_dmpdz2e-3_noQBO/MJO_PC_INDEX.nc';
%file_mjopc0='./WACCMSC_CONTROL_noQBO/MJO_PC_INDEX.nc';

%file1='./SSW_WA1Jan2D35_EQ1M40_5.nc';
%file0='./SSW_WA1Jan2D35_FREE.nc';
%file1='./SSW_WA1Jan35_EQ1M40lon60180_5.nc';
%file0='./SSW_WA1Jan35_FREE.nc';
%file1='./SSW_WA1Jan35_EQ1M40_5.nc';
%file0='./SSW_WA1Jan35_FREE.nc';
%file1='./SSW_WA1Jan35_EQ1M40_2.nc';
%file0='./SSW_WA1Jan35_FREE.nc';
%file1='./SSW_WA1Jan35_EQ1M40lon60180_30_2.nc';
%file0='./SSW_WA1Jan35_EQ1M40_0.nc';

file1='./SSW_WA1Jan35_EQ1M40lon60180_5.nc';
file0='./SSW_WA1Jan35_EQ1M40_0.nc';
%file1='./SSW_IDEAL_WA1Jan35_EQ1M40lon180300_5_2.nc';
%file0='./SSW_WA1Jan35_EQ1M40_0.nc';
%file1='./SSW_IDEAL_WA1Jan35_EQ1M40lon30060_5_2.nc';
%file0='./SSW_WA1Jan35_EQ1M40_0.nc';

lag=0; %only for ~IDEAL
Ndim=3;
var='U';
offset=-4; 
if ~IDEAL
    offset=-1;
end

SSW_WA1Jan_EQ1M40=double(squeeze(ncread(file1,var)));
SSW_WA1Jan_FREE=double(squeeze(ncread(file0,var)));
t1=ncread(file1,'time');

t1=t1(1478:1478+50*365);
SSW_WA1Jan_EQ1M40=SSW_WA1Jan_EQ1M40(1478:1478+50*365);

%t1=t1-offset;
t1=t1-offset-37*40;
t0=ncread(file0,'time');
%t0=t0-offset
t0=t0-offset-37*40;


i_beg=find(t0==t1(1));
i_end=find(t0==t1(length(t1)));
t0=t0(i_beg:i_end);

t=t1;
SSW_WA1Jan_FREE=SSW_WA1Jan_FREE(i_beg:i_end);

%% Time Series
figure1=figure(20);
clf;
% Create axes
axes1 = axes('Parent',figure1,'FontSize',16);
xlim([min(t),max(t)])
ylim([-20.,80.])
box(axes1,'on');
hold(axes1,'on');

plot(t,SSW_WA1Jan_FREE,'b-')
plot(t,SSW_WA1Jan_EQ1M40,'r-')
hold on
%plot(0)
plot(t,SSW_WA1Jan_FREE.*0,'k--')

%peak finding
if Ndim==3
[Umin_EQ1M40,pSSW_EQ1M40]=findpeaks(-SSW_WA1Jan_EQ1M40,'MinPeakDistance',100);
if ~IDEAL
    pSSW_EQ1M40=pSSW_EQ1M40(Umin_EQ1M40>0&(month(t(pSSW_EQ1M40))==1|month(t(pSSW_EQ1M40))==2|month(t(pSSW_EQ1M40))==12));
else
    pSSW_EQ1M40=pSSW_EQ1M40(Umin_EQ1M40>0);
end
[Umin_FREE,pSSW_FREE]=findpeaks(-SSW_WA1Jan_FREE,'MinPeakDistance',100);
if ~IDEAL
    pSSW_FREE=pSSW_FREE(Umin_FREE>0&(month(t(pSSW_FREE))==1|month(t(pSSW_FREE))==2|month(t(pSSW_FREE))==12));
else
    pSSW_FREE=pSSW_FREE(Umin_FREE>0);
end
else
SSWsmt_WA1Jan_EQ1M40=smooth(SSW_WA1Jan_EQ1M40,5);
SSWsmt_WA1Jan_FREE=smooth(SSW_WA1Jan_FREE,5);
SSWsmtdiff_WA1Jan_EQ1M40=diff(SSWsmt_WA1Jan_EQ1M40);
SSWsmtdiff_WA1Jan_FREE=diff(SSWsmt_WA1Jan_FREE);
[Slopemin_EQ1M40,pSSW_EQ1M40]=findpeaks(-SSWsmtdiff_WA1Jan_EQ1M40,'MinPeakDistance',100);
[Slopemin_FREE,pSSW_FREE]=findpeaks(-SSWsmtdiff_WA1Jan_FREE,'MinPeakDistance',100);
pSSW_EQ1M40=pSSW_EQ1M40(Slopemin_EQ1M40>5.);
pSSW_FREE=pSSW_FREE(Slopemin_FREE>5.);
end
    disp(strcat('Using peak automatic diagose, MJO forced #SSW:',sprintf('%d',length(pSSW_EQ1M40)),'/',sprintf('%d',length(SSW_WA1Jan_EQ1M40)/365),' Years'))
    disp(strcat('Using peak automatic diagose, CONTROL #SSW:',sprintf('%d',length(pSSW_FREE)),'/',sprintf('%d',length(SSW_WA1Jan_FREE)/365),' Years'))
%    hold on
%if length(pSSW_EQ1M40)>0
%    plot(pSSW_EQ1M40+t(1),-20.,'r^','markerfacecolor',[1 0 0])
%end
%if length(pSSW_FREE)>0 
%plot(pSSW_FREE+t(1),-20.,'b^')
%end


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

if Ndim==3
p1_Urev=find(SSW_WA1Jan_EQ1M40<0&conseason);
nssw1=0;
l1_SSW=p1_Urev.*0;
r1_SSW=l1_SSW;
if (length(p1_Urev)>1)
    nssw1=1;
    l1_SSW(nssw1)=p1_Urev(1);
    r1_SSW(nssw1)=p1_Urev(1);
    for i_urev=2:length(p1_Urev)
        if (p1_Urev(i_urev)-p1_Urev(i_urev-1)<20)
            r1_SSW(nssw1)=p1_Urev(i_urev);
        else
            nssw1=nssw1+1;
            l1_SSW(nssw1)=p1_Urev(i_urev);
            r1_SSW(nssw1)=p1_Urev(i_urev);
        end
    end
    l1_SSW=l1_SSW(1:nssw1);
    r1_SSW=r1_SSW(1:nssw1);
end

p0_Urev=find(SSW_WA1Jan_FREE<0&conseason);
nssw0=0;
l0_SSW=p0_Urev.*0;
r0_SSW=l0_SSW;
if (length(p0_Urev)>1)
    nssw0=1;
    l0_SSW(nssw0)=p0_Urev(1);
    r0_SSW(nssw0)=p0_Urev(1);
    for i_urev=2:length(p0_Urev)
        if (p0_Urev(i_urev)-p0_Urev(i_urev-1)<20)
            r0_SSW(nssw0)=p0_Urev(i_urev);
        else
            nssw0=nssw0+1;
            l0_SSW(nssw0)=p0_Urev(i_urev);
            r0_SSW(nssw0)=p0_Urev(i_urev);
        end
    end
    l0_SSW=l0_SSW(1:nssw0);
    r0_SSW=r0_SSW(1:nssw0);
end

%USE WYKANG DIAGNOSTIC to IDENTIFY SSW POSTIION
%pSSW_EQ1M40=l1_SSW;
%pSSW_FREE=l0_SSW;

%output and plot peak triangles
disp(strcat('MJO forced #SSW:',sprintf('%d',nssw1),'/',sprintf('%d',length(SSW_WA1Jan_EQ1M40)/365),' Years'))
disp(strcat('CONTROL #SSW:',sprintf('%d',nssw0),'/',sprintf('%d',length(SSW_WA1Jan_FREE)/365),' Years'))
end

hold on
if length(pSSW_EQ1M40)>0
    plot(pSSW_EQ1M40+t(1),-20.,'r^','markerfacecolor',[1 0 0])
end
if length(pSSW_FREE)>0 
plot(pSSW_FREE+t(1),-20.,'b^')
end

if ~x_day_year
    xlabel('Days','FontSize',16)
else
    xlabel('Year','FontSize',16)
    datetick('x','yy')
end

ylabel('U 70N 10mb (m/s)','FontSize',16)
set(gcf, 'Position', [0 0 900 300]);

%% histogram plot for SSW events at different MJO phases
tSSW_EQ1M40=pSSW_EQ1M40+t(1)-1;
figure(21)
clf
if ~IDEAL
  phaCnt       = [202.5, 247.5, 292.5, 337.5,  22.5,  67.5, 112.5, 157.5];
  nPhase      = length(phaCnt);
  phaBnd      = 180.0/nPhase;
  angBnd      = [phaCnt;phaCnt];
  angBnd(1,:) = phaCnt-phaBnd;
  angBnd(2,:) = phaCnt+phaBnd;
end

if IDEAL
    phase1=ceil((mod(t+20,40)+1)/5); %+20 is new
    phase0=ceil((mod(t+20,40)+1)/5); %+20 is new
else
    pc1_EQ1M40=ncread(file_mjopc1,'PC1');
    pc2_EQ1M40=ncread(file_mjopc1,'PC2');
    indx_EQ1M40=ncread(file_mjopc1,'MJO_INDEX');
    timemjo_EQ1M40=ncread(file_mjopc1,'time');
    need1=find(timemjo_EQ1M40>=min(t)-lag&timemjo_EQ1M40<=max(t)-lag);
    timemjo_EQ1M40=timemjo_EQ1M40(need1);
    indx_EQ1M40=indx_EQ1M40(need1);
    pc1_EQ1M40=pc1_EQ1M40(need1);
    pc2_EQ1M40=pc2_EQ1M40(need1);
    
    pc1_FREE=ncread(file_mjopc0,'PC1');
    pc2_FREE=ncread(file_mjopc0,'PC2');
    indx_FREE=ncread(file_mjopc0,'MJO_INDEX');
    timemjo_FREE=ncread(file_mjopc0,'time');
    need1=find(timemjo_FREE>=min(t)-lag&timemjo_FREE<=max(t)-lag);
    timemjo_FREE=timemjo_FREE(need1);
    indx_FREE=indx_FREE(need1);
    pc1_FREE=pc1_FREE(need1);
    pc2_FREE=pc2_FREE(need1);
    
    ang_EQ1M40         = atan2(pc2_EQ1M40,pc1_EQ1M40)*180./pi;
    ang_EQ1M40(ang_EQ1M40<0)     = ang_EQ1M40(ang_EQ1M40<0) + 360;
    phase1=t.*0;
    for iphase=1:nPhase
        phase1(ang_EQ1M40>angBnd(1,iphase)&ang_EQ1M40<angBnd(2,iphase)&indx_EQ1M40>indx0)=iphase;
    end
    
    ang_FREE = atan2(pc2_FREE,pc1_FREE)*180./pi;
    ang_FREE(ang_FREE<0)  = ang_FREE(ang_FREE<0) + 360;
    phase0=t.*0;
    for iphase=1:nPhase
        phase0(ang_FREE>angBnd(1,iphase)&ang_FREE<angBnd(2,iphase)&indx_FREE>indx0)=iphase;
    end
end
dist1=phase1(pSSW_EQ1M40);
dist0=phase0(pSSW_FREE);

if IDEAL
    hist1=histogram(dist1);    
    set(hist1,'FaceColor',[1 0.4 0]);
    xlim([0.5 8.5])
    hold on
    %hist0=histogram(dist0);
    %set(hist0,'FaceColor','none','EdgeColor','k','LineWidth',3);
    %plot([1:8],repmat(length(dist0)/8,8),'k*')
    plot([0.5,8.5],[length(dist0)/8,length(dist0)/8],'k-','LineWidth',3)
    [count1,ed1]=histcounts(phase1,[0.5:1:8.5]);
    count1_normssw=count1.*length(dist1)./length(phase1);
    [count0,ed0]=histcounts(phase0,[0.5:1:8.5]);
    count0_normssw=count0.*length(dist0)./length(phase0);
    plot([1:8],count0_normssw,'k*')
    plot([1:8],count1_normssw,'r*')
else
    hist1=histogram(dist1);    
    set(hist1,'FaceColor',[1 0.4 0]);
    xlim([-0.5 8.5])
    hold on
    hist0=histogram(dist0);
    set(hist0,'FaceColor','none','EdgeColor','k','LineWidth',3);
    %[sswcount0,ed0]=histcounts(dist0,[-0.5:1:8.5]);
    %plot([0:8],sswcount0,'k*')
    [count1,ed1]=histcounts(phase1,[-0.5:1:8.5]);
    count1_normssw=count1.*length(dist1)./length(phase1);
    [count0,ed0]=histcounts(phase0,[-0.5:1:8.5]);
    count0_normssw=count0.*length(dist0)./length(phase0);
    plot([0:8],count0_normssw,'k*')
    plot([0:8],count1_normssw,'r*')
end

set(gca, 'FontSize',18)
xlabel('MJO Phase')
ylabel('# of SSW events')
%bar([phase0,phase1])

[f1,x1]=ecdf(phase1);
[h1,p1,ks1,cv1] = kstest(dist1,'CDF',[x1(2:end),f1(2:end)]);
if p1<0.05
    disp(strcat('MJO forced model Significantly different from uniform distribution. p=',sprintf('%3.2f',p1)));
else
    disp(strcat('MJO forced model Not significantly different from uniform distribution. p=',sprintf('%3.2f',p1)));
end
[f0,x0]=ecdf(phase0);
[h0,p0,ks0,cv0] = kstest(dist0,'CDF',[x0(2:end),f0(2:end)]);
if p0<0.05
    disp(strcat('CONTROL model Significantly different from uniform distribution. p=',sprintf('%3.2f',p0)));
else
    disp(strcat('CONTROL model Not significantly different from uniform distribution. p=',sprintf('%3.2f',p0)));
end
%th=text(5,11,strcat('CONTROL:    p=',sprintf('%3.2f',p0)));
%set(th,'Color','k','FontSize',18)
th=text(0.45,0.9,strcat('P-Value for Uniform Dist: ',sprintf('%3.2f',p1)),'Units','normalized','FontSize',18,'Color','r');
%,'FontSize',18,'Position',[0.4 0.8],'Units','normalized'
%histogram plot of U wind
if ~IDEAL
    indDJF=find(month(t)==1|month(t)==12|month(t)==2);
    SSW_WA1Jan_EQ1M40=SSW_WA1Jan_EQ1M40(indDJF);
    SSW_WA1Jan_FREE=SSW_WA1Jan_FREE(indDJF);
end
    
%% U wind distribution catagorized by MJO phases
figure(22)
clf
edges=[0:5:60];
edges=[-20,edges];
hist3=histogram(SSW_WA1Jan_EQ1M40,edges);
hold on
set(hist3,'FaceColor',[1 0 0],'EdgeColor','none','Orientation','horizontal')
hist4=histogram(SSW_WA1Jan_FREE,edges);
set(hist4,'FaceColor','none','EdgeColor',[0 0 1],'Orientation','horizontal')
xlabel('#')
ylabel('U wind (m/s)')
ylim([-20 60])
set(gca,'FontSize',16)
set(gcf, 'Position', [850 15 300 285]);

%bar3 plot of U wind
figure(23)
clf
counts1=zeros(8,length(edges)-1);
counts0=zeros(8,length(edges)-1);
for bin=1:length(edges)-1
    temp1=find(SSW_WA1Jan_EQ1M40>edges(bin)&SSW_WA1Jan_EQ1M40<edges(bin+1));
    for phase=1:8
        counts1(phase,bin)=sum(phase1(temp1)==phase);
    end
end
for bin=1:length(edges)-1
    temp0=find(SSW_WA1Jan_FREE>edges(bin)&SSW_WA1Jan_FREE<edges(bin+1));
    %temp0=ceil((mod(temp0+t(1)-1+20,40)+1)/5);% +20 new
    for phase=1:8
        counts0(phase,bin)=sum(phase0(temp0)==phase);
    end
end
%bar3(edges(2:end),counts1')
barh((edges(1:end-1)+edges(2:end))./2,counts1','stacked')
hold on
barh((edges(1:end-1)+edges(2:end))./2,sum(counts0),'stacked','LineWidth',2,'FaceColor','none')
barh((edges(1:end-1)+edges(2:end))./2,mean(counts0),'stacked','LineWidth',2,'FaceColor','none','EdgeColor',[0 1 0])
set(gca,'FontSize',16)
xlabel('#')
ylabel('U (m/s)')
set(gcf, 'Position', [850 15 300 285]);

figure(24)
clf
bar3(edges(2:end),[mean(counts0)',counts1'])

