%% clc;close all;clear; mainpath = 'C:\Users\l1989\Desktop\PolarVortex\Epm\'; cd(mainpath); load('Epm2011-2019final.mat');load('EpvN2_2011-2019final.mat'); timesort_datetime = datetime(timesort,'convertfrom','datenum'); timesortTotal = timesort; n = 1; l = size(gwpedErr,1); for nYr = 2011:1:2019 for nMon = 1:1:12 index = (timesort_datetime.Year == nYr & timesort_datetime.Month == nMon); if n == 1 mean30to50EpmTotal = nanmean(gwped(:,index),1); ErrorMean30to50EpmTotal = sqrt(sum(gwpedErr(:,index).^2,1))/l; else mean30to50EpmTotal = [mean30to50EpmTotal nanmean(gwped(:,index),1)]; ErrorMean30to50EpmTotal = [ErrorMean30to50EpmTotal sqrt(sum(gwpedErr(:,index).^2,1))/l]; end n = n+1; end end dn_start = datenum(2010,12,31); TimeSeries9yrs = datenum(2011,1,1):1:datenum(2019,12,31); dtSeries9yrs = datetime(TimeSeries9yrs,'convertfrom','datenum'); TimeSeries9yrs = TimeSeries9yrs - dn_start; %% fit each year (for plotting) % eq = fittype('a + b*cos(2*pi/365*(x-c))+d*cos(4*pi/365*(x-e))'); % ft = cell(9,1); % dataMx = cell(9,1); % dataMy = cell(9,1); % for nYr = 2011:1:2019 % index = find(timesort_datetime.Year == nYr); % dataMx{nYr-2010} = timesortTotal(index); % dataMy{nYr-2010} = mean30to50EpmTotal(index)'; % ft{nYr-2010} = fit(timesortTotal(index),mean30to50EpmTotal(index)',eq); % i = i+1; % end ft = cell(9,1); dataMx = cell(9,1); dataMy = cell(9,1); % get weight AnnualSegCnt = zeros(9,1); MonthlySegCnt = zeros(9,12); ftWeight = cell(9,1); for nYr = 2011:1:2019 AnnualSegCnt(nYr-2010) = length(find(timesort_datetime.Year == nYr)~=0); tempWeight = []; for nMon = 1:1:12 MonthlySegCnt(nYr-2010,nMon) = length(find(timesort_datetime.Year == nYr & timesort_datetime.Month == nMon)~=0); tempWeight{nMon} = ones(MonthlySegCnt(nYr-2010,nMon),1)*(1./MonthlySegCnt(nYr-2010,nMon)); end ftWeight{nYr-2010} = [tempWeight{1};tempWeight{2};tempWeight{3};tempWeight{4};tempWeight{5};tempWeight{6};tempWeight{7};tempWeight{8};tempWeight{9};tempWeight{10};tempWeight{11};tempWeight{12}]; end for nYr = 2011:1:2019 index = find(timesort_datetime.Year == nYr); dataMx{nYr-2010} = timesortTotal(index)-dn_start; % temp = day(datetime(dataMx{nYr-2010},'convertfrom','datenum'),'dayofyear'); dataMy{nYr-2010} = mean30to50EpmTotal(index)'; dataN2{nYr-2010} = N2(index)'; dataEpv{nYr-2010} = Epv(index)'; beta0 = [1,1,180,1,180]; ft{nYr-2010} = fitnlm(dataMx{nYr-2010},mean30to50EpmTotal(index)','y ~ b1 + b2.*cos((2*pi/365).*(x1-b3))+b4.*cos((4*pi/365).*(x1-b5))',beta0); ftN2{nYr-2010} = fitnlm(dataMx{nYr-2010},N2(index)','y ~ b1 + b2.*cos((2*pi/365).*(x1-b3))+b4.*cos((4*pi/365).*(x1-b5))',beta0); ftEpv{nYr-2010} = fitnlm(dataMx{nYr-2010},Epv(index)','y ~ b1 + b2.*cos((2*pi/365).*(x1-b3))+b4.*cos((4*pi/365).*(x1-b5))',beta0); i = i+1; end %% fitting 9 years temp = timesortTotal-dn_start; ft9yrs = fitnlm(temp,mean30to50EpmTotal','y ~ b1 + b2.*cos((2*pi/365).*(x1-b3))+b4.*cos((4*pi/365).*(x1-b5))',beta0); ft9yrsN2 = fitnlm(temp,N2','y ~ b1 + b2.*cos((2*pi/365).*(x1-b3))+b4.*cos((4*pi/365).*(x1-b5))',beta0); ft9yrsEpv = fitnlm(temp,Epv','y ~ b1 + b2.*cos((2*pi/365).*(x1-b3))+b4.*cos((4*pi/365).*(x1-b5))',beta0); ft911 = fitnlm(temp,mean30to50EpmTotal','y ~ b1 + b2.*cos(b3*x1-b4)',[1,1,600,1]); %% Epm annual and uniform fitting - figure 1(a)(b) scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3),scrsz(4)*0.7]); subplot('position',[0.07,0.58,0.9,0.33]); p1 = errorbar(timesortTotal-dn_start,mean30to50EpmTotal,ErrorMean30to50EpmTotal,'*b'); ylim([0,20]);xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'YMinorTick','on');set(gca,'TickDir','out');set(gca,'TickLength',[0.006 0.0025]); ylabel('E_{pm} (J/kg)','FontSize',20); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',15); text(datenum('20110101','yyyymmdd')-dn_start+110,-4.8,'2011','fontsize',20);text(datenum('20120101','yyyymmdd')-dn_start+110,-4.8,'2012','fontsize',20); text(datenum('20130101','yyyymmdd')-dn_start+110,-4.8,'2013','fontsize',20);text(datenum('20140101','yyyymmdd')-dn_start+110,-4.8,'2014','fontsize',20); text(datenum('20150101','yyyymmdd')-dn_start+110,-4.8,'2015','fontsize',20);text(datenum('20160101','yyyymmdd')-dn_start+110,-4.8,'2016','fontsize',20); text(datenum('20170101','yyyymmdd')-dn_start+110,-4.8,'2017','fontsize',20);text(datenum('20180101','yyyymmdd')-dn_start+110,-4.8,'2018','fontsize',20); text(datenum('20190101','yyyymmdd')-dn_start+110,-4.8,'2019','fontsize',20); for nYr = 2011:1:2019 hold on; dnAll = (datenum(nYr,1,1)-dn_start):1:(datenum(nYr,12,31)-dn_start); p2 = plot(dnAll,ft{nYr-2010}.Coefficients.Estimate(1) + ft{nYr-2010}.Coefficients.Estimate(2).*cos((2*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(3)))+ft{nYr-2010}.Coefficients.Estimate(4).*cos((4*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(5))), 'r','linewidth',2); end set(gca,'fontsize',16);title('Altitude-mean Gravity Wave Epm from 30 to 50 km'); % monthly mean Epm % for nYr = 2011:1:2019 % for month = 1:1:12 % ind = find(timesort_datetime.Year == nYr & timesort_datetime.Month == month); % EpmMonthlyMean(month, nYr-2010) = mean(mean30to50EpmTotal(ind)); % scatter(datenum(nYr,month,15)-dn_start,EpmMonthlyMean(month, nYr-2010),80); % end % end subplot('position',[0.07,0.12,0.9,0.33]); errorbar(timesortTotal-dn_start,mean30to50EpmTotal,ErrorMean30to50EpmTotal,'*b');hold on; ylim([0,20]);xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'YMinorTick','on');set(gca,'TickDir','out');set(gca,'TickLength',[0.006 0.0025]); ylabel('E_{pm} (J/kg)','FontSize',20); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',15); text(datenum('20110101','yyyymmdd')-dn_start+110,-4.8,'2011','fontsize',20);text(datenum('20120101','yyyymmdd')-dn_start+110,-4.8,'2012','fontsize',20); text(datenum('20130101','yyyymmdd')-dn_start+110,-4.8,'2013','fontsize',20);text(datenum('20140101','yyyymmdd')-dn_start+110,-4.8,'2014','fontsize',20); text(datenum('20150101','yyyymmdd')-dn_start+110,-4.8,'2015','fontsize',20);text(datenum('20160101','yyyymmdd')-dn_start+110,-4.8,'2016','fontsize',20); text(datenum('20170101','yyyymmdd')-dn_start+110,-4.8,'2017','fontsize',20);text(datenum('20180101','yyyymmdd')-dn_start+110,-4.8,'2018','fontsize',20); text(datenum('20190101','yyyymmdd')-dn_start+110,-4.8,'2019','fontsize',20); plot(TimeSeries9yrs, ft9yrs.Coefficients.Estimate(1)+ft9yrs.Coefficients.Estimate(2).*cos((2*pi/365).*(TimeSeries9yrs-ft9yrs.Coefficients.Estimate(3)))+ft9yrs.Coefficients.Estimate(4).*cos((4*pi/365).*(TimeSeries9yrs-ft9yrs.Coefficients.Estimate(5))),'r','linewidth',2); set(gca,'fontsize',16); %% Epv annual and uniform fitting - figure 1(c)(d) scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3),scrsz(4)*0.7]); subplot('position',[0.07,0.58,0.9,0.33]); errorbar(timesortTotal-dn_start,Epv,EpvErr,'*b'); xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'YMinorTick','on');set(gca,'TickDir','out');set(gca,'TickLength',[0.006 0.0025]); ylabel('E_{pv} (J/m^3)','FontSize',20);ylim([0 0.05]); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',16); text(datenum('20110101','yyyymmdd')-dn_start+110,-0.012,'2011','fontsize',20);text(datenum('20120101','yyyymmdd')-dn_start+110,-0.012,'2012','fontsize',20); text(datenum('20130101','yyyymmdd')-dn_start+110,-0.012,'2013','fontsize',20);text(datenum('20140101','yyyymmdd')-dn_start+110,-0.012,'2014','fontsize',20); text(datenum('20150101','yyyymmdd')-dn_start+110,-0.012,'2015','fontsize',20);text(datenum('20160101','yyyymmdd')-dn_start+110,-0.012,'2016','fontsize',20); text(datenum('20170101','yyyymmdd')-dn_start+110,-0.012,'2017','fontsize',20);text(datenum('20180101','yyyymmdd')-dn_start+110,-0.012,'2018','fontsize',20); text(datenum('20190101','yyyymmdd')-dn_start+110,-0.012,'2019','fontsize',20);hold on; for nYr = 2011:1:2019 hold on; dnAll = (datenum(nYr,1,1)-dn_start):1:(datenum(nYr,12,31)-dn_start); p2 = plot(dnAll,ftEpv{nYr-2010}.Coefficients.Estimate(1) + ftEpv{nYr-2010}.Coefficients.Estimate(2).*cos((2*pi/365).*(dnAll-ftEpv{nYr-2010}.Coefficients.Estimate(3)))+ftEpv{nYr-2010}.Coefficients.Estimate(4).*cos((4*pi/365).*(dnAll-ftEpv{nYr-2010}.Coefficients.Estimate(5))), 'r','linewidth',2); end title('Altitude-mean Gravity Wave Epv from 30 to 50 km'); subplot('position',[0.07,0.12,0.9,0.33]); errorbar(timesortTotal-dn_start,Epv,EpvErr,'*b'); xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'YMinorTick','on');set(gca,'TickDir','out');set(gca,'TickLength',[0.006 0.0025]); ylabel('E_{pv} (J/m^3)','FontSize',20);ylim([0 0.05]); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',16); text(datenum('20110101','yyyymmdd')-dn_start+110,-0.012,'2011','fontsize',20);text(datenum('20120101','yyyymmdd')-dn_start+110,-0.012,'2012','fontsize',20); text(datenum('20130101','yyyymmdd')-dn_start+110,-0.012,'2013','fontsize',20);text(datenum('20140101','yyyymmdd')-dn_start+110,-0.012,'2014','fontsize',20); text(datenum('20150101','yyyymmdd')-dn_start+110,-0.012,'2015','fontsize',20);text(datenum('20160101','yyyymmdd')-dn_start+110,-0.012,'2016','fontsize',20); text(datenum('20170101','yyyymmdd')-dn_start+110,-0.012,'2017','fontsize',20);text(datenum('20180101','yyyymmdd')-dn_start+110,-0.012,'2018','fontsize',20); text(datenum('20190101','yyyymmdd')-dn_start+110,-0.012,'2019','fontsize',20);hold on; set(gca,'fontsize',16); plot(TimeSeries9yrs, ft9yrsEpv.Coefficients.Estimate(1)+ft9yrsEpv.Coefficients.Estimate(2).*cos((2*pi/365).*(TimeSeries9yrs-ft9yrsEpv.Coefficients.Estimate(3)))+ft9yrsEpv.Coefficients.Estimate(4).*cos((4*pi/365).*(TimeSeries9yrs-ft9yrsEpv.Coefficients.Estimate(5))),'r','linewidth',2); %% N2 - figure 2 f = figure('position',[1,300,scrsz(3),scrsz(4)*0.7]); subplot('position',[0.07,0.58,0.9,0.33]); errorbar(timesortTotal-dn_start,N2*10^4,N2Err*10^4,'*b'); xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'YMinorTick','on');set(gca,'TickDir','out');set(gca,'fontsize',20);set(gca,'TickLength',[0.006 0.0025]); ylabel('N^2 (10^{-4} S^{-2})','FontSize',20);ylim([3 6]);yticks([3 4 5 6]); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',15); text(datenum('20110101','yyyymmdd')-dn_start+110,2.3,'2011','fontsize',20);text(datenum('20120101','yyyymmdd')-dn_start+110,2.3,'2012','fontsize',20); text(datenum('20130101','yyyymmdd')-dn_start+110,2.3,'2013','fontsize',20);text(datenum('20140101','yyyymmdd')-dn_start+110,2.3,'2014','fontsize',20); text(datenum('20150101','yyyymmdd')-dn_start+110,2.3,'2015','fontsize',20);text(datenum('20160101','yyyymmdd')-dn_start+110,2.3,'2016','fontsize',20); text(datenum('20170101','yyyymmdd')-dn_start+110,2.3,'2017','fontsize',20);text(datenum('20180101','yyyymmdd')-dn_start+110,2.3,'2018','fontsize',20); text(datenum('20190101','yyyymmdd')-dn_start+110,2.3,'2019','fontsize',20); for nYr = 2011:1:2019 hold on; dnAll = (datenum(nYr,1,1)-dn_start):1:(datenum(nYr,12,31)-dn_start); p2 = plot(dnAll,(ftN2{nYr-2010}.Coefficients.Estimate(1) + ftN2{nYr-2010}.Coefficients.Estimate(2).*cos((2*pi/365).*(dnAll-ftN2{nYr-2010}.Coefficients.Estimate(3)))+ftN2{nYr-2010}.Coefficients.Estimate(4).*cos((4*pi/365).*(dnAll-ftN2{nYr-2010}.Coefficients.Estimate(5))))*10^4, 'r','linewidth',2); end title('Altitude-mean Gravity Wave Epv from 30 to 50 km','fontsize',24); subplot('position',[0.07,0.12,0.9,0.33]); errorbar(timesortTotal-dn_start,N2*10^4,N2Err*10^4,'*b');hold on; ylim([3,6]);yticks([3 4 5 6]);xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'YMinorTick','on');set(gca,'TickDir','out');set(gca,'TickLength',[0.006 0.0025]); ylabel('N^2 (10^{-4} S^{-2})','FontSize',20);set(gca,'fontsize',20); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',15); text(datenum('20110101','yyyymmdd')-dn_start+110,2.3,'2011','fontsize',20);text(datenum('20120101','yyyymmdd')-dn_start+110,2.3,'2012','fontsize',20); text(datenum('20130101','yyyymmdd')-dn_start+110,2.3,'2013','fontsize',20);text(datenum('20140101','yyyymmdd')-dn_start+110,2.3,'2014','fontsize',20); text(datenum('20150101','yyyymmdd')-dn_start+110,2.3,'2015','fontsize',20);text(datenum('20160101','yyyymmdd')-dn_start+110,2.3,'2016','fontsize',20); text(datenum('20170101','yyyymmdd')-dn_start+110,2.3,'2017','fontsize',20);text(datenum('20180101','yyyymmdd')-dn_start+110,2.3,'2018','fontsize',20); text(datenum('20190101','yyyymmdd')-dn_start+110,2.3,'2019','fontsize',20); plot(TimeSeries9yrs, (ft9yrsN2.Coefficients.Estimate(1)+ft9yrsN2.Coefficients.Estimate(2).*cos((2*pi/365).*(TimeSeries9yrs-ft9yrsN2.Coefficients.Estimate(3)))+ft9yrsN2.Coefficients.Estimate(4).*cos((4*pi/365).*(TimeSeries9yrs-ft9yrsN2.Coefficients.Estimate(5))))*10^4,'r','linewidth',2); %% composite fitting - figure 3 n = 1; CompositeEpm = cell(365,1); CompositeN2 = cell(365,1); CompositeEpv = cell(365,1); for nMon = 1:1:12 for nDay = 1:1:eomday(2011,nMon) Ind = timesort_datetime.Month == nMon & timesort_datetime.Day == nDay; CompositeEpm{n} = mean30to50EpmTotal(Ind); CompositeN2{n} = N2(Ind)'; CompositeEpv{n} = Epv(Ind)'; n = n+1; end end % take mean for every 5 days n = 1; binDay = 5; yData = zeros(ceil(365/binDay),1); xData = zeros(ceil(365/binDay),1); yN2Data = zeros(ceil(365/binDay),1); yEpvData = zeros(ceil(365/binDay),1); for nCnt = 1:binDay:365 yData(n) = nanmean([CompositeEpm{1+(n-1)*binDay} CompositeEpm{2+(n-1)*binDay} CompositeEpm{3+(n-1)*binDay} CompositeEpm{4+(n-1)*binDay} CompositeEpm{5+(n-1)*binDay}]); yN2Data(n) = nanmean([CompositeN2{1+(n-1)*binDay} CompositeN2{2+(n-1)*binDay} CompositeN2{3+(n-1)*binDay} CompositeN2{4+(n-1)*binDay} CompositeN2{5+(n-1)*binDay}]); yEpvData(n) = nanmean([CompositeEpv{1+(n-1)*binDay} CompositeEpv{2+(n-1)*binDay} CompositeEpv{3+(n-1)*binDay} CompositeEpv{4+(n-1)*binDay} CompositeEpv{5+(n-1)*binDay}]); xData(n) = nCnt+2; n = n+1; end % fitting ftComposite = fitnlm(xData,yData','y ~ b1 + b2.*cos((2*pi/365).*(x1-b3))+b4.*cos((4*pi/365).*(x1-b5))',beta0); ftN2Composite = fitnlm(xData,yN2Data','y ~ b1 + b2.*cos((2*pi/365).*(x1-b3))+b4.*cos((4*pi/365).*(x1-b5))',beta0); ftEpvComposite = fitnlm(xData,yEpvData','y ~ b1 + b2.*cos((2*pi/365).*(x1-b3))+b4.*cos((4*pi/365).*(x1-b5))',beta0); dn = 1:1:365; UniformFit = ftComposite.Coefficients.Estimate(1)+ftComposite.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ftComposite.Coefficients.Estimate(3)))+ftComposite.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ftComposite.Coefficients.Estimate(5))); figure('position',[1,300,scrsz(3)*0.28,scrsz(4)]); subplot('position',[0.158,0.705,0.8,0.26]); plot(xData,yData,'*b','MarkerSize',8);hold on; xlim([1 365]);ylim([0 12]); set(gca,'YMinorTick','on');set(gca,'TickDir','out');set(gca,'box','on'); set(gca,'xtick',[1,32,60,91,121,152,182,213,244,275,305,335],'XTickLabel',{}); set(gca,'ytick',[0,2,4,6,8,10,12,14,16],'fontsize',18); text([1,32,60,91,121,152,182,213,244,275,305,335]+10,[-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1]+0.4,['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D'],'fontsize',18); ylabel('Epm (J/kg)','FontSize',18);%xlabel('Month','fontsize',20); dn = 1:1:365; UniformFit = ftComposite.Coefficients.Estimate(1)+ftComposite.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ftComposite.Coefficients.Estimate(3)))+ftComposite.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ftComposite.Coefficients.Estimate(5))); plot(1:1:365,UniformFit,'r','linewidth',2); %legend('E_p_m', 'Fitting','fontsize',16); title('9-year Composite Epm (30-50 km)'); subplot('position',[0.158,0.385,0.8,0.26]); plot(xData,yEpvData*10^2,'*b','MarkerSize',8);hold on; xlim([1 365]);ylim([0 3]); set(gca,'YMinorTick','on');set(gca,'TickDir','out');set(gca,'box','on'); set(gca,'xtick',[1,32,60,91,121,152,182,213,244,275,305,335],'XTickLabel',{}); set(gca,'ytick',[0,1,2,3],'fontsize',18); text([1,32,60,91,121,152,182,213,244,275,305,335]+10,[-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1]+0.83,['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D'],'fontsize',18); ylabel('Epv (J/m^3)','FontSize',18);%xlabel('Month','fontsize',20); dn = 1:1:365; UniformFit = ftEpvComposite.Coefficients.Estimate(1)+ftEpvComposite.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ftEpvComposite.Coefficients.Estimate(3)))+ftEpvComposite.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ftEpvComposite.Coefficients.Estimate(5))); plot(1:1:365,UniformFit*10^2,'r','linewidth',2); %legend('E_p_v', 'Fitting','fontsize',18); title('9-year Composite Epv (30-50 km)'); subplot('position',[0.158,0.055,0.8,0.26]); plot(xData,yN2Data*10^4,'*b','MarkerSize',8);hold on; xlim([1 365]);ylim([3 6]); set(gca,'YMinorTick','on');set(gca,'TickDir','out');set(gca,'box','on'); set(gca,'xtick',[1,32,60,91,121,152,182,213,244,275,305,335],'XTickLabel',{}); set(gca,'ytick',[3 4 5 6],'fontsize',18); text([1,32,60,91,121,152,182,213,244,275,305,335]+10,[-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1]+3.83,['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D'],'fontsize',18); ylabel('N^2 (10^{-4} S^{-2})','FontSize',18);xlabel('Month','fontsize',20); dn = 1:1:365; UniformFit = ftN2Composite.Coefficients.Estimate(1)+ftN2Composite.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ftN2Composite.Coefficients.Estimate(3)))+ftN2Composite.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ftN2Composite.Coefficients.Estimate(5))); plot(1:1:365,UniformFit*10^4,'r','linewidth',2); title('9-year Composite N^2 (30-50 km)'); % legend('N^2', 'Fitting','fontsize',16); %% seven steps mean and 1 point step SevenStepMean = zeros(1,size(mean30to50EpmTotal,2)); for i = 1:size(mean30to50EpmTotal,2) if i < 4 SevenStepMean(1,i) = mean(mean30to50EpmTotal(1:(i+3))); elseif i > size(mean30to50EpmTotal,2)-3 SevenStepMean(1,i) = mean(mean30to50EpmTotal((i-3):end)); else SevenStepMean(1,i) = mean(mean30to50EpmTotal((i-3):(i+3))); end end f = figure(2); f.WindowState = 'maximized'; plot(timesortTotal,SevenStepMean,'-r'); ylim([0,20]);xlim([datenum('20110101','yyyymmdd'),datenum('20191231','yyyymmdd')]); set(gca,'YMinorTick','on');set(gca,'TickDir','out'); xlabel('Month','FontSize',14);ylabel('E_{pm} (J/kg)','FontSize',20); set(gca,'xtick',[datenum('20110101','yyyymmdd'),datenum('20110201','yyyymmdd'),datenum('20110301','yyyymmdd'),... datenum('20110401','yyyymmdd'),datenum('20110501','yyyymmdd'),datenum('20110601','yyyymmdd'),... datenum('20110701','yyyymmdd'),datenum('20110801','yyyymmdd'),datenum('20110901','yyyymmdd'),... datenum('20111001','yyyymmdd'),datenum('20111101','yyyymmdd'),datenum('20111201','yyyymmdd'),... datenum('20120101','yyyymmdd'),datenum('20120201','yyyymmdd'),datenum('20120301','yyyymmdd'),... datenum('20120401','yyyymmdd'),datenum('20120501','yyyymmdd'),datenum('20120601','yyyymmdd'),... datenum('20120701','yyyymmdd'),datenum('20120801','yyyymmdd'),datenum('20120901','yyyymmdd'),... datenum('20121001','yyyymmdd'),datenum('20121101','yyyymmdd'),datenum('20121201','yyyymmdd'),... datenum('20130101','yyyymmdd'),datenum('20130201','yyyymmdd'),datenum('20130301','yyyymmdd'),... datenum('20130401','yyyymmdd'),datenum('20130501','yyyymmdd'),datenum('20130601','yyyymmdd'),... datenum('20130701','yyyymmdd'),datenum('20130801','yyyymmdd'),datenum('20130901','yyyymmdd'),... datenum('20131001','yyyymmdd'),datenum('20131101','yyyymmdd'),datenum('20131201','yyyymmdd'),... datenum('20140101','yyyymmdd'),datenum('20140201','yyyymmdd'),datenum('20140301','yyyymmdd'),... datenum('20140401','yyyymmdd'),datenum('20140501','yyyymmdd'),datenum('20140601','yyyymmdd'),... datenum('20140701','yyyymmdd'),datenum('20140801','yyyymmdd'),datenum('20140901','yyyymmdd'),... datenum('20141001','yyyymmdd'),datenum('20141101','yyyymmdd'),datenum('20141201','yyyymmdd'),... datenum('20150101','yyyymmdd'),datenum('20150201','yyyymmdd'),datenum('20150301','yyyymmdd'),... datenum('20150401','yyyymmdd'),datenum('20150501','yyyymmdd'),datenum('20150601','yyyymmdd'),... datenum('20150701','yyyymmdd'),datenum('20150801','yyyymmdd'),datenum('20150901','yyyymmdd'),... datenum('20151001','yyyymmdd'),datenum('20151101','yyyymmdd'),datenum('20151201','yyyymmdd'),... datenum('20160101','yyyymmdd'),datenum('20160201','yyyymmdd'),datenum('20160301','yyyymmdd'),... datenum('20160401','yyyymmdd'),datenum('20160501','yyyymmdd'),datenum('20160601','yyyymmdd'),... datenum('20160701','yyyymmdd'),datenum('20160801','yyyymmdd'),datenum('20160901','yyyymmdd'),... datenum('20161001','yyyymmdd'),datenum('20161101','yyyymmdd'),datenum('20161201','yyyymmdd'),... datenum('20170101','yyyymmdd'),datenum('20170201','yyyymmdd'),datenum('20170301','yyyymmdd'),... datenum('20170401','yyyymmdd'),datenum('20170501','yyyymmdd'),datenum('20170601','yyyymmdd'),... datenum('20170701','yyyymmdd'),datenum('20170801','yyyymmdd'),datenum('20170901','yyyymmdd'),... datenum('20171001','yyyymmdd'),datenum('20171101','yyyymmdd'),datenum('20171201','yyyymmdd'),... datenum('20180101','yyyymmdd'),datenum('20180201','yyyymmdd'),datenum('20180301','yyyymmdd'),... datenum('20180401','yyyymmdd'),datenum('20180501','yyyymmdd'),datenum('20180601','yyyymmdd'),... datenum('20180701','yyyymmdd'),datenum('20180801','yyyymmdd'),datenum('20180901','yyyymmdd'),... datenum('20181001','yyyymmdd'),datenum('20181101','yyyymmdd'),datenum('20181201','yyyymmdd'),... datenum('20190101','yyyymmdd'),datenum('20190201','yyyymmdd'),datenum('20190301','yyyymmdd'),... datenum('20190401','yyyymmdd'),datenum('20190501','yyyymmdd'),datenum('20190601','yyyymmdd'),... datenum('20190701','yyyymmdd'),datenum('20190801','yyyymmdd'),datenum('20190901','yyyymmdd'),... datenum('20191001','yyyymmdd'),datenum('20191101','yyyymmdd'),datenum('20191201','yyyymmdd'),],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',15); text(datenum('20110101','yyyymmdd'),-1.5,'2011');text(datenum('20120101','yyyymmdd'),-1.5,'2012'); text(datenum('20130101','yyyymmdd'),-1.5,'2013');text(datenum('20140101','yyyymmdd'),-1.5,'2014'); text(datenum('20150101','yyyymmdd'),-1.5,'2015');text(datenum('20160101','yyyymmdd'),-1.5,'2016'); text(datenum('20170101','yyyymmdd'),-1.5,'2017');text(datenum('20180101','yyyymmdd'),-1.5,'2018'); text(datenum('20190101','yyyymmdd'),-1.5,'2019'); ylabel('Epm ( J/kg )'); %% plot distance from polar vortex edge to McMurdo path = 'C:\Users\l1989\Desktop\PolarVortex\FigureofPolarVortex\merra2_at_mcmurdo_20100101-20161231.nc'; PolarVortexDates = ncread(path,'DATES'); PolarVortexAlt = ncread(path, 'ALTITUDE'); PolarVortexDis = ncread(path, 'EDGE_DISTANCE'); PolarVortexDatenum = datenum(num2str(PolarVortexDates),'yyyymmdd'); AltIndex = 9; yyaxis right; xValues = datenum('20110101','yyyymmdd'):1:datenum('20191231','yyyymmdd'); SmthDis = smooth(PolarVortexDis(:,AltIndex),7); SmthDis(isnan(PolarVortexDis(:,AltIndex))) = nan; plot(PolarVortexDatenum,SmthDis);hold on; plot(xValues, zeros(1,length(xValues)),'--','LineWidth',1); ylabel('polar vortex edge distance ( km )'); %% add title title('60 km'); %% each month mean Epm (2011-2019) EpmMonthlyMean = zeros(12,9); EpmMonthlyMean_std = zeros(12,9); EpmMonthlyMean_SE = zeros(12,9); for year = 2011:2019 for month = 1:1:12 ind = find(timesort_datetime.Year == year & timesort_datetime.Month == month); EpmMonthlyMean(month, year-2010) = mean(mean30to50EpmTotal(ind)); EpmMonthlyMean_std(month, year-2010) = std(mean30to50EpmTotal(ind)); EpmMonthlyMean_SE(month, year-2010) = std(mean30to50EpmTotal(ind))/sqrt(length(mean30to50EpmTotal(ind))); EpmMonthlyMean_TimeSort(month, year-2010) = datetime(year,month,15); textTime(month, year-2010) = year*10000+month*100+15; end end WinterMonth567Mean = zeros(9,1); WinterMonth567Mean_std = zeros(9,1); WinterMonth567Mean_SE = zeros(9,1); WinterMonth5678Mean = zeros(9,1); WinterMonth5678Mean_std = zeros(9,1); WinterMonth5678Mean_SE = zeros(9,1); for year = 2011:2019 ind1 = find(timesort_datetime.Year == year & timesort_datetime.Month == 5); ind2 = find(timesort_datetime.Year == year & timesort_datetime.Month == 6); ind3 = find(timesort_datetime.Year == year & timesort_datetime.Month == 7); ind = [ind1;ind2;ind3]; WinterMonth567Mean(year-2010) = mean(mean30to50EpmTotal(ind)); WinterMonth567Mean_std(year-2010) = std(mean30to50EpmTotal(ind)); WinterMonth567Mean_SE(year-2010) = std(mean30to50EpmTotal(ind))/sqrt(length(mean30to50EpmTotal(ind))); end for year = 2011:2019 ind1 = find(timesort_datetime.Year == year & timesort_datetime.Month == 5); ind2 = find(timesort_datetime.Year == year & timesort_datetime.Month == 6); ind3 = find(timesort_datetime.Year == year & timesort_datetime.Month == 7); ind4 = find(timesort_datetime.Year == year & timesort_datetime.Month == 8); ind = [ind1;ind2;ind3;ind4]; WinterMonth5678Mean(year-2010) = mean(mean30to50EpmTotal(ind)); WinterMonth5678Mean_std(year-2010) = std(mean30to50EpmTotal(ind)); WinterMonth5678Mean_SE(year-2010) = std(mean30to50EpmTotal(ind))/sqrt(length(mean30to50EpmTotal(ind))); end MonthlyMeanThenWinterMean567 = zeros(9,1); MonthlyMeanThenWinterMean567_std = zeros(9,1); MonthlyMeanThenWinterMean567_SE = zeros(9,1); for year = 2011:2019 if year == 2018% kick out the 2018 June MonthlyMeanThenWinterMean567(year-2010) = mean(EpmMonthlyMean(5:2:7,year-2010)); MonthlyMeanThenWinterMean567_std(year-2010) = std(EpmMonthlyMean(5:2:7,year-2010)); MonthlyMeanThenWinterMean567_SE(year-2010) = std(EpmMonthlyMean(5:2:7,year-2010))/sqrt(length(EpmMonthlyMean(5:2:7,year-2010))); else MonthlyMeanThenWinterMean567(year-2010) = mean(EpmMonthlyMean(5:7,year-2010)); MonthlyMeanThenWinterMean567_std(year-2010) = std(EpmMonthlyMean(5:7,year-2010)); MonthlyMeanThenWinterMean567_SE(year-2010) = std(EpmMonthlyMean(5:7,year-2010))/sqrt(length(EpmMonthlyMean(5:7,year-2010))); end end MonthlyMeanThenWinterMean5678 = zeros(9,1); MonthlyMeanThenWinterMean5678_std = zeros(9,1); MonthlyMeanThenWinterMean5678_SE = zeros(9,1); for year = 2011:2019 if year == 2018% kick out the 2018 June MonthlyMeanThenWinterMean5678(year-2010) = mean(EpmMonthlyMean([5,7,8],year-2010)); MonthlyMeanThenWinterMean5678_std(year-2010) = std(EpmMonthlyMean([5,7,8],year-2010)); MonthlyMeanThenWinterMean5678_SE(year-2010) = std(EpmMonthlyMean([5,7,8],year-2010))/sqrt(length(EpmMonthlyMean([5,7,8],year-2010))); else MonthlyMeanThenWinterMean5678(year-2010) = mean(EpmMonthlyMean(5:8,year-2010)); MonthlyMeanThenWinterMean5678_std(year-2010) = std(EpmMonthlyMean(5:8,year-2010)); MonthlyMeanThenWinterMean5678_SE(year-2010) = std(EpmMonthlyMean(5:8,year-2010))/sqrt(length(EpmMonthlyMean(5:8,year-2010))); end end WinterMonth567Mean_NoJun2018 = zeros(9,1); WinterMonth567Mean_std_NoJun2018 = zeros(9,1); WinterMonth567Mean_SE_NoJun2018 = zeros(9,1); WinterMonth5678Mean_NoJun2018 = zeros(9,1); WinterMonth5678Mean_std_NoJun2018 = zeros(9,1); WinterMonth5678Mean_SE_NoJun2018 = zeros(9,1); for year = 2011:2019 if year == 2018 ind1 = find(timesort_datetime.Year == year & timesort_datetime.Month == 5); ind3 = find(timesort_datetime.Year == year & timesort_datetime.Month == 7); ind = [ind1;ind3]; else ind1 = find(timesort_datetime.Year == year & timesort_datetime.Month == 5); ind2 = find(timesort_datetime.Year == year & timesort_datetime.Month == 6); ind3 = find(timesort_datetime.Year == year & timesort_datetime.Month == 7); ind = [ind1;ind2;ind3]; end WinterMonth567Mean_NoJun2018(year-2010) = mean(mean30to50EpmTotal(ind)); WinterMonth567Mean_std_NoJun2018(year-2010) = std(mean30to50EpmTotal(ind)); WinterMonth567Mean_SE_NoJun2018(year-2010) = std(mean30to50EpmTotal(ind))/sqrt(length(mean30to50EpmTotal(ind))); end for year = 2011:2019 if year == 2018 ind1 = find(timesort_datetime.Year == year & timesort_datetime.Month == 5); ind3 = find(timesort_datetime.Year == year & timesort_datetime.Month == 7); ind4 = find(timesort_datetime.Year == year & timesort_datetime.Month == 8); ind = [ind1;ind3;ind4]; else ind1 = find(timesort_datetime.Year == year & timesort_datetime.Month == 5); ind2 = find(timesort_datetime.Year == year & timesort_datetime.Month == 6); ind3 = find(timesort_datetime.Year == year & timesort_datetime.Month == 7); ind4 = find(timesort_datetime.Year == year & timesort_datetime.Month == 8); ind = [ind1;ind2;ind3;ind4]; end WinterMonth5678Mean_NoJun2018(year-2010) = mean(mean30to50EpmTotal(ind)); WinterMonth5678Mean_std_NoJun2018(year-2010) = std(mean30to50EpmTotal(ind)); WinterMonth5678Mean_SE_NoJun2018(year-2010) = std(mean30to50EpmTotal(ind))/sqrt(length(mean30to50EpmTotal(ind))); end %% plot Singapore data load('SingarporeData.mat'); SingarporeTimeSeries = datetime(SingarporeData(2:end,2),SingarporeData(2:end,1),15.*ones(length(SingarporeData(2:end,1)),1)); SingarporeDateNum = datenum(SingarporeTimeSeries); f = figure(1); f.WindowState = 'maximized'; subplot(2,1,1); errorbar(timesortTotal-dn_start,mean30to50EpmTotal,ErrorMean30to50EpmTotal,'*b'); ylim([0,20]);xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'YMinorTick','on');set(gca,'TickDir','out'); xlabel('Month','FontSize',14);ylabel('E_{pm} (J/kg)','FontSize',20); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',15); text(datenum('20110101','yyyymmdd')-dn_start,-1.5,'2011');text(datenum('20120101','yyyymmdd')-dn_start,-1.5,'2012'); text(datenum('20130101','yyyymmdd')-dn_start,-1.5,'2013');text(datenum('20140101','yyyymmdd')-dn_start,-1.5,'2014'); text(datenum('20150101','yyyymmdd')-dn_start,-1.5,'2015');text(datenum('20160101','yyyymmdd')-dn_start,-1.5,'2016'); text(datenum('20170101','yyyymmdd')-dn_start,-1.5,'2017');text(datenum('20180101','yyyymmdd')-dn_start,-1.5,'2018'); text(datenum('20190101','yyyymmdd')-dn_start,-1.5,'2019'); for nYr = 2011:1:2019 hold on; index = find(timesort_datetime.Year == nYr); dnAll = (datenum(nYr,1,1)-dn_start):1:(datenum(nYr,12,31)-dn_start); plot(dnAll,ft{nYr-2010}.Coefficients.Estimate(1) + ft{nYr-2010}.Coefficients.Estimate(2).*cos((2*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(3)))+ft{nYr-2010}.Coefficients.Estimate(4).*cos((4*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(5))), 'r','LineWidth',2); end subplot(2,1,2); plot(SingarporeDateNum-dn_start, mean(SingarporeData(2:end,7),2),'m','LineWidth',4); xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'YMinorTick','on');set(gca,'TickDir','out'); xlabel('Month','FontSize',14);ylabel('Wind Speed (m/s)','FontSize',20); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',15); text(datenum('20110101','yyyymmdd')-dn_start,-22,'2011');text(datenum('20120101','yyyymmdd')-dn_start,-22,'2012'); text(datenum('20130101','yyyymmdd')-dn_start,-22,'2013');text(datenum('20140101','yyyymmdd')-dn_start,-22,'2014'); text(datenum('20150101','yyyymmdd')-dn_start,-22,'2015');text(datenum('20160101','yyyymmdd')-dn_start,-22,'2016'); text(datenum('20170101','yyyymmdd')-dn_start,-22,'2017');text(datenum('20180101','yyyymmdd')-dn_start,-22,'2018'); text(datenum('20190101','yyyymmdd')-dn_start,-22,'2019'); % hold on; % plot(SingarporeDateNum,SingarporeData(2:end,5),'LineWidth',2); % plot(SingarporeDateNum,SingarporeData(2:end,7),'LineWidth',2); % plot(SingarporeDateNum,SingarporeData(2:end,9),'LineWidth',2); set(gca,'Clipping','Off');refline([0 0]); line([735062-dn_start 735062-dn_start],[-20,76],'LineStyle','--','color','black','LineWidth',3); line([736150-dn_start 736150-dn_start],[-20,76],'LineStyle','--','color','black','LineWidth',3); line([737255-dn_start 737255-dn_start],[-20,76],'LineStyle','--','color','black','LineWidth',3);ylim([-20 20]); legend('mean of 30,40,50,70,80 hPa');%,'80hPa','50hPa','30hPa','xxx' %% MERRA2 zonal wind and meridional wind load('F:\New folder\MERRA2\ZonalMeanUV(-90to0S).mat'); StartDate = [2011 1]; EndDate = [2019 12]; %% calc monthly mean zonal mean wind - MERRA2 i = 1; T = datetime(StartDate(1),StartDate(2),1):datetime(EndDate(1),EndDate(2),31); for yr = StartDate(1):EndDate(1) for mon = 1:1:12 MonInd{i} = find(T.Year == yr & T.Month == mon); TMon(i) = datetime(yr,mon,1); i = i+1; end end load('F:\New folder\MERRA2\MERRA2_ModelLevel 72.txt'); alt = atmospalt(MERRA2_ModelLevel_72(:,2).*100)./1000; for nCnt = 1:length(ZonalMeanU) ZonalMeanU{nCnt} = ZonalMeanU{nCnt}(end-160:end,:); ZonalMeanV{nCnt} = ZonalMeanV{nCnt}(end-160:end,:); end MERRA2Pressure = 50; SingarporePressure = 50; % altRg = 30:30; latRg = 161:161; [~,altRg] = min(abs(MERRA2_ModelLevel_72(:,2)-MERRA2Pressure)); [~,SingInd] = min(abs(SingarporePressure-SingarporeData(1,1:12))); etrZonalMeanU = zeros(1, length(ZonalMeanU)); etrZonalMeanV = zeros(1, length(ZonalMeanU)); for nCnt = 1:1:length(ZonalMeanU) etrZonalMeanU(nCnt) = mean(ZonalMeanU{nCnt}(latRg,altRg),'all'); stdZonalMeanU(nCnt) = std(ZonalMeanU{nCnt}(latRg,altRg),0,'all'); etrZonalMeanV(nCnt) = mean(ZonalMeanV{nCnt}(latRg,altRg),'all'); stdZonalMeanV(nCnt) = std(ZonalMeanV{nCnt}(latRg,altRg),0,'all'); end f=figure(4); f.WindowState = 'maximized'; plot(1:1:3287,etrZonalMeanU);hold on; xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'YMinorTick','on');set(gca,'TickDir','out'); ylabel('Wind Speed (m/s)','FontSize',20); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',15); text(datenum('20110101','yyyymmdd')-dn_start,-46,'2011');text(datenum('20120101','yyyymmdd')-dn_start,-46,'2012'); text(datenum('20130101','yyyymmdd')-dn_start,-46,'2013');text(datenum('20140101','yyyymmdd')-dn_start,-46,'2014'); text(datenum('20150101','yyyymmdd')-dn_start,-46,'2015');text(datenum('20160101','yyyymmdd')-dn_start,-46,'2016'); text(datenum('20170101','yyyymmdd')-dn_start,-46,'2017');text(datenum('20180101','yyyymmdd')-dn_start,-46,'2018'); text(datenum('20190101','yyyymmdd')-dn_start,-46,'2019'); plot(SingarporeDateNum-dn_start, mean(SingarporeData(2:end,SingInd),2),'r','LineWidth',2); plot(SingarporeDateNum-dn_start,zeros(length(SingarporeDateNum),1),'--'); % plot(SingarporeDateNum-dn_start, mean(SingarporeData(2:end,SingInd+1),2),'--r','LineWidth',2); legend({['MERRA2 ' num2str(MERRA2Pressure) 'hPa'],['Singarpore ' num2str(SingarporePressure) 'hPa']}); %% interp to even altitude range ( 80:-0.2:0 km) interpVerticalZonalMeanU = cell(1,length(ZonalMeanU)); interpVerticalZonalMeanV = cell(1,length(ZonalMeanV)); interpAlt = 80:-0.2:0; for nDay = 1:length(ZonalMeanU) for nLat = 1:size(ZonalMeanU{1},1) xd = interp1(alt,ZonalMeanU{nDay}(nLat,:),interpAlt); interpVerticalZonalMeanU{nDay}(nLat,:) = xd; xd = interp1(alt,ZonalMeanV{nDay}(nLat,:),interpAlt); interpVerticalZonalMeanV{nDay}(nLat,:) = xd; end end for nCont = 1:length(MonInd) SumU = zeros(length(lat),length(lev)); SumV = zeros(length(lat),length(lev)); for nCt = 1:length(MonInd{nCont}) SumU = SumU+ZonalMeanU{MonInd{nCont}(nCt)}(end-160:end,:); SumV = SumV+ZonalMeanV{MonInd{nCont}(nCt)}(end-160:end,:); end MonZonalMeanU{nCont} = SumU./length(MonInd{nCont}); MonZonalMeanV{nCont} = SumV./length(MonInd{nCont}); end for nCont = 1:length(MonInd) SumU = zeros(length(lat),length(interpAlt)); SumV = zeros(length(lat),length(interpAlt)); for nCt = 1:length(MonInd{nCont}) SumU = SumU+interpVerticalZonalMeanU{MonInd{nCont}(nCt)}; SumV = SumV+interpVerticalZonalMeanU{MonInd{nCont}(nCt)}; end interpMonZonalMeanU{nCont} = SumU./length(MonInd{nCont}); interpMonZonalMeanV{nCont} = SumV./length(MonInd{nCont}); end %% different method to quantify filter effect % zonal mean wind altRg = 30:30; latRg = 161:161; etrMonZonalMeanU = zeros(1,length(MonZonalMeanU)); etrMonZonalMeanV = zeros(1,length(MonZonalMeanV)); MonVecSum = zeros(1,length(MonZonalMeanU)); for nCnt = 1:1:length(MonZonalMeanU) etrMonZonalMeanU(nCnt) = mean(MonZonalMeanU{nCnt}(latRg,altRg),'all'); stdMonZonalMeanU(nCnt) = std(MonZonalMeanU{nCnt}(latRg,altRg),0,'all'); etrMonZonalMeanV(nCnt) = mean(MonZonalMeanV{nCnt}(latRg,altRg),'all'); stdMonZonalMeanV(nCnt) = std(MonZonalMeanV{nCnt}(latRg,altRg),0,'all'); MonVecSum(nCnt) = sqrt(etrMonZonalMeanU(nCnt).^2+etrMonZonalMeanV(nCnt).^2); end etrZonalMeanU = zeros(1, length(ZonalMeanU)); etrZonalMeanV = zeros(1, length(ZonalMeanU)); for nCnt = 1:1:length(ZonalMeanU) etrZonalMeanU(nCnt) = mean(ZonalMeanU{nCnt}(latRg,altRg),'all'); stdZonalMeanU(nCnt) = std(ZonalMeanU{nCnt}(latRg,altRg),0,'all'); etrZonalMeanV(nCnt) = mean(ZonalMeanV{nCnt}(latRg,altRg),'all'); stdZonalMeanV(nCnt) = std(ZonalMeanV{nCnt}(latRg,altRg),0,'all'); end % interped interpAltRg = 151:376; etrInterpMonZonalMeanU = zeros(1,length(MonZonalMeanU)); etrInterpMonZonalMeanV = zeros(1,length(MonZonalMeanV)); for nCnt = 1:1:length(MonZonalMeanU) etrInterpMonZonalMeanU(nCnt) = mean(interpMonZonalMeanU{nCnt}(latRg,interpAltRg),'all'); stdInterpMonZonalMeanU(nCnt) = std(interpMonZonalMeanU{nCnt}(latRg,interpAltRg),0,'all'); etrInterpMonZonalMeanV(nCnt) = mean(interpMonZonalMeanV{nCnt}(latRg,interpAltRg),'all'); stdInterpMonZonalMeanV(nCnt) = std(interpMonZonalMeanV{nCnt}(latRg,interpAltRg),0,'all'); end interpZonalMeanU = zeros(1, length(ZonalMeanU)); interpZonalMeanV = zeros(1, length(ZonalMeanU)); for nCnt = 1:1:length(ZonalMeanU) etrInterpZonalMeanU(nCnt) = nanmean(interpVerticalZonalMeanU{nCnt}(latRg,interpAltRg),'all'); stdInterpZonalMeanU(nCnt) = std(interpVerticalZonalMeanU{nCnt}(latRg,interpAltRg),0,'all'); etrInterpZonalMeanV(nCnt) = nanmean(interpVerticalZonalMeanV{nCnt}(latRg,interpAltRg),'all'); stdInterpZonalMeanV(nCnt) = std(interpVerticalZonalMeanV{nCnt}(latRg,interpAltRg),0,'all'); end QuantifyMethod = 1; switch QuantifyMethod case 1 %monthly mean zonal wind YData = etrMonZonalMeanU; case 2 %monthly mean meridional wind YData = etrMonZonalMeanV; case 3 YData = MonVecSum; case 4 YData = stdZonalMeanU; case 5 YData = stdMonZonalMeanU; case 6 YData = smooth(stdZonalMeanU,7,'rlowess'); case 7 YData = etrInterpZonalMeanU; case 8 YData = stdInterpZonalMeanU; case 9 YData = smooth(stdInterpZonalMeanU,7,'rlowess'); case 10 YData = etrInterpMonZonalMeanU; case 11 YData = stdInterpMonZonalMeanU; end dn_start = datenum(StartDate(1),StartDate(2),1)-1; hold on; yyaxis right; switch QuantifyMethod case 4 p3 = plot(datenum(T)-dn_start,YData); case 6 p3 = plot(datenum(T)-dn_start,YData); case 7 p3 = plot(datenum(T)-dn_start,YData); case 8 p3 = plot(datenum(T)-dn_start,YData); case 9 p3 = plot(datenum(T)-dn_start,YData); otherwise p3 = plot(datenum(TMon)-dn_start+15,YData); end % mark Jun Jul Aug i=1; for yr = StartDate(1):EndDate(1) for mon = 1:1:12 switch mon case 6 hold on; JunM = YData(i); scatter(datenum(TMon(i))-dn_start+15,YData(i),80,'filled'); case 7 hold on; JulM = YData(i); scatter(datenum(TMon(i))-dn_start+15,YData(i),80,'filled'); case 8 hold on; AugM = YData(i); % scatter(datenum(TMon(i))-dn_start+15,YData(i),80,'filled'); case 9 hold on; % scatter(datenum(TMon(i))-dn_start+15,YData(i),80,'filled'); SepM = YData(i); JJASM = (JunM+JulM+AugM)/3; % scatter(datenum(TMon(i))-dn_start-4,JJASM,200,'filled','d'); end i=i+1; end end %% add title and legend title(['Altitude ' num2str(interpAlt(interpAltRg(end))) ' to ' num2str(interpAlt(interpAltRg((1)))) ' km, latitude -80 to -60 S']); title(['Altitude ' num2str(alt(altRg(end))) ' to ' num2str(alt(altRg(1))) ' km, latitude -80 to -60 S']); legend([p1 p2 p3],{'Lidar Epm','Fitting Epm','std 7day Smooth U'}); %% Epm anomaly - figure 4(a)(b) clear absEpmAnomaly; clear absMonEpmAnomaly; clear relEpmAnomaly; clear relMonEpmAnomaly; clear SingarporeWinData; clear uniformfit; EpmMonthlyMean = zeros(12,9); for year = 2011:2019 for month = 1:1:12 ind = find(timesort_datetime.Year == year & timesort_datetime.Month == month); EpmMonthlyMean(month, year-2010) = mean(mean30to50EpmTotal(ind)); EpmMonthlyMean_TimeSort(month, year-2010) = datetime(year,month,15); end end method =1; % method 1: using annual fitting curve as minuend % method 2: using monthly mean as minuend uniformfitMethod = 2;% 1 - 9 years uniformfit, 2 - composite fit % take mean for 9 years Epm in the same day of a year n = 1; CompositeEpm = cell(365,1); for nMon = 1:1:12 for nDay = 1:1:eomday(2011,nMon) Ind = timesort_datetime.Month == nMon & timesort_datetime.Day == nDay; CompositeEpm{n} = mean30to50EpmTotal(Ind); n = n+1; end end % take mean for every 5 days n = 1; binDay = 5; yData = zeros(365/binDay,1); xData = zeros(365/binDay,1); for nCnt = 1:binDay:365 yData(n) = mean([CompositeEpm{1+(n-1)*binDay} CompositeEpm{2+(n-1)*binDay} CompositeEpm{3+(n-1)*binDay} CompositeEpm{4+(n-1)*binDay} CompositeEpm{5+(n-1)*binDay}]); xData(n) = nCnt+2; n = n+1; end % fitting ftComposite = fitnlm(xData,yData','y ~ b1 + b2.*cos((2*pi/365).*(x1-b3))+b4.*cos((4*pi/365).*(x1-b5))',beta0); % calc Epm anomaly figure(1); switch method case 1 for nYr = 2011:1:2019 dnAll = (datenum(nYr,1,1)-dn_start):1:(datenum(nYr,12,31)-dn_start); AnnualFit = ft{nYr-2010}.Coefficients.Estimate(1) + ft{nYr-2010}.Coefficients.Estimate(2).*cos((2*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(3)))+ft{nYr-2010}.Coefficients.Estimate(4).*cos((4*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(5))); if uniformfitMethod == 1 Ind = dtSeries9yrs.Year == nYr; dn = TimeSeries9yrs(Ind); UniformFit = ft9yrs.Coefficients.Estimate(1)+ft9yrs.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ft9yrs.Coefficients.Estimate(3)))+ft9yrs.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ft9yrs.Coefficients.Estimate(5))); elseif uniformfitMethod == 2 Ind = dtSeries9yrs.Year == nYr; dn = TimeSeries9yrs(Ind); UniformFit = ftComposite.Coefficients.Estimate(1)+ftComposite.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ftComposite.Coefficients.Estimate(3)))+ftComposite.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ftComposite.Coefficients.Estimate(5))); end plot(dnAll, AnnualFit);hold on; plot(dn, UniformFit,'r');hold on; absEpmAnomaly{nYr-2010} = AnnualFit - UniformFit; % plot(dn,absEpmAnomaly{nYr-2010},'b');hold on; relEpmAnomaly{nYr-2010} = absEpmAnomaly{nYr-2010}./UniformFit; % plot(dn,relEpmAnomaly{nYr-2010},'g');hold on; end n = 1; for nYr = 2011:1:2019 for nMon = 1:1:12 dtEachyear = datetime(nYr,1,1):1:datetime(nYr,12,31); Ind = find(dtEachyear.Year == nYr & dtEachyear.Month == nMon); absMonEpmAnomaly(n) = mean(absEpmAnomaly{nYr-2010}(Ind)); relMonEpmAnomaly(n) = mean(relEpmAnomaly{nYr-2010}(Ind)); dtMonTimeseries9yrs(n) = datetime(nYr,nMon,15); n = n+1; end end dnMonTimeseries9yrs = datenum(dtMonTimeseries9yrs)-dn_start; % scatter(dnMonTimeseries9yrs,absMonEpmAnomaly,'*');hold on; % scatter(dnMonTimeseries9yrs,relMonEpmAnomaly,'O'); case 2 n = 1; for nYr = 2011:1:2019 if uniformfitMethod == 1 Ind = dtSeries9yrs.Year == nYr; dn = TimeSeries9yrs(Ind); UniformFit = ft9yrs.Coefficients.Estimate(1)+ft9yrs.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ft9yrs.Coefficients.Estimate(3)))+ft9yrs.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ft9yrs.Coefficients.Estimate(5))); elseif uniformfitMethod == 2 Ind = dtSeries9yrs.Year == nYr; dn = TimeSeries9yrs(Ind); UniformFit = ftComposite.Coefficients.Estimate(1)+ftComposite.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ftComposite.Coefficients.Estimate(3)))+ftComposite.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ftComposite.Coefficients.Estimate(5))); end for nMon = 1:1:12 dtEachyear = datetime(nYr,1,1):1:datetime(nYr,12,31); IndYM = find(dtEachyear.Year == nYr & dtEachyear.Month == nMon); UniformMonFit(n) = mean(UniformFit(IndYM)); absMonEpmAnomaly(n) = EpmMonthlyMean(nMon,nYr-2010) - UniformMonFit(n); relMonEpmAnomaly(n) = absMonEpmAnomaly(n)/UniformMonFit(n); scatter(datenum(nYr,nMon,15)-dn_start,absMonEpmAnomaly(n),'*b');hold on; scatter(datenum(nYr,nMon,15)-dn_start,relMonEpmAnomaly(n),'o');hold on; n = n+1; end end dnMonTimeseries9yrs = datenum(dtMonTimeseries9yrs)-dn_start; plot(dnMonTimeseries9yrs, absMonEpmAnomaly);hold on; plot(dnMonTimeseries9yrs, relMonEpmAnomaly);hold on; plot(dnMonTimeseries9yrs, UniformMonFit); ylim([-5 5]); end plot(dnMonTimeseries9yrs, zeros(1,size(dnMonTimeseries9yrs,2)),'--'); ylim([-15 15]); % get MERRA2 QBO Index load('C:\Users\l1989\Desktop\PolarVortex\FigureofPolarVortex\MERRA2QBOIndex.mat'); PressureVal = [40 50]; %unit: hPa [PressureInd,~] = find(lev == PressureVal); [StartDateInd,~] = find(dtTimeSeriesQBO == datetime(2011,1,1)); [EndDateInd,~] = find(dtTimeSeriesQBO == datetime(2019,11,31)); QBOIndex = mean(ZonalMeanWind(StartDateInd:EndDateInd,PressureInd),2); % QBOIndex = WindShear; % extract May to Aug n=1; Ind = 0; for nYr = 2011:1:2019 for nMon = 1:1:12 if nMon >=6 && nMon <=7 && nYr ~= 2016 absWinMonEpmAnomaly(n) = absMonEpmAnomaly(nMon+(nYr-2011)*12); relWinMonEpmAnomaly(n) = relMonEpmAnomaly(nMon+(nYr-2011)*12); % Ind = find(SingarporeTimeSeries.Year == nYr & SingarporeTimeSeries.Month == nMon); % SingarporeWinData(n) = mean(SingarporeData(Ind+1,7:8),2); xxData(n) = datetime(nYr,nMon,15); else absWinMonEpmAnomaly(n) = NaN; relWinMonEpmAnomaly(n) = NaN; xxData(n) = datetime(nYr,nMon,15); % SingarporeWinData(n) = NaN; end n = n+1; end end % set 2018 June to Nan absWinMonEpmAnomaly(90) = NaN; relWinMonEpmAnomaly(90) = NaN; % MMon = [6 7]; % for nYr = 2011:1:2019 % IndMM = find(xxData.Year == nYr & xxData.Month >= MMon(1) & xxData.Month <= MMon(end)); % absXMonEpmAnmaly(nYr-2011,:) = % end xQBOdata = dtTimeSeriesQBO(StartDateInd:EndDateInd); scrsz = get(0, 'ScreenSize'); fig = figure('position',[1,300,scrsz(3)*0.95,scrsz(4)*0.6]); subplot('position',[0.05,0.51,0.89,0.37]); % left_color = [0 0 0]; % right_color = [0 0 0]; % set(fig,'defaultAxesColorOrder',[left_color; right_color]); yyaxis left p1 = scatter(datenum(xxData)-dn_start, absWinMonEpmAnomaly,100,'d','filled');hold on; plot(ones(21,1)*(datenum('20120101','yyyymmdd')-dn_start),-10:10,'-','color',[0.8 0.8 0.8]);hold on; plot(ones(21,1)*(datenum('20130101','yyyymmdd')-dn_start),-10:10,'-','color',[0.8 0.8 0.8]);hold on; plot(ones(21,1)*(datenum('20140101','yyyymmdd')-dn_start),-10:10,'-','color',[0.8 0.8 0.8]);hold on; plot(ones(21,1)*(datenum('20150101','yyyymmdd')-dn_start),-10:10,'-','color',[0.8 0.8 0.8]);hold on; plot(ones(21,1)*(datenum('20160101','yyyymmdd')-dn_start),-10:10,'-','color',[0.8 0.8 0.8]);hold on; plot(ones(21,1)*(datenum('20170101','yyyymmdd')-dn_start),-10:10,'-','color',[0.8 0.8 0.8]);hold on; plot(ones(21,1)*(datenum('20180101','yyyymmdd')-dn_start),-10:10,'-','color',[0.8 0.8 0.8]);hold on; plot(ones(21,1)*(datenum('20190101','yyyymmdd')-dn_start),-10:10,'-','color',[0.8 0.8 0.8]);hold on; scatter(dnMonTimeseries9yrs(66:67),absMonEpmAnomaly(66:67),100,'d','filled','markerfacecolor','k');hold on; scatter(dnMonTimeseries9yrs(90),absMonEpmAnomaly(90),100,'d','filled','markerfacecolor','k'); % plot(1:1:108, zeros(1,length(1:1:108)),'--','LineWidth',1); set(gca,'fontsize',18); ylim([-4 4]);ylabel('E_{pm} Anomaly (J/kg)');set(gca,'box','on');set(gca,'YMinorTick','on');set(gca,'yticklabel',{'';'-2';'0';'2';'4'},'fontsize',18); text(datenum('20110101','yyyymmdd')-dn_start+110,-4.8,'2011','fontsize',16);text(datenum('20120101','yyyymmdd')-dn_start+110,-4.8,'2012','fontsize',16); text(datenum('20130101','yyyymmdd')-dn_start+110,-4.8,'2013','fontsize',16);text(datenum('20140101','yyyymmdd')-dn_start+110,-4.8,'2014','fontsize',16); text(datenum('20150101','yyyymmdd')-dn_start+110,-4.8,'2015','fontsize',16);text(datenum('20160101','yyyymmdd')-dn_start+110,-4.8,'2016','fontsize',16); text(datenum('20170101','yyyymmdd')-dn_start+110,-4.8,'2017','fontsize',16);text(datenum('20180101','yyyymmdd')-dn_start+110,-4.8,'2018','fontsize',16); text(datenum('20190101','yyyymmdd')-dn_start+110,-4.8,'2019','fontsize',16); yyaxis right smthQBOIndex = smooth(QBOIndex,7,'rlowess'); p2 = plot(datenum(dtTimeSeriesQBO(StartDateInd):days(1):dtTimeSeriesQBO(EndDateInd))-dn_start, smthQBOIndex,'linewidth',3);hold on; plot(datenum(xQBOdata)-dn_start,zeros(length(xQBOdata),1),'--','linewidth',1.5);hold on; [IndJJ1,~] = find(xQBOdata.Month == 6 & xQBOdata.Day == 15); [IndJJ2,~] = find(xQBOdata.Month == 7 & xQBOdata.Day == 15); p3 = scatter(datenum(xQBOdata(IndJJ1))-dn_start,smthQBOIndex(IndJJ1),100,'linewidth',2);hold on; scatter(datenum(xQBOdata(IndJJ2))-dn_start,smthQBOIndex(IndJJ2),100,'linewidth',2);hold on; % hold on;plot(1:1:108, SingarporeData(2:end,7)); set(gca,'fontsize',18); ylabel('Zonal Mean Zonal Wind (m/s)');set(gca,'YMinorTick','on'); ylim([-30 30]);xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',15); legend([p1 p2 p3],{'June-July Epm anomaly',['at ' num2str(PressureVal) ' hPa'],'June-July 15th'}); title(['(a)June-July Epm Anomaly versus QBO Index (' num2str(PressureVal) ' hPa)']); subplot('position',[0.05,0.11,0.89,0.4]); p1 = errorbar(timesortTotal-dn_start,mean30to50EpmTotal,ErrorMean30to50EpmTotal,'*b');hold on; plot(ones(21,1)*(datenum('20120101','yyyymmdd')-dn_start),0:20,'-','color',[0.8 0.8 0.8]);hold on; plot(ones(21,1)*(datenum('20130101','yyyymmdd')-dn_start),0:20,'-','color',[0.8 0.8 0.8]);hold on; plot(ones(21,1)*(datenum('20140101','yyyymmdd')-dn_start),0:20,'-','color',[0.8 0.8 0.8]);hold on; plot(ones(21,1)*(datenum('20150101','yyyymmdd')-dn_start),0:20,'-','color',[0.8 0.8 0.8]);hold on; plot(ones(21,1)*(datenum('20160101','yyyymmdd')-dn_start),0:20,'-','color',[0.8 0.8 0.8]);hold on; plot(ones(21,1)*(datenum('20170101','yyyymmdd')-dn_start),0:20,'-','color',[0.8 0.8 0.8]);hold on; plot(ones(21,1)*(datenum('20180101','yyyymmdd')-dn_start),0:20,'-','color',[0.8 0.8 0.8]);hold on; plot(ones(21,1)*(datenum('20190101','yyyymmdd')-dn_start),0:20,'-','color',[0.8 0.8 0.8]);hold on; ylim([0,20]);xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'YMinorTick','on');set(gca,'TickDir','out');set(gca,'TickLength',[0.006 0.0025]); ylabel('E_{pm} (J/kg)','FontSize',20); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',15); text(datenum('20110101','yyyymmdd')-dn_start+110,-4,'2011','fontsize',18);text(datenum('20120101','yyyymmdd')-dn_start+110,-4,'2012','fontsize',18); text(datenum('20130101','yyyymmdd')-dn_start+110,-4,'2013','fontsize',18);text(datenum('20140101','yyyymmdd')-dn_start+110,-4,'2014','fontsize',18); text(datenum('20150101','yyyymmdd')-dn_start+110,-4,'2015','fontsize',18);text(datenum('20160101','yyyymmdd')-dn_start+110,-4,'2016','fontsize',18); text(datenum('20170101','yyyymmdd')-dn_start+110,-4,'2017','fontsize',18);text(datenum('20180101','yyyymmdd')-dn_start+110,-4,'2018','fontsize',18); text(datenum('20190101','yyyymmdd')-dn_start+110,-4,'2019','fontsize',18); for nYr = 2011:1:2019 hold on; dnAll = (datenum(nYr,1,1)-dn_start):1:(datenum(nYr,12,31)-dn_start); p2 = plot(dnAll,ft{nYr-2010}.Coefficients.Estimate(1) + ft{nYr-2010}.Coefficients.Estimate(2).*cos((2*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(3)))+ft{nYr-2010}.Coefficients.Estimate(4).*cos((4*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(5))), 'r','linewidth',2); end % corr QBO index and Epm % take QBO index mean % [IndJ1,~] = find(xQBOdata.Month == 6 | xQBOdata.Day >= 15); % [IndJ2,~] = find(xQBOdata.Month == 7 | xQBOdata.Day <= 15); % QBOIndex1 = QBOIndex(IndJ); % QBOIndexMean = []; % for nYr = 2011:1:2019 % IndTmp = find(xQBOdata(IndJ).Year == nYr); % Ind6 = IndTmp(1:30); % Ind7 = IndTmp(31:61); % if nYr == 2016 % QBOIndexMean(nYr-2010,1) = NaN; % QBOIndexMean(nYr-2010,2) = NaN; % else % QBOIndexMean(nYr-2010,1) = mean(QBOIndex1(Ind6),'all'); % QBOIndexMean(nYr-2010,2) = mean(QBOIndex1(Ind7),'all'); % end % end % QBOIndexMean = reshape(QBOIndexMean',[1,size(QBOIndexMean,1)*size(QBOIndexMean,2)]); % take single date [IndJ1,~] = find(xQBOdata.Month == 6); [IndJ2,~] = find(xQBOdata.Month == 7); QBOIndex1 = QBOIndex(IndJ1); QBOIndex2 = QBOIndex(IndJ2); QBOIndexMean = []; for nYr = 2011:1:2019 IndTmp1 = find(xQBOdata(IndJ1).Year == nYr); IndTmp2 = find(xQBOdata(IndJ2).Year == nYr); if nYr == 2016 QBOIndexMean(nYr-2010,1) = NaN; QBOIndexMean(nYr-2010,2) = NaN; elseif nYr == 2018 QBOIndexMean(nYr-2010,1) = NaN; QBOIndexMean(nYr-2010,2) = mean(QBOIndex2(IndTmp2),'all'); else QBOIndexMean(nYr-2010,1) = mean(QBOIndex1(IndTmp1),'all'); QBOIndexMean(nYr-2010,2) = mean(QBOIndex2(IndTmp2),'all'); end end QBOIndexMean = reshape(QBOIndexMean',[1,size(QBOIndexMean,1)*size(QBOIndexMean,2)]); IndNaN = ~isnan(absWinMonEpmAnomaly); Ind2NaN = ~isnan(QBOIndexMean); [Rabs Pabs]=corrcoef(absWinMonEpmAnomaly(IndNaN),QBOIndexMean(Ind2NaN)); [Rrel Prel]=corrcoef(relWinMonEpmAnomaly(IndNaN),QBOIndexMean(Ind2NaN)); p = polyfit(QBOIndexMean(Ind2NaN),absWinMonEpmAnomaly(IndNaN),1); myFit = NonLinearModel.fit(QBOIndexMean(Ind2NaN),absWinMonEpmAnomaly(IndNaN), 'y ~ b0 + b1*x', [0,1]); f = polyval(p,QBOIndexMean(Ind2NaN)); scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3)*0.39,scrsz(4)*0.55]); subplot('position',[0.12,0.16,0.78,0.7]); plot(QBOIndexMean(Ind2NaN),absWinMonEpmAnomaly(IndNaN),'LineStyle','none','Marker','o','MarkerSize',8,'MarkerFaceColor','b');hold on; plot(QBOIndexMean(Ind2NaN),f,'-r','LineWidth',2); xlim([-26 16]);set(gca,'fontsize',18); set(gca,'box','on');xlabel('QBO Index (m/s)');ylabel('Epm Anomaly (J/kg)'); title({'June-July Epm Anomaly vs. QBO Index (40-50 hPa)'; ['Correlation Coeff = ' num2str(Rabs(1,2)) ' P = ' num2str(1-Pabs(1,2))] }); %% plot Epm versus QBO index (40 50hPa) scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3)*0.95,scrsz(4)*0.7]); subplot('position',[0.055,0.12,0.88,0.8]); yyaxis left p1 = errorbar(timesortTotal-dn_start,mean30to50EpmTotal,ErrorMean30to50EpmTotal,'*b');hold on; for nYr = 2011:1:2019 hold on; dnAll = (datenum(nYr,1,1)-dn_start):1:(datenum(nYr,12,31)-dn_start); p2 = plot(dnAll,ft{nYr-2010}.Coefficients.Estimate(1) + ft{nYr-2010}.Coefficients.Estimate(2).*cos((2*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(3)))+ft{nYr-2010}.Coefficients.Estimate(4).*cos((4*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(5))), '-r','LineWidth',2); end xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); ylim([0 20]); set(gca,'TickDir','out');set(gca,'Ycolor','b'); ylabel('Epm (J/kg)','fontsize',18); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',15); text(datenum('20110101','yyyymmdd')-dn_start+110,-2,'2011','fontsize',18);text(datenum('20120101','yyyymmdd')-dn_start+110,-2,'2012','fontsize',18); text(datenum('20130101','yyyymmdd')-dn_start+110,-2,'2013','fontsize',18);text(datenum('20140101','yyyymmdd')-dn_start+110,-2,'2014','fontsize',18); text(datenum('20150101','yyyymmdd')-dn_start+110,-2,'2015','fontsize',18);text(datenum('20160101','yyyymmdd')-dn_start+110,-2,'2016','fontsize',18); text(datenum('20170101','yyyymmdd')-dn_start+110,-2,'2017','fontsize',18);text(datenum('20180101','yyyymmdd')-dn_start+110,-2,'2018','fontsize',18); text(datenum('20190101','yyyymmdd')-dn_start+110,-2,'2019','fontsize',18); yyaxis right p3 = plot(datenum(dtTimeSeriesQBO(StartDateInd):days(1):dtTimeSeriesQBO(EndDateInd))-dn_start, smooth(QBOIndex,7,'rlowess'),'LineWidth',2,'color',[0 0.5 0]);hold on; plot(datenum(dtTimeSeriesQBO(StartDateInd):days(1):dtTimeSeriesQBO(EndDateInd))-dn_start,zeros(1,length(dtTimeSeriesQBO(StartDateInd):days(1):dtTimeSeriesQBO(EndDateInd))),'--'); ylabel('Zonal Mean Zonal Wind ( m/s )','fontsize',18); ylim([-70 30]);set(gca,'Ycolor',[0 0.5 0]); legend([p1 p2 p3],{'Epm';'Annual fit';'QBO Index at 40 50 hPa'},'fontsize',18); legend('boxoff'); set(gca,'ticklength',[0.008 0.0025]); title('Epm versus QBO Index','position',[(datenum('20150701','yyyymmdd')-dn_start),20.5],'fontsize',18); %% plot Epm(7 points smooth) versus QBO index (wind shear[50-70hPa]) PressureLevel = [50 70]; [PresInd,~] = find(lev == PressureLevel); WindShear = ZonalMeanWind(StartDateInd:EndDateInd,PresInd(1))-ZonalMeanWind(StartDateInd:EndDateInd,PresInd(2)); scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3)*0.95,scrsz(4)*0.6]); subplot('position',[0.05,0.11,0.89,0.82]); yyaxis left hold on; xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); ylabel('Epm (J/kg)'); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',15); text(datenum('20110101','yyyymmdd')-dn_start+110,-1,'2011','fontsize',16);text(datenum('20120101','yyyymmdd')-dn_start+110,-1,'2012','fontsize',16); text(datenum('20130101','yyyymmdd')-dn_start+110,-1,'2013','fontsize',16);text(datenum('20140101','yyyymmdd')-dn_start+110,-1,'2014','fontsize',16); text(datenum('20150101','yyyymmdd')-dn_start+110,-1,'2015','fontsize',16);text(datenum('20160101','yyyymmdd')-dn_start+110,-1,'2016','fontsize',16); text(datenum('20170101','yyyymmdd')-dn_start+110,-1,'2017','fontsize',16);text(datenum('20180101','yyyymmdd')-dn_start+110,-1,'2018','fontsize',16); text(datenum('20190101','yyyymmdd')-dn_start+110,-1,'2019','fontsize',16); yyaxis right p3 = plot(datenum(dtTimeSeriesQBO(StartDateInd):days(1):dtTimeSeriesQBO(EndDateInd))-dn_start, WindShear);hold on; plot(datenum(dtTimeSeriesQBO(StartDateInd):days(1):dtTimeSeriesQBO(EndDateInd))-dn_start,zeros(1,length(dtTimeSeriesQBO(StartDateInd):days(1):dtTimeSeriesQBO(EndDateInd))),'--'); ylabel('Zonal Mean Zonal Wind shear ( m/s per 20hPa )'); ylim([-20 20]); legend([p1 p2 p3],{'7 points Smooth Epm';'June and July';'Wind Shear between 50-70 hPa'}); title('Epm versus QBO Index'); %% %% Epm anomaly versus QBO Index (wind shear 50-70hPa) clear absEpmAnomaly; clear absMonEpmAnomaly; clear relEpmAnomaly; clear relMonEpmAnomaly; clear SingarporeWinData; clear uniformfit; EpmMonthlyMean = zeros(12,9); for year = 2011:2019 for month = 1:1:12 ind = find(timesort_datetime.Year == year & timesort_datetime.Month == month); EpmMonthlyMean(month, year-2010) = mean(mean30to50EpmTotal(ind)); EpmMonthlyMean_TimeSort(month, year-2010) = datetime(year,month,15); end end method =1; % method 1: using annual fitting curve as minuend % method 2: using monthly mean as minuend uniformfitMethod = 2;% 1 - 9 years uniformfit, 2 - composite fit % take mean for 9 years Epm in the same day of a year n = 1; CompositeEpm = cell(365,1); for nMon = 1:1:12 for nDay = 1:1:eomday(2011,nMon) Ind = timesort_datetime.Month == nMon & timesort_datetime.Day == nDay; CompositeEpm{n} = mean30to50EpmTotal(Ind); n = n+1; end end % take mean for every 5 days n = 1; binDay = 5; yData = zeros(365/binDay,1); xData = zeros(365/binDay,1); for nCnt = 1:binDay:365 yData(n) = mean([CompositeEpm{1+(n-1)*binDay} CompositeEpm{2+(n-1)*binDay} CompositeEpm{3+(n-1)*binDay} CompositeEpm{4+(n-1)*binDay} CompositeEpm{5+(n-1)*binDay}]); xData(n) = nCnt+2; n = n+1; end % fitting ftComposite = fitnlm(xData,yData','y ~ b1 + b2.*cos((2*pi/365).*(x1-b3))+b4.*cos((4*pi/365).*(x1-b5))',beta0); % calc Epm anomaly figure(1); switch method case 1 for nYr = 2011:1:2019 dnAll = (datenum(nYr,1,1)-dn_start):1:(datenum(nYr,12,31)-dn_start); AnnualFit = ft{nYr-2010}.Coefficients.Estimate(1) + ft{nYr-2010}.Coefficients.Estimate(2).*cos((2*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(3)))+ft{nYr-2010}.Coefficients.Estimate(4).*cos((4*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(5))); if uniformfitMethod == 1 Ind = dtSeries9yrs.Year == nYr; dn = TimeSeries9yrs(Ind); UniformFit = ft9yrs.Coefficients.Estimate(1)+ft9yrs.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ft9yrs.Coefficients.Estimate(3)))+ft9yrs.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ft9yrs.Coefficients.Estimate(5))); elseif uniformfitMethod == 2 Ind = dtSeries9yrs.Year == nYr; dn = TimeSeries9yrs(Ind); UniformFit = ftComposite.Coefficients.Estimate(1)+ftComposite.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ftComposite.Coefficients.Estimate(3)))+ftComposite.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ftComposite.Coefficients.Estimate(5))); end plot(dnAll, AnnualFit);hold on; plot(dn, UniformFit,'r');hold on; absEpmAnomaly{nYr-2010} = AnnualFit - UniformFit; % plot(dn,absEpmAnomaly{nYr-2010},'b');hold on; relEpmAnomaly{nYr-2010} = absEpmAnomaly{nYr-2010}./UniformFit; % plot(dn,relEpmAnomaly{nYr-2010},'g');hold on; end n = 1; for nYr = 2011:1:2019 for nMon = 1:1:12 dtEachyear = datetime(nYr,1,1):1:datetime(nYr,12,31); Ind = find(dtEachyear.Year == nYr & dtEachyear.Month == nMon); absMonEpmAnomaly(n) = mean(absEpmAnomaly{nYr-2010}(Ind)); relMonEpmAnomaly(n) = mean(relEpmAnomaly{nYr-2010}(Ind)); dtMonTimeseries9yrs(n) = datetime(nYr,nMon,15); n = n+1; end end dnMonTimeseries9yrs = datenum(dtMonTimeseries9yrs)-dn_start; % scatter(dnMonTimeseries9yrs,absMonEpmAnomaly,'*');hold on; % scatter(dnMonTimeseries9yrs,relMonEpmAnomaly,'O'); case 2 n = 1; for nYr = 2011:1:2019 if uniformfitMethod == 1 Ind = dtSeries9yrs.Year == nYr; dn = TimeSeries9yrs(Ind); UniformFit = ft9yrs.Coefficients.Estimate(1)+ft9yrs.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ft9yrs.Coefficients.Estimate(3)))+ft9yrs.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ft9yrs.Coefficients.Estimate(5))); elseif uniformfitMethod == 2 Ind = dtSeries9yrs.Year == nYr; dn = TimeSeries9yrs(Ind); UniformFit = ftComposite.Coefficients.Estimate(1)+ftComposite.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ftComposite.Coefficients.Estimate(3)))+ftComposite.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ftComposite.Coefficients.Estimate(5))); end for nMon = 1:1:12 dtEachyear = datetime(nYr,1,1):1:datetime(nYr,12,31); IndYM = find(dtEachyear.Year == nYr & dtEachyear.Month == nMon); UniformMonFit(n) = mean(UniformFit(IndYM)); absMonEpmAnomaly(n) = EpmMonthlyMean(nMon,nYr-2010) - UniformMonFit(n); relMonEpmAnomaly(n) = absMonEpmAnomaly(n)/UniformMonFit(n); scatter(datenum(nYr,nMon,15)-dn_start,absMonEpmAnomaly(n),'*b');hold on; scatter(datenum(nYr,nMon,15)-dn_start,relMonEpmAnomaly(n),'o');hold on; n = n+1; end end dnMonTimeseries9yrs = datenum(dtMonTimeseries9yrs)-dn_start; plot(dnMonTimeseries9yrs, absMonEpmAnomaly);hold on; plot(dnMonTimeseries9yrs, relMonEpmAnomaly);hold on; plot(dnMonTimeseries9yrs, UniformMonFit); ylim([-5 5]); end plot(dnMonTimeseries9yrs, zeros(1,size(dnMonTimeseries9yrs,2)),'--'); ylim([-15 15]); % get MERRA2 QBO Index load('C:\Users\l1989\Desktop\PolarVortex\FigureofPolarVortex\MERRA2QBOIndex.mat'); PressureVal = [40 50]; %unit: hPa [PressureInd,~] = find(lev == PressureVal); [StartDateInd,~] = find(dtTimeSeriesQBO == datetime(2011,1,1)); [EndDateInd,~] = find(dtTimeSeriesQBO == datetime(2019,11,31)); QBOIndex = WindShear; % extract May to Aug n=1; Ind = 0; for nYr = 2011:1:2019 for nMon = 1:1:12 if nMon >=6 && nMon <=7 && nYr ~= 2016 absWinMonEpmAnomaly(n) = absMonEpmAnomaly(nMon+(nYr-2011)*12); relWinMonEpmAnomaly(n) = relMonEpmAnomaly(nMon+(nYr-2011)*12); % Ind = find(SingarporeTimeSeries.Year == nYr & SingarporeTimeSeries.Month == nMon); % SingarporeWinData(n) = mean(SingarporeData(Ind+1,7:8),2); xxData(n) = datetime(nYr,nMon,15); else absWinMonEpmAnomaly(n) = NaN; relWinMonEpmAnomaly(n) = NaN; xxData(n) = datetime(nYr,nMon,15); % SingarporeWinData(n) = NaN; end n = n+1; end end % set 2018 June to Nan absWinMonEpmAnomaly(90) = NaN; relWinMonEpmAnomaly(90) = NaN; % MMon = [6 7]; % for nYr = 2011:1:2019 % IndMM = find(xxData.Year == nYr & xxData.Month >= MMon(1) & xxData.Month <= MMon(end)); % absXMonEpmAnmaly(nYr-2011,:) = % end xQBOdata = dtTimeSeriesQBO(StartDateInd:EndDateInd); figure('position',[1,300,scrsz(3)*0.95,scrsz(4)*0.6]); subplot('position',[0.05,0.11,0.89,0.82]); p1 = scatter(datenum(xxData)-dn_start, absWinMonEpmAnomaly,'filled');hold on; % plot(1:1:108, zeros(1,length(1:1:108)),'--','LineWidth',1); ylim([-4 4]);ylabel('E_{pm} Anomaly (J/kg)');set(gca,'box','on'); text(datenum('20110101','yyyymmdd')-dn_start+110,-4.8,'2011','fontsize',16);text(datenum('20120101','yyyymmdd')-dn_start+110,-4.8,'2012','fontsize',16); text(datenum('20130101','yyyymmdd')-dn_start+110,-4.8,'2013','fontsize',16);text(datenum('20140101','yyyymmdd')-dn_start+110,-4.8,'2014','fontsize',16); text(datenum('20150101','yyyymmdd')-dn_start+110,-4.8,'2015','fontsize',16);text(datenum('20160101','yyyymmdd')-dn_start+110,-4.8,'2016','fontsize',16); text(datenum('20170101','yyyymmdd')-dn_start+110,-4.8,'2017','fontsize',16);text(datenum('20180101','yyyymmdd')-dn_start+110,-4.8,'2018','fontsize',16); text(datenum('20190101','yyyymmdd')-dn_start+110,-4.8,'2019','fontsize',16); yyaxis right p2 = plot(datenum(dtTimeSeriesQBO(StartDateInd):days(1):dtTimeSeriesQBO(EndDateInd))-dn_start, QBOIndex);hold on; plot(datenum(xQBOdata)-dn_start,zeros(length(xQBOdata),1),'--');hold on; [IndJJ1,~] = find(xQBOdata.Month == 6 & xQBOdata.Day == 15); [IndJJ2,~] = find(xQBOdata.Month == 7 & xQBOdata.Day == 15); p3 = scatter(datenum(xQBOdata(IndJJ1))-dn_start,QBOIndex(IndJJ1),'filled');hold on; scatter(datenum(xQBOdata(IndJJ2))-dn_start,QBOIndex(IndJJ2),'filled'); % hold on;plot(1:1:108, SingarporeData(2:end,7)); ylabel(['Zonal Mean Zonal Wind shear (m/s per 20hPa)']); ylim([-30 30]);xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';],'FontSize',15); legend([p1 p2 p3],{'June-July Epm anomaly',['Wind Shear between 50-70 hPa'],'June-July 15th'}); title('June-July Epm Anomaly versus QBO Index'); % corr QBO index and Epm % take QBO index mean % [IndJ1,~] = find(xQBOdata.Month == 6 | xQBOdata.Day >= 15); % [IndJ2,~] = find(xQBOdata.Month == 7 | xQBOdata.Day <= 15); % QBOIndex1 = QBOIndex(IndJ); % QBOIndexMean = []; % for nYr = 2011:1:2019 % IndTmp = find(xQBOdata(IndJ).Year == nYr); % Ind6 = IndTmp(1:30); % Ind7 = IndTmp(31:61); % if nYr == 2016 % QBOIndexMean(nYr-2010,1) = NaN; % QBOIndexMean(nYr-2010,2) = NaN; % else % QBOIndexMean(nYr-2010,1) = mean(QBOIndex1(Ind6),'all'); % QBOIndexMean(nYr-2010,2) = mean(QBOIndex1(Ind7),'all'); % end % end % QBOIndexMean = reshape(QBOIndexMean',[1,size(QBOIndexMean,1)*size(QBOIndexMean,2)]); % take single date [IndJ1,~] = find(xQBOdata.Month == 6 | xQBOdata.Day == 15); [IndJ2,~] = find(xQBOdata.Month == 7 | xQBOdata.Day == 15); QBOIndex1 = QBOIndex(IndJ1); QBOIndex2 = QBOIndex(IndJ2); QBOIndexMean = []; for nYr = 2011:1:2019 IndTmp1 = find(xQBOdata(IndJ1).Year == nYr); IndTmp2 = find(xQBOdata(IndJ2).Year == nYr); if nYr == 2016 QBOIndexMean(nYr-2010,1) = NaN; QBOIndexMean(nYr-2010,2) = NaN; elseif nYr == 2018 QBOIndexMean(nYr-2010,1) = NaN; QBOIndexMean(nYr-2010,2) = mean(QBOIndex2(IndTmp2),'all'); else QBOIndexMean(nYr-2010,1) = mean(QBOIndex1(IndTmp1),'all'); QBOIndexMean(nYr-2010,2) = mean(QBOIndex2(IndTmp2),'all'); end end QBOIndexMean = reshape(QBOIndexMean',[1,size(QBOIndexMean,1)*size(QBOIndexMean,2)]); IndNaN = ~isnan(absWinMonEpmAnomaly); Ind2NaN = ~isnan(QBOIndexMean); [Rabs Pabs]=corrcoef(absWinMonEpmAnomaly(IndNaN),QBOIndexMean(Ind2NaN)); [Rrel Prel]=corrcoef(relWinMonEpmAnomaly(IndNaN),QBOIndexMean(Ind2NaN)); figure(3); scatter(QBOIndexMean(Ind2NaN),absWinMonEpmAnomaly(IndNaN),'filled'); set(gca,'box','on');xlabel('QBO Index (m/s)');ylabel('Epm Anomaly (J/kg)'); title({'June-July Epm Anomaly Correlate With QBO Index '; ['coef = ' num2str(Rabs(1,2)) ' P = ' num2str(1-Pabs(1,2))] }); %% Epm and Epm anomaly versus Polar vortex edge distance - figure 4(b)(d) clear absEpmAnomaly; clear absMonEpmAnomaly; clear relEpmAnomaly; clear relMonEpmAnomaly; clear SingarporeWinData; clear uniformfit; EpmMonthlyMean = zeros(12,9); for year = 2011:2019 for month = 1:1:12 ind = find(timesort_datetime.Year == year & timesort_datetime.Month == month); EpmMonthlyMean(month, year-2010) = mean(mean30to50EpmTotal(ind)); EpmMonthlyMean_TimeSort(month, year-2010) = datetime(year,month,15); end end method =1; % method 1: using annual fitting curve as minuend % method 2: using monthly mean as minuend uniformfitMethod = 2;% 1 - 9 years uniformfit, 2 - composite fit % take mean for 9 years Epm in the same day of a year n = 1; CompositeEpm = cell(365,1); for nMon = 1:1:12 for nDay = 1:1:eomday(2011,nMon) Ind = timesort_datetime.Month == nMon & timesort_datetime.Day == nDay; CompositeEpm{n} = mean30to50EpmTotal(Ind); n = n+1; end end % take mean for every 5 days n = 1; binDay = 5; yData = zeros(365/binDay,1); xData = zeros(365/binDay,1); for nCnt = 1:binDay:365 yData(n) = mean([CompositeEpm{1+(n-1)*binDay} CompositeEpm{2+(n-1)*binDay} CompositeEpm{3+(n-1)*binDay} CompositeEpm{4+(n-1)*binDay} CompositeEpm{5+(n-1)*binDay}]); xData(n) = nCnt+2; n = n+1; end % fitting ftComposite = fitnlm(xData,yData','y ~ b1 + b2.*cos((2*pi/365).*(x1-b3))+b4.*cos((4*pi/365).*(x1-b5))',beta0); % calc Epm anomaly figure(1); switch method case 1 for nYr = 2011:1:2019 dnAll = (datenum(nYr,1,1)-dn_start):1:(datenum(nYr,12,31)-dn_start); AnnualFit = ft{nYr-2010}.Coefficients.Estimate(1) + ft{nYr-2010}.Coefficients.Estimate(2).*cos((2*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(3)))+ft{nYr-2010}.Coefficients.Estimate(4).*cos((4*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(5))); if uniformfitMethod == 1 Ind = dtSeries9yrs.Year == nYr; dn = TimeSeries9yrs(Ind); UniformFit = ft9yrs.Coefficients.Estimate(1)+ft9yrs.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ft9yrs.Coefficients.Estimate(3)))+ft9yrs.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ft9yrs.Coefficients.Estimate(5))); elseif uniformfitMethod == 2 Ind = dtSeries9yrs.Year == nYr; dn = TimeSeries9yrs(Ind); UniformFit = ftComposite.Coefficients.Estimate(1)+ftComposite.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ftComposite.Coefficients.Estimate(3)))+ftComposite.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ftComposite.Coefficients.Estimate(5))); end % plot(dnAll, AnnualFit);hold on; % plot(dn, UniformFit,'r');hold on; absEpmAnomaly{nYr-2010} = AnnualFit - UniformFit; plot(dn,absEpmAnomaly{nYr-2010},'b');hold on; relEpmAnomaly{nYr-2010} = absEpmAnomaly{nYr-2010}./UniformFit; % plot(dn,relEpmAnomaly{nYr-2010},'g');hold on; end n = 1; for nYr = 2011:1:2019 for nMon = 1:1:12 dtEachyear = datetime(nYr,1,1):1:datetime(nYr,12,31); Ind = find(dtEachyear.Year == nYr & dtEachyear.Month == nMon); absMonEpmAnomaly(n) = mean(absEpmAnomaly{nYr-2010}(Ind)); relMonEpmAnomaly(n) = mean(relEpmAnomaly{nYr-2010}(Ind)); maxMonEpmAnomaly(n) = max(absEpmAnomaly{nYr-2010}(Ind),[],'all'); dtMonTimeseries9yrs(n) = datetime(nYr,nMon,15); n = n+1; end end dnMonTimeseries9yrs = datenum(dtMonTimeseries9yrs)-dn_start; % scatter(dnMonTimeseries9yrs,absMonEpmAnomaly,'*');hold on; % scatter(dnMonTimeseries9yrs,relMonEpmAnomaly,'O'); case 2 n = 1; for nYr = 2011:1:2019 if uniformfitMethod == 1 Ind = dtSeries9yrs.Year == nYr; dn = TimeSeries9yrs(Ind); UniformFit = ft9yrs.Coefficients.Estimate(1)+ft9yrs.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ft9yrs.Coefficients.Estimate(3)))+ft9yrs.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ft9yrs.Coefficients.Estimate(5))); elseif uniformfitMethod == 2 Ind = dtSeries9yrs.Year == nYr; dn = TimeSeries9yrs(Ind); UniformFit = ftComposite.Coefficients.Estimate(1)+ftComposite.Coefficients.Estimate(2).*cos((2*pi/365).*(dn-ftComposite.Coefficients.Estimate(3)))+ftComposite.Coefficients.Estimate(4).*cos((4*pi/365).*(dn-ftComposite.Coefficients.Estimate(5))); end for nMon = 1:1:12 dtEachyear = datetime(nYr,1,1):1:datetime(nYr,12,31); IndYM = find(dtEachyear.Year == nYr & dtEachyear.Month == nMon); UniformMonFit(n) = mean(UniformFit(IndYM)); absMonEpmAnomaly(n) = EpmMonthlyMean(nMon,nYr-2010) - UniformMonFit(n); relMonEpmAnomaly(n) = absMonEpmAnomaly(n)/UniformMonFit(n); scatter(datenum(nYr,nMon,15)-dn_start,absMonEpmAnomaly(n),'*b');hold on; scatter(datenum(nYr,nMon,15)-dn_start,relMonEpmAnomaly(n),'o');hold on; n = n+1; end end dnMonTimeseries9yrs = datenum(dtMonTimeseries9yrs)-dn_start; plot(dnMonTimeseries9yrs, absMonEpmAnomaly);hold on; plot(dnMonTimeseries9yrs, relMonEpmAnomaly);hold on; plot(dnMonTimeseries9yrs, UniformMonFit); ylim([-5 5]); end path = 'C:\Users\l1989\Desktop\PolarVortex\FigureofPolarVortex\'; %Reading data % read new version data if ~exist('PolarVortexEdge', 'var') ttt = struct('DATE',{},'ALTITUDE',{},'EDGE_DISTANCE',{},'VORTEX_EDGE_WIND',{},'U',{},'V',{},'P',{},'T',{}); PolarVortexEdge = struct2cell(ttt); PolarVortexEdge.DATE = datetime(num2str(ncread([path 'merra2_at_mcmurdo_19800101-20200131.nc'],'DATES')),'InputFormat','yyyyMMdd'); PolarVortexEdge.ALTITUDE = ncread([path 'merra2_at_mcmurdo_19800101-20200131.nc'],'ALTITUDE'); PolarVortexEdge.EDGE_DISTANCE = ncread([path 'merra2_at_mcmurdo_19800101-20200131.nc'],'EDGE_DISTANCE'); PolarVortexEdge.VORTEX_EDGE_WIND = ncread([path 'merra2_at_mcmurdo_19800101-20200131.nc'],'VORTEX_EDGE_WIND'); PolarVortexEdge.U = ncread([path 'merra2_at_mcmurdo_19800101-20200131.nc'],'U'); PolarVortexEdge.V = ncread([path 'merra2_at_mcmurdo_19800101-20200131.nc'],'V'); PolarVortexEdge.P = ncread([path 'merra2_at_mcmurdo_19800101-20200131.nc'],'P'); PolarVortexEdge.T = ncread([path 'merra2_at_mcmurdo_19800101-20200131.nc'],'T'); end % alititude value, unit: km alt = 45; % time range StartDate = datetime(2011,1,1); EndDate = datetime(2019,12,31); durYear = ceil(years(EndDate-StartDate)); [StartDateInd,~] = find(PolarVortexEdge.DATE == StartDate); [EndDateInd,~] = find(PolarVortexEdge.DATE == EndDate); [~,altInd] = min(abs(PolarVortexEdge.ALTITUDE-alt)); n = 1; MonVortexEdgeDis = []; for nYr = StartDate.Year:1:EndDate.Year for nMon = 1:1:12 [DateInd,~] = find(PolarVortexEdge.DATE.Year == nYr & PolarVortexEdge.DATE.Month == nMon); MonVortexEdgeDis(n) = nanmean(PolarVortexEdge.EDGE_DISTANCE(DateInd,altInd)); n = n+1; end end set(0,'defaultfigurecolor','w'); scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3)*0.95,scrsz(4)*0.7]); subplot('position',[0.055,0.12,0.87,0.8]); yyaxis left % zoffset = 1; % t = hgtransform('Matrix',makehgtform('translate',0,0,zoffset)); p1 = errorbar(timesortTotal-dn_start,mean30to50EpmTotal,ErrorMean30to50EpmTotal,'*b','markersize',8); for nYr = 2011:1:2019 hold on; dnAll = (datenum(nYr,1,1)-dn_start):1:(datenum(nYr,12,31)-dn_start); p2 = plot(dnAll,ft{nYr-2010}.Coefficients.Estimate(1) + ft{nYr-2010}.Coefficients.Estimate(2).*cos((2*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(3)))+ft{nYr-2010}.Coefficients.Estimate(4).*cos((4*pi/365).*(dnAll-ft{nYr-2010}.Coefficients.Estimate(5))), '-r','LineWidth',2); end ax = ancestor(p1, 'axes'); ax.XAxis.FontSize = 15; %yrule = ax.YAxis;yrule.FontSize = 18; xlim([datenum('20110101','yyyymmdd')-dn_start,datenum('20191231','yyyymmdd')-dn_start]); ylim([0 20]); set(gca,'TickDir','out');set(gca,'Ycolor','b');set(gca,'box','on'); ylabel('Epm (J/kg)','fontsize',18); set(gca,'xtick',[datenum('20110101','yyyymmdd')-dn_start,datenum('20110201','yyyymmdd')-dn_start,datenum('20110301','yyyymmdd')-dn_start,... datenum('20110401','yyyymmdd')-dn_start,datenum('20110501','yyyymmdd')-dn_start,datenum('20110601','yyyymmdd')-dn_start,... datenum('20110701','yyyymmdd')-dn_start,datenum('20110801','yyyymmdd')-dn_start,datenum('20110901','yyyymmdd')-dn_start,... datenum('20111001','yyyymmdd')-dn_start,datenum('20111101','yyyymmdd')-dn_start,datenum('20111201','yyyymmdd')-dn_start,... datenum('20120101','yyyymmdd')-dn_start,datenum('20120201','yyyymmdd')-dn_start,datenum('20120301','yyyymmdd')-dn_start,... datenum('20120401','yyyymmdd')-dn_start,datenum('20120501','yyyymmdd')-dn_start,datenum('20120601','yyyymmdd')-dn_start,... datenum('20120701','yyyymmdd')-dn_start,datenum('20120801','yyyymmdd')-dn_start,datenum('20120901','yyyymmdd')-dn_start,... datenum('20121001','yyyymmdd')-dn_start,datenum('20121101','yyyymmdd')-dn_start,datenum('20121201','yyyymmdd')-dn_start,... datenum('20130101','yyyymmdd')-dn_start,datenum('20130201','yyyymmdd')-dn_start,datenum('20130301','yyyymmdd')-dn_start,... datenum('20130401','yyyymmdd')-dn_start,datenum('20130501','yyyymmdd')-dn_start,datenum('20130601','yyyymmdd')-dn_start,... datenum('20130701','yyyymmdd')-dn_start,datenum('20130801','yyyymmdd')-dn_start,datenum('20130901','yyyymmdd')-dn_start,... datenum('20131001','yyyymmdd')-dn_start,datenum('20131101','yyyymmdd')-dn_start,datenum('20131201','yyyymmdd')-dn_start,... datenum('20140101','yyyymmdd')-dn_start,datenum('20140201','yyyymmdd')-dn_start,datenum('20140301','yyyymmdd')-dn_start,... datenum('20140401','yyyymmdd')-dn_start,datenum('20140501','yyyymmdd')-dn_start,datenum('20140601','yyyymmdd')-dn_start,... datenum('20140701','yyyymmdd')-dn_start,datenum('20140801','yyyymmdd')-dn_start,datenum('20140901','yyyymmdd')-dn_start,... datenum('20141001','yyyymmdd')-dn_start,datenum('20141101','yyyymmdd')-dn_start,datenum('20141201','yyyymmdd')-dn_start,... datenum('20150101','yyyymmdd')-dn_start,datenum('20150201','yyyymmdd')-dn_start,datenum('20150301','yyyymmdd')-dn_start,... datenum('20150401','yyyymmdd')-dn_start,datenum('20150501','yyyymmdd')-dn_start,datenum('20150601','yyyymmdd')-dn_start,... datenum('20150701','yyyymmdd')-dn_start,datenum('20150801','yyyymmdd')-dn_start,datenum('20150901','yyyymmdd')-dn_start,... datenum('20151001','yyyymmdd')-dn_start,datenum('20151101','yyyymmdd')-dn_start,datenum('20151201','yyyymmdd')-dn_start,... datenum('20160101','yyyymmdd')-dn_start,datenum('20160201','yyyymmdd')-dn_start,datenum('20160301','yyyymmdd')-dn_start,... datenum('20160401','yyyymmdd')-dn_start,datenum('20160501','yyyymmdd')-dn_start,datenum('20160601','yyyymmdd')-dn_start,... datenum('20160701','yyyymmdd')-dn_start,datenum('20160801','yyyymmdd')-dn_start,datenum('20160901','yyyymmdd')-dn_start,... datenum('20161001','yyyymmdd')-dn_start,datenum('20161101','yyyymmdd')-dn_start,datenum('20161201','yyyymmdd')-dn_start,... datenum('20170101','yyyymmdd')-dn_start,datenum('20170201','yyyymmdd')-dn_start,datenum('20170301','yyyymmdd')-dn_start,... datenum('20170401','yyyymmdd')-dn_start,datenum('20170501','yyyymmdd')-dn_start,datenum('20170601','yyyymmdd')-dn_start,... datenum('20170701','yyyymmdd')-dn_start,datenum('20170801','yyyymmdd')-dn_start,datenum('20170901','yyyymmdd')-dn_start,... datenum('20171001','yyyymmdd')-dn_start,datenum('20171101','yyyymmdd')-dn_start,datenum('20171201','yyyymmdd')-dn_start,... datenum('20180101','yyyymmdd')-dn_start,datenum('20180201','yyyymmdd')-dn_start,datenum('20180301','yyyymmdd')-dn_start,... datenum('20180401','yyyymmdd')-dn_start,datenum('20180501','yyyymmdd')-dn_start,datenum('20180601','yyyymmdd')-dn_start,... datenum('20180701','yyyymmdd')-dn_start,datenum('20180801','yyyymmdd')-dn_start,datenum('20180901','yyyymmdd')-dn_start,... datenum('20181001','yyyymmdd')-dn_start,datenum('20181101','yyyymmdd')-dn_start,datenum('20181201','yyyymmdd')-dn_start,... datenum('20190101','yyyymmdd')-dn_start,datenum('20190201','yyyymmdd')-dn_start,datenum('20190301','yyyymmdd')-dn_start,... datenum('20190401','yyyymmdd')-dn_start,datenum('20190501','yyyymmdd')-dn_start,datenum('20190601','yyyymmdd')-dn_start,... datenum('20190701','yyyymmdd')-dn_start,datenum('20190801','yyyymmdd')-dn_start,datenum('20190901','yyyymmdd')-dn_start,... datenum('20191001','yyyymmdd')-dn_start,datenum('20191101','yyyymmdd')-dn_start,datenum('20191201','yyyymmdd')-dn_start,],'xticklabel',... ['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';... 'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';'J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D';]); text(datenum('20110101','yyyymmdd')-dn_start+110,-2,'2011','fontsize',18);text(datenum('20120101','yyyymmdd')-dn_start+110,-2,'2012','fontsize',18); text(datenum('20130101','yyyymmdd')-dn_start+110,-2,'2013','fontsize',18);text(datenum('20140101','yyyymmdd')-dn_start+110,-2,'2014','fontsize',18); text(datenum('20150101','yyyymmdd')-dn_start+110,-2,'2015','fontsize',18);text(datenum('20160101','yyyymmdd')-dn_start+110,-2,'2016','fontsize',18); text(datenum('20170101','yyyymmdd')-dn_start+110,-2,'2017','fontsize',18);text(datenum('20180101','yyyymmdd')-dn_start+110,-2,'2018','fontsize',18); text(datenum('20190101','yyyymmdd')-dn_start+110,-2,'2019','fontsize',18); yyaxis right xData = datenum(PolarVortexEdge.DATE(StartDateInd):days(1):PolarVortexEdge.DATE(EndDateInd))-dn_start; yData = PolarVortexEdge.EDGE_DISTANCE(StartDateInd:EndDateInd,altInd)/1000; % temp = smooth(yData,7,'rlowess'); % IndNan = ~isnan(yData); % [Ind,~] = find(IndNan == 1); p3 = plot(xData,yData,'LineWidth',1.5);hold on; % p3 = % plot(xData,smoothdata(yData,'rlowess',7,'omitnan'),'LineWidth',1.5);hold % on;% 7 points smooth ax = gca; % set(p1,'parent',ax); [~,IndApr] = find(dtMonTimeseries9yrs.Month == 4); [~,IndMay] = find(dtMonTimeseries9yrs.Month == 5); [~,IndJunExclude] = find(dtMonTimeseries9yrs.Month == 6 & dtMonTimeseries9yrs.Year ~= 2018);%& dtMonTimeseries9yrs.Year ~= 2018 [~,IndJulExclude] = find(dtMonTimeseries9yrs.Month == 7 & dtMonTimeseries9yrs.Year ~= 2018); [~,IndJun2018] = find(dtMonTimeseries9yrs.Month == 6 & dtMonTimeseries9yrs.Year == 2018); [~,IndJun] = find(dtMonTimeseries9yrs.Month == 6);%& dtMonTimeseries9yrs.Year ~= 2018 [~,IndJul] = find(dtMonTimeseries9yrs.Month == 7); [~,IndJJ] = find(dtMonTimeseries9yrs.Month == 7 | dtMonTimeseries9yrs.Month == 6); % IndT = [Ind1 Ind2]; p4 = scatter(dnMonTimeseries9yrs(IndJJ),MonVortexEdgeDis(IndJJ)/1000,'filled','MarkerFaceColor','g');hold on; % set(gca, 'SortMethod', 'depth'); ylim([-12 4]);ylabel('Edge Distance (10^3 km)'); title(['Epm versus Polar Vortex Edge Distance @' num2str(alt) ' km'],'position',[(datenum('20150701','yyyymmdd')-dn_start),20.6],'fontsize',20); % legend([p1 p2 p3 p4],{'Epm';'Annual Fit';'Edge Distance';'June - July'}); % legend('boxoff'); n = 1; for nYr = StartDate.Year:1:EndDate.Year [DateInd1,~] = find(PolarVortexEdge.DATE.Year == nYr & PolarVortexEdge.DATE.Month == 5 & PolarVortexEdge.DATE.Month >= 10); [DateInd2,~] = find(PolarVortexEdge.DATE.Year == nYr & PolarVortexEdge.DATE.Month == 6 & PolarVortexEdge.DATE.Month <= 10); IndTT = [DateInd1 DateInd2]; VortexEdgeDis(n) = nanmean(PolarVortexEdge.EDGE_DISTANCE(IndTT,altInd)); n = n+1; end scrsz = get(0, 'ScreenSize'); figure('position',[1,300,scrsz(3)*0.39,scrsz(4)*0.55]); subplot('position',[0.13,0.16,0.78,0.7]); [R1,P1] = corrcoef(MonVortexEdgeDis(IndJun)/1000,absMonEpmAnomaly(IndJun)); myFit = NonLinearModel.fit(MonVortexEdgeDis(IndJun)/1000,absMonEpmAnomaly(IndJun), 'y ~ b0 + b1*x', [0,1]); [p,S]= polyfit(MonVortexEdgeDis(IndJun)/1000,absMonEpmAnomaly(IndJun),1); [f,delta]= polyval(p,MonVortexEdgeDis(IndJun)/1000,S); s1 = plot(MonVortexEdgeDis(IndJun)/1000,f,'-r','LineWidth',2);hold on; [R2,P2] = corrcoef(MonVortexEdgeDis(IndJunExclude)/1000,absMonEpmAnomaly(IndJunExclude)); myFit1 = NonLinearModel.fit(MonVortexEdgeDis(IndJunExclude)/1000,absMonEpmAnomaly(IndJunExclude), 'y ~ b0 + b1*x', [0,1]); [p1,S1]= polyfit(MonVortexEdgeDis(IndJunExclude)/1000,absMonEpmAnomaly(IndJunExclude),1); [f1,delta1]= polyval(p1,(min(MonVortexEdgeDis(IndJunExclude)):1:max(MonVortexEdgeDis(IndJunExclude)))/1000,S1); s2 = plot([min(MonVortexEdgeDis(IndJunExclude)):1:max(MonVortexEdgeDis(IndJunExclude))]/1000,f1,'--r','LineWidth',2);hold on;set(gca,'fontsize',18); plot(MonVortexEdgeDis(IndJun)/1000,absMonEpmAnomaly(IndJun),'LineStyle','none','Marker','o','MarkerSize',8,'MarkerFaceColor','b');hold on; plot(MonVortexEdgeDis(IndJun2018)/1000,absMonEpmAnomaly(IndJun2018),'linewidth',1.5,'Marker','o','MarkerSize',15,'MarkerFaceColor','none'); xlabel('Vortex Edge Dis (10^3 km)');ylabel('Epm Anomaly (J/kg)');xlim([-4.6 -2.8]); % legend([s1 s2],{['corr = ' num2str(R1(1,2)) 'ConfLev = ' num2str(1-P1(1,2))];['corr = ' num2str(R2(1,2)) 'ConfLev = ' num2str(1-P2(1,2)) ' - without 2018']}); title({'June Epm Anomaly vs. June Vortex Edge Distance';['Correlation Coeff = ' num2str(roundn(R1(1,2),-2)) ' / ' num2str(roundn(R2(1,2),-2)) ' P = ' num2str(roundn(1-P1(1,2),-2)) ' / ' num2str(roundn(1-P2(1,2),-2))]}); % title({'July Epm Anomaly corr with June Polar Vortex Edge Dis';['corr = ' num2str(R1(1,2)) 'ConfLev = ' num2str(1-P1(1,2))]}); figure(4); [R3,P3] = corrcoef([MonVortexEdgeDis(IndJun)/1000 MonVortexEdgeDis(IndJun)/1000],[absMonEpmAnomaly(IndJun) absMonEpmAnomaly(IndJul)]); scatter(MonVortexEdgeDis(IndJun)/1000,absMonEpmAnomaly(IndJun),'b');hold on; scatter(MonVortexEdgeDis(IndJun)/1000,absMonEpmAnomaly(IndJul),'r'); title({'June - July Epm Anomaly corr with June Polar Vortex Edge Dis';['corr = ' num2str(R3(1,2)) 'ConfLev = ' num2str(1-P3(1,2))]});