%% %For D23 data p1x=647; p1y=938; % 130.8948 32.8112 p2x=938; p2y=737; % 131.0501 32.9225 p3x=827; p3y=892; % 130.9908 32.8367 p4x=499; p4y=915; % 130.8158 32.8239 p5x=493; p5y=999; % 130.8126 32.7774 p6x=389; p6y=1101; % 130.7571 32.7209 p7x=185; p7y=1021; % 130.6482 32.7652 p8x=758; p8y=1185; % 130.9540 32.6744 figure,surf(RefK, RnfK,intD23_29502960);shading flat; view(0,90);... axis equal;colormap jet;caxis([-10 10]);hold on; scatter3( RefK(p1y, p1x), RnfK(p1y, p1x), 100, 'filled', '^k'); scatter3(RefK(p2y, p2x), RnfK(p2y, p2x), 100, 'filled', '^k'); scatter3(RefK(p3y, p3x), RnfK(p3y, p3x), 100, 'filled', '^k'); scatter3( RefK(p4y, p4x), RnfK(p4y, p4x), 100, 'filled', '^k'); scatter3( RefK(p5y, p5x), RnfK(p5y, p5x), 100, 'filled', '^k'); scatter3( RefK(p6y, p6x), RnfK(p6y, p6x), 100, 'filled', '^k'); scatter3( RefK(p7y, p7x), RnfK(p7y, p7x), 100, 'filled', '^k'); timeD23=[ (datenum('2016-05-16', 'yyyy-mm-dd')-datenum('2016-04-18', 'yyyy-mm-dd'))/2+datenum('2016-04-18', 'yyyy-mm-dd'); (datenum('2016-06-13', 'yyyy-mm-dd')-datenum('2016-05-16', 'yyyy-mm-dd'))/2+datenum('2016-05-16', 'yyyy-mm-dd'); (datenum('2016-07-11', 'yyyy-mm-dd')-datenum('2016-06-13', 'yyyy-mm-dd'))/2+datenum('2016-06-13', 'yyyy-mm-dd'); (datenum('2016-08-08', 'yyyy-mm-dd')-datenum('2016-07-11', 'yyyy-mm-dd'))/2+datenum('2016-07-11', 'yyyy-mm-dd'); (datenum('2016-09-19', 'yyyy-mm-dd')-datenum('2016-08-08', 'yyyy-mm-dd'))/2+datenum('2016-08-08', 'yyyy-mm-dd'); (datenum('2016-11-14', 'yyyy-mm-dd')-datenum('2016-09-19', 'yyyy-mm-dd'))/2+datenum('2016-09-19', 'yyyy-mm-dd'); (datenum('2017-11-13', 'yyyy-mm-dd')-datenum('2016-11-14', 'yyyy-mm-dd'))/2+datenum('2016-11-14', 'yyyy-mm-dd'); (datenum('2018-03-05', 'yyyy-mm-dd')-datenum('2017-11-13', 'yyyy-mm-dd'))/2+datenum('2017-11-13', 'yyyy-mm-dd'); (datenum('2018-08-20', 'yyyy-mm-dd')-datenum('2018-03-05', 'yyyy-mm-dd'))/2+datenum('2018-03-05', 'yyyy-mm-dd'); (datenum('2018-10-15', 'yyyy-mm-dd')-datenum('2018-08-20', 'yyyy-mm-dd'))/2+datenum('2018-08-20', 'yyyy-mm-dd');]; timeD23dif=timeD23-736436; err_ran=15; %% p1data=[ D23int20160418_20160516a(p1y, p1x); D23int20160516_20160613a(p1y, p1x); D23int20160613_20160711a(p1y, p1x); D23int20160711_20160808a(p1y, p1x); D23int20160808_20160919a(p1y, p1x); D23int20160919_20161114a(p1y, p1x); D23int20161114_20171113a(p1y, p1x); D23int20171113_20180305a(p1y, p1x); D23int20180305_20180820a(p1y, p1x); D23int20180820_20181015a(p1y, p1x)]; p1cum=cumsum(p1data); p1_err=[ std(std((D23int20160418_20160516a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((D23int20160516_20160613a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((D23int20160613_20160711a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((D23int20160711_20160808a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((D23int20160808_20160919a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((D23int20160919_20161114a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((D23int20161114_20171113a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((D23int20171113_20180305a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((D23int20180305_20180820a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((D23int20180820_20181015a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran))));]; t=timeD23dif;y1=p1cum; p1_fittype1=fittype('a*log(1+(t/tau))', 'dependent', {'y1'}, 'independent', {'t'}, ... 'coefficients', {'a', 'tau'}); %p1_fittype2=fittype('a*log(1+(t/tau))+b*exp(-t/taug)+c', 'dependent', {'y1'}, 'independent', {'t'}, ... % 'coefficients', {'a', 'tau', 'b','taug', 'c'}); %startPoints=[-3 0.1 -3 1 0.1];ExcludeP=[100 1000 100 1000]; myfit11=fit(t, y1, p1_fittype1, 'Lower',[-100 0.000001] , 'Upper', [100 1000] ) %myfit12=fit(t, y1, p1_fittype2, 'Lower',[-3 0.1 -3 0.1 -100] , 'Upper', [100 1000 100 1000 100] ) SSresid1=sum((myfit11(t)-p1cum).^2); %SSresid2=sum((myfit12(t)-p1cum).^2); SStotal=(length(t)-1)*var(p1cum); rsq11=1-SSresid1/SStotal time=linspace(1,1000, 1000); fit11=myfit11(time); timeP=linspace(736436, 737425,1000); %rsq12=1-SSresid2/SStotal %figure, %plot(myfit11, 'r', t, y1);hold on;xlim([5 1000]); %plot(myfit12, 'b', t, y1); %plot(myfit12,'b', t, y1); %% %Point at Aso volcano p2data=[ D23int20160418_20160516a(p2y, p2x); D23int20160516_20160613a(p2y, p2x); D23int20160613_20160711a(p2y, p2x); D23int20160711_20160808a(p2y, p2x); D23int20160808_20160919a(p2y, p2x); D23int20160919_20161114a(p2y, p2x); D23int20161114_20171113a(p2y, p2x); D23int20171113_20180305a(p2y, p2x); D23int20180305_20180820a(p2y, p2x); D23int20180820_20181015a(p2y, p2x)]; p2cum=cumsum(p2data); p2_err=[ std(std((D23int20160418_20160516a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); std(std((D23int20160516_20160613a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); std(std((D23int20160613_20160711a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); std(std((D23int20160711_20160808a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); std(std((D23int20160808_20160919a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); std(std((D23int20160919_20161114a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); std(std((D23int20161114_20171113a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); std(std((D23int20171113_20180305a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); std(std((D23int20180305_20180820a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); std(std((D23int20180820_20181015a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran))));]; t=timeD23dif;y2=p2cum; p2_fittype1=fittype('a*log(1+(t/tau))', 'dependent', {'y1'}, 'independent', {'t'}, ... 'coefficients', {'a', 'tau'}); %p1_fittype2=fittype('a*log(1+(t/tau))+b*exp(-t/taug)+c', 'dependent', {'y1'}, 'independent', {'t'}, ... % 'coefficients', {'a', 'tau', 'b','taug', 'c'}); %startPoints=[-3 0.1 -3 1 0.1];ExcludeP=[100 1000 100 1000]; myfit21=fit(t, y2, p2_fittype1, 'Lower',[-100 0.000001] , 'Upper', [100 1000] ) %myfit22=fit(t, y1, p1_fittype2, 'Lower',[-3 0.1 -3 0.1 -100] , 'Upper', [100 1000 100 1000 100] ) SSresid21=sum((myfit21(t)-p2cum).^2); SStotal2=(length(t)-1)*var(p2cum); rsq21=1-SSresid21/SStotal2 time=linspace(1,1000, 1000); fit21=myfit21(time); % figure,scatter(timeD23, p2cum, 54, 'r', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2);ylim([0 9]);grid on;hold on; % plot([736436 736436], [-10 10], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% %Point at north flank p3data=[ D23int20160418_20160516a(p3y, p3x); D23int20160516_20160613a(p3y, p3x); D23int20160613_20160711a(p3y, p3x); D23int20160711_20160808a(p3y, p3x); D23int20160808_20160919a(p3y, p3x); D23int20160919_20161114a(p3y, p3x); D23int20161114_20171113a(p3y, p3x); D23int20171113_20180305a(p3y, p3x); D23int20180305_20180820a(p3y, p3x); D23int20180820_20181015a(p3y, p3x)]; p3cum=cumsum(p3data); p3_err=[ std(std((D23int20160418_20160516a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((D23int20160516_20160613a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((D23int20160613_20160711a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((D23int20160711_20160808a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((D23int20160808_20160919a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((D23int20160919_20161114a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((D23int20161114_20171113a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((D23int20171113_20180305a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((D23int20180305_20180820a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((D23int20180820_20181015a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran))));]; t=timeD23dif;y3=p3cum; p3_fittype1=fittype('a*log(1+(t/tau))', 'dependent', {'y3'}, 'independent', {'t'}, ... 'coefficients', {'a', 'tau'}); p3_fittype2=fittype('b*t+c', 'dependent', {'y3'}, 'independent', {'t'}, ... 'coefficients', {'b', 'c'}); %startPoints=[-3 0.1 -3 1 0.1];ExcludeP=[100 1000 100 1000]; myfit31=fit(t, y3, p3_fittype1, 'Lower',[-100 0.000001] , 'Upper', [100 1000] ) myfit32=fit(t, y3, p3_fittype2, 'Lower',[-30 -100] , 'Upper', [1000 1000] ) SSresid31=sum((myfit31(t)-p3cum).^2); SSresid32=sum((myfit32(t)-p3cum).^2); SStotal3=(length(t)-1)*var(p3cum); rsq31=1-SSresid31/SStotal3 rsq32=1-SSresid32/SStotal3 time=linspace(1,1000, 1000); fit31=myfit31(time); fit32=myfit32(time); % figure,scatter(timeD23, p3cum, 54, 'r', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2);ylim([-9 0]);grid on;hold on; % plot([736436 736436], [-10 10], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% %Point at south flank p4data=[ D23int20160418_20160516a(p4y, p4x); D23int20160516_20160613a(p4y, p4x); D23int20160613_20160711a(p4y, p4x); D23int20160711_20160808a(p4y, p4x); D23int20160808_20160919a(p4y, p4x); D23int20160919_20161114a(p4y, p4x); D23int20161114_20171113a(p4y, p4x); D23int20171113_20180305a(p4y, p4x); D23int20180305_20180820a(p4y, p4x); D23int20180820_20181015a(p4y, p4x)]; p4cum=cumsum(p4data); p4_err=[ std(std((D23int20160418_20160516a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((D23int20160516_20160613a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((D23int20160613_20160711a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((D23int20160711_20160808a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((D23int20160808_20160919a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((D23int20160919_20161114a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((D23int20161114_20171113a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((D23int20171113_20180305a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((D23int20180305_20180820a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((D23int20180820_20181015a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran))));]; t=timeD23dif;y4=p4cum; p4_fittype1=fittype('a*log(1+(t/tau))', 'dependent', {'y4'}, 'independent', {'t'}, ... 'coefficients', {'a', 'tau'}); %p1_fittype4=fittype('a*log(1+(t/tau))+b*exp(-t/taug)+c', 'dependent', {'y1'}, 'independent', {'t'}, ... % 'coefficients', {'a', 'tau', 'b','taug', 'c'}); %startPoints=[-3 0.1 -3 1 0.1];ExcludeP=[100 1000 100 1000]; myfit41=fit(t, y4, p4_fittype1, 'Lower',[-300 0.000001] , 'Upper', [100 1000] ) %myfit42=fit(t, y1, p1_fittype4, 'Lower',[-3 0.1 -3 0.1 -100] , 'Upper', [100 1000 100 1000 100] ) SSresid41=sum((myfit41(t)-p4cum).^2); SStotal4=(length(t)-1)*var(p4cum); rsq41=1-SSresid41/SStotal4 time=linspace(1,1000, 1000); fit41=myfit41(time); % figure,scatter(timeD23, p4cum, 54, 'r', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2);ylim([-9 0]);grid on;hold on; % plot(timeP, fit41); % plot([736436 736436], [-10 10], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% %Point at fault junction p5data=[ D23int20160418_20160516a(p5y, p5x); D23int20160516_20160613a(p5y, p5x); D23int20160613_20160711a(p5y, p5x); D23int20160711_20160808a(p5y, p5x); D23int20160808_20160919a(p5y, p5x); D23int20160919_20161114a(p5y, p5x); D23int20161114_20171113a(p5y, p5x); D23int20171113_20180305a(p5y, p5x); D23int20180305_20180820a(p5y, p5x); D23int20180820_20181015a(p5y, p5x)]; p5cum=cumsum(p5data); p5_err=[ std(std((D23int20160418_20160516a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((D23int20160516_20160613a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((D23int20160613_20160711a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((D23int20160711_20160808a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((D23int20160808_20160919a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((D23int20160919_20161114a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((D23int20161114_20171113a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((D23int20171113_20180305a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((D23int20180305_20180820a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((D23int20180820_20181015a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran))));]; t=timeD23dif;y5=p5cum; p5_fittype1=fittype('a*log(1+(t/tau))', 'dependent', {'y5'}, 'independent', {'t'}, ... 'coefficients', {'a', 'tau'}); % p5_fittype2=fittype('a*log(1+(t/tau))+b*exp(-t/taug)', 'dependent', {'y1'}, 'independent', {'t'}, ... % 'coefficients', {'a', 'tau', 'b','taug'}); %startPoints=[-3 0.1 -3 1 0.1];ExcludeP=[100 1000 100 1000]; myfit51=fit(t, y5, p5_fittype1, 'Lower',[-100 0.000001 ] , 'Upper', [100 1000 ] ) % myfit52=fit(t, y5, p5_fittype2, 'Lower',[0 0.000001 0 0.000001 ] , 'Upper', [10000 1000 10000 100000 ] ) SSresid51=sum((myfit51(t)-p5cum).^2); %SSresid52=sum((myfit52(t)-p5cum).^2); SStotal5=(length(t)-1)*var(p5cum); rsq51=1-SSresid51/SStotal5 %rsq52=1-SSresid52/SStotal5 time=linspace(1,1000, 1000); fit51=myfit51(time); %fit52=myfit52(time); % figure,scatter(timeD23, p5cum, 54, 'r', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2);ylim([0 9]);grid on;hold on; % plot(timeP, fit51, 'r');hold on; % plot(timeP, fit52, 'b');hold on; % plot([736436 736436], [-10 10], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% %Point at Kumamoto delta p6data=[ D23int20160418_20160516a(p6y, p6x); D23int20160516_20160613a(p6y, p6x); D23int20160613_20160711a(p6y, p6x); D23int20160711_20160808a(p6y, p6x); D23int20160808_20160919a(p6y, p6x); D23int20160919_20161114a(p6y, p6x); D23int20161114_20171113a(p6y, p6x); D23int20171113_20180305a(p6y, p6x); D23int20180305_20180820a(p6y, p6x); D23int20180820_20181015a(p6y, p6x)]; p6cum=cumsum(p6data); p6_err=[ std(std((D23int20160418_20160516a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((D23int20160516_20160613a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((D23int20160613_20160711a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((D23int20160711_20160808a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((D23int20160808_20160919a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((D23int20160919_20161114a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((D23int20161114_20171113a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((D23int20171113_20180305a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((D23int20180305_20180820a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((D23int20180820_20181015a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran))));]; t=timeD23dif;y6=p6cum; p6_fittype1=fittype('a*log(1+(t/tau))', 'dependent', {'y6'}, 'independent', {'t'}, ... 'coefficients', {'a', 'tau'}); p6_fittype2=fittype('b*t+c', 'dependent', {'y1'}, 'independent', {'t'}, ... 'coefficients', {'b', 'c'}); %startPoints=[-3 0.1 -3 1 0.1];ExcludeP=[100 1000 100 1000]; myfit61=fit(t, y6, p6_fittype1, 'Lower',[-100 0.000001 ] , 'Upper', [100 1000 ] ) myfit62=fit(t, y6, p6_fittype2, 'Lower',[-100 -10] , 'Upper', [100 10 ] ) SSresid61=sum((myfit61(t)-p6cum).^2); SSresid62=sum((myfit62(t)-p6cum).^2); SStotal6=(length(t)-1)*var(p6cum); rsq61=1-SSresid61/SStotal6 rsq62=1-SSresid62/SStotal6 time=linspace(1,1000, 1000); fit61=myfit61(time); fit62=myfit62(time); % % figure,scatter(timeD23, p6cum, 54, 'r', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2);ylim([-9 0]);grid on;hold on; % plot(timeP, fit61, 'r');hold on; % %plot(timeP, fit62, 'b');hold on; % plot([736436 736436], [-10 10], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% %Point at Kumamoto city p7data=[ D23int20160418_20160516a(p7y, p7x); D23int20160516_20160613a(p7y, p7x); D23int20160613_20160711a(p7y, p7x); D23int20160711_20160808a(p7y, p7x); D23int20160808_20160919a(p7y, p7x); D23int20160919_20161114a(p7y, p7x); D23int20161114_20171113a(p7y, p7x); D23int20171113_20180305a(p7y, p7x); D23int20180305_20180820a(p7y, p7x); D23int20180820_20181015a(p7y, p7x)]; p7cum=cumsum(p7data); p7_err=[ std(std((D23int20160418_20160516a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((D23int20160516_20160613a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((D23int20160613_20160711a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((D23int20160711_20160808a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((D23int20160808_20160919a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((D23int20160919_20161114a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((D23int20161114_20171113a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((D23int20171113_20180305a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((D23int20180305_20180820a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((D23int20180820_20181015a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran))));]; t=timeD23dif;y7=p7cum; p7_fittype1=fittype('a*log(1+(t/tau))', 'dependent', {'y7'}, 'independent', {'t'}, ... 'coefficients', {'a', 'tau'}); p7_fittype2=fittype('b*t+c', 'dependent', {'y7'}, 'independent', {'t'}, ... 'coefficients', {'b','c'}); %startPoints=[-3 0.1 -3 1 0.1];ExcludeP=[100 1000 100 1000]; myfit71=fit(t, y7, p7_fittype1, 'Lower',[-30 0.000001] , 'Upper', [10 10000 ] ) myfit72=fit(t, y7, p7_fittype2, 'Lower',[-30 -30 ] , 'Upper', [100 100 ], 'Start', [0 0] ) SSresid71=sum((myfit71(t)-p7cum).^2); SSresid72=sum((myfit72(t)-p7cum).^2); SStotal7=(length(t)-1)*var(p7cum); rsq71=1-SSresid71/SStotal7 rsq72=1-SSresid72/SStotal7 time=linspace(1,1000, 1000); fit71=myfit71(time); %fit72=myfit72(time); % figure,scatter(timeD23, p7cum, 54, 'r', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2);ylim([0 9]);grid on;hold on; % plot(timeP, fit71, 'r');hold on; % %plot(timeP, fit72, 'b');hold on; % plot([736436 736436], [-10 10], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% p8data=[ D23int20160418_20160516a(p8y, p8x); D23int20160516_20160613a(p8y, p8x); D23int20160613_20160711a(p8y, p8x); D23int20160711_20160808a(p8y, p8x); D23int20160808_20160919a(p8y, p8x); D23int20160919_20161114a(p8y, p8x); D23int20161114_20171113a(p8y, p8x); D23int20171113_20180305a(p8y, p8x); D23int20180305_20180820a(p8y, p8x); D23int20180820_20181015a(p8y, p8x)]; p8cum=cumsum(p8data); p8_err=[ std(std((D23int20160418_20160516a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((D23int20160516_20160613a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((D23int20160613_20160711a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((D23int20160711_20160808a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((D23int20160808_20160919a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((D23int20160919_20161114a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((D23int20161114_20171113a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((D23int20171113_20180305a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((D23int20180305_20180820a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((D23int20180820_20181015a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran))));]; t=timeD23dif;y8=p8cum; p8_fittype1=fittype('a*log(1+(t/tau))', 'dependent', {'y8'}, 'independent', {'t'}, ... 'coefficients', {'a', 'tau'}); p8_fittype2=fittype('b*t+c', 'dependent', {'y8'}, 'independent', {'t'}, ... 'coefficients', {'b', 'c'}); %startPoints=[-3 0.1 -3 1 0.1];ExcludeP=[100 1000 100 1000]; myfit81=fit(t, y8, p8_fittype1, 'Lower',[-30 0.000001 ] , 'Upper', [10 10000 ] ) myfit82=fit(t, y8, p8_fittype2, 'Lower',[-30 -30 ] , 'Upper', [30 30 ] ) SSresid81=sum((myfit81(t)-p8cum).^2); SSresid82=sum((myfit82(t)-p8cum).^2); SStotal8=(length(t)-1)*var(p8cum); rsq81=1-SSresid81/SStotal8 rsq82=1-SSresid82/SStotal8 time=linspace(1,1000, 1000); fit81=myfit81(time); fit82=myfit82(time); % figure,scatter(timeD23, p8cum, 54, 'r', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2);ylim([-9 3]);grid on;hold on; % plot(timeP, fit81, 'r');hold on; % plot(timeP, fit82, 'b');hold on; % plot([736436 736436], [-10 10], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% figure, subplot(4,2,1); errorbar(timeD23, p1cum, p1_err, 'ko', 'LineWidth', 2);ylim([-3 11]);grid on;hold on; plot(timeP, fit11, 'r', 'LineWidth', 3);hold on; scatter(timeD23, p1cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); subplot(4,2,3); errorbar(timeD23, p2cum, p2_err, 'ko', 'LineWidth', 2);ylim([-3 11]);grid on;hold on; plot(timeP, fit21, 'r', 'LineWidth', 3);hold on; scatter(timeD23, p2cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); subplot(4,2,5); errorbar(timeD23, p3cum, p3_err, 'ko', 'LineWidth', 2);ylim([-11 3]);grid on;hold on; plot(timeP, fit31, 'r', 'LineWidth', 3);hold on; plot(timeP, fit32, 'b', 'LineWidth', 1);hold on; scatter(timeD23, p3cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); subplot(4,2,7); errorbar(timeD23, p4cum, p4_err, 'ko', 'LineWidth', 2);ylim([-11 3]);grid on;hold on; plot(timeP, fit41, 'r', 'LineWidth', 3);hold on; scatter(timeD23, p4cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); subplot(4,2,2); errorbar(timeD23, p5cum, p5_err, 'ko', 'LineWidth', 2);ylim([-3 11]);grid on;hold on; plot(timeP, fit51, 'r', 'LineWidth', 3);hold on; scatter(timeD23, p5cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); subplot(4,2,4); errorbar(timeD23, p6cum, p6_err, 'ko', 'LineWidth', 2);ylim([-11 3]);grid on;hold on; plot(timeP, fit61, 'r', 'LineWidth', 3);hold on; plot(timeP, fit62, 'b', 'LineWidth', 1);hold on; scatter(timeD23, p6cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); subplot(4,2,6); errorbar(timeD23, p7cum, p7_err, 'ko', 'LineWidth', 2);ylim([-3 11]);grid on;hold on; plot(timeP, fit71, 'r', 'LineWidth', 3);hold on; scatter(timeD23, p7cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); subplot(4,2,8); errorbar(timeD23, p8cum, p8_err, 'ko', 'LineWidth', 2);ylim([-11 3]);grid on;hold on; plot(timeP, fit81, 'r', 'LineWidth', 3);hold on; plot(timeP, fit82, 'b', 'LineWidth', 1);hold on; scatter(timeD23, p8cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); %% rsqd=[rsq11; rsq21; rsq31; rsq41; rsq51; rsq61; rsq62; rsq71; rsq72]; rsqdtag=['rsqd11='; 'rsqd21='; 'rsqd31='; 'rsqd41='; 'rsqd51='; 'rsqd61='; 'rsqd62='; 'rsqd71='; 'rsqd72=']; disp([rsqdtag num2str(rsqd)]); disp(myfit11); disp(myfit21); disp(myfit31); disp(myfit41); disp(myfit51); disp(myfit61); disp(myfit71); %% [tv,pvd11]=ttest(myfit11(t)-p1cum) [tv,pvd21]=ttest(myfit21(t)-p2cum) [tv,pvd31]=ttest(myfit31(t)-p3cum) [tv,pvd41]=ttest(myfit41(t)-p4cum) [tv,pvd51]=ttest(myfit51(t)-p5cum) [tv,pvd61]=ttest(myfit61(t)-p6cum) [tv,pvd71]=ttest(myfit71(t)-p7cum) [tv,pvd81]=ttest(myfit81(t)-p8cum) disp(tv); %%EOF %% %For A131 data p1x=647; p1y=938; % 130.8948 32.8112 p2x=938; p2y=737; % 131.0501 32.9225 p3x=827; p3y=892; % 130.9908 32.8367 p4x=499; p4y=915; % 130.8158 32.8239 p5x=493; p5y=999; % 130.8126 32.7774 p6x=389; p6y=1101; % 130.7571 32.7209 p7x=185; p7y=1021; % 130.6482 32.7652 p8x=758; p8y=1185; % 130.9540 32.6744 timeA131=[ (datenum('2016-05-05', 'yyyy-mm-dd')-datenum('2016-04-21', 'yyyy-mm-dd'))/2+datenum('2016-05-05', 'yyyy-mm-dd'); (datenum('2016-07-19', 'yyyy-mm-dd')-datenum('2016-05-05', 'yyyy-mm-dd'))/2+datenum('2016-07-19', 'yyyy-mm-dd'); (datenum('2016-12-06', 'yyyy-mm-dd')-datenum('2016-07-19', 'yyyy-mm-dd'))/2+datenum('2016-12-06', 'yyyy-mm-dd'); (datenum('2017-03-14', 'yyyy-mm-dd')-datenum('2016-12-06', 'yyyy-mm-dd'))/2+datenum('2017-03-14', 'yyyy-mm-dd'); (datenum('2018-02-13', 'yyyy-mm-dd')-datenum('2017-03-14', 'yyyy-mm-dd'))/2+datenum('2018-02-13', 'yyyy-mm-dd'); (datenum('2018-07-03', 'yyyy-mm-dd')-datenum('2018-02-13', 'yyyy-mm-dd'))/2+datenum('2018-07-03', 'yyyy-mm-dd'); (datenum('2018-08-28', 'yyyy-mm-dd')-datenum('2018-07-03', 'yyyy-mm-dd'))/2+datenum('2018-08-28', 'yyyy-mm-dd');]; timeA131dif=timeA131-736436; err_ran=15; %% %Point at a detached block between Futagawa and Idenokuchi faults pa1data=[ A131int20160426_20160510a(p1y, p1x); A131int20160510_20160719a(p1y, p1x); A131int20160719_20161206a(p1y, p1x); A131int20161206_20170314a(p1y, p1x); A131int20170314_20180213a(p1y, p1x); A131int20180213_20180703a(p1y, p1x); A131int20180703_20180828a(p1y, p1x)]; pa1cum=cumsum(pa1data); pa1_err=[ std(std((A131int20160426_20160510a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((A131int20160510_20160719a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((A131int20160719_20161206a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((A131int20161206_20170314a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((A131int20170314_20180213a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((A131int20180213_20180703a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran)))); std(std((A131int20180703_20180828a(p1y-err_ran:p1y+err_ran, p1x-err_ran:p1x+err_ran))))]; ta=timeA131dif;ya1=pa1cum; pa1_fittype1=fittype('a*log(1+(ta/tau))', 'dependent', {'ya1'}, 'independent', {'ta'}, ... 'coefficients', {'a', 'tau'}); pa1_fittype2=fittype('b*ta+c', 'dependent', {'ya1'}, 'independent', {'ta'}, ... 'coefficients', {'b', 'c'}); %startPoints=[-3 0.1 -3 1 0.1];ExcludeP=[100 1000 100 1000]; myfita11=fit(ta, ya1, pa1_fittype1, 'Lower',[-100 0.000001] , 'Upper', [100 1000] ) myfita12=fit(ta, ya1, pa1_fittype2, 'Lower',[-100 -100] , 'Upper', [1000 1000] ) SSresida11=sum((myfita11(ta)-pa1cum).^2); SSresida12=sum((myfita12(ta)-pa1cum).^2); SStotala1=(length(ta)-1)*var(pa1cum); rsqa11=1-SSresida11/SStotala1 rsqa12=1-SSresida12/SStotala1 timea=linspace(1,1000, 1000); fita11=myfita11(timea); fita12=myfita12(timea); timeaP=linspace(736436, 737425,1000); %rsq12=1-SSresid2/SStotal %figure, %plot(myfit11, 'r', t, y1);hold on;xlim([5 1000]); %plot(myfit12, 'b', t, y1); %plot(myfit12,'b', t, y1); % figure, % errorbar(timeA131, pa1cum, pa1_err, 'ko', 'LineWidth', 2); % ylim([-9 1]);grid on;hold on; % plot(timeaP, fita11, 'r');hold on; % scatter(timeA131, pa1cum, 54, 'filled','r', 'MarkerEdgeColor', 'k', 'LineWidth', 2); % plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% %Point at Aso volcano % pa2data=[ % A131int20160426_20160510a(p2y, p2x); % A131int20160510_20160719a(p2y, p2x); % A131int20160719_20161206a(p2y, p2x); % A131int20161206_20170314a(p2y, p2x); % A131int20170314_20180213a(p2y, p2x); % A131int20180213_20180703a(p2y, p2x); % A131int20180703_20180828a(p2y, p2x)]; % pa2cum=cumsum(pa2data); % % pa2_err=[ % std(std((A131int20160426_20160510a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); % std(std((A131int20160510_20160719a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); % std(std((A131int20160719_20161206a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); % std(std((A131int20161206_20170314a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); % std(std((A131int20170314_20180213a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); % std(std((A131int20180213_20180703a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran)))); % std(std((A131int20180703_20180828a(p2y-err_ran:p2y+err_ran, p2x-err_ran:p2x+err_ran))))]; % % ta=timeA131dif;ya2=pa2cum; % pa2_fittype1=fittype('a*log(1+(ta/tau))', 'dependent', {'ya2'}, 'independent', {'ta'}, ... % 'coefficients', {'a', 'tau'}); % %p2_fittype2=fittype('a*log(1+(t/tau))+b*exp(-t/taug)+c', 'dependent', {'y1'}, 'independent', {'t'}, ... % % 'coefficients', {'a', 'tau', 'b','taug', 'c'}); % %startPoints=[-3 0.1 -3 1 0.1];ExcludeP=[100 1000 100 1000]; % % myfita21=fit(ta, ya2, pa2_fittype1, 'Lower',[-30 0.000001] , 'Upper', [100 1000] ) % %myfit12=fit(t, y1, p2_fittype2, 'Lower',[-3 0.1 -3 0.1 -100] , 'Upper', [100 1000 100 1000 100] ) % SSresida2=sum((myfita21(t)-pa2cum).^2); % %SSresid2=sum((myfit12(t)-p2cum).^2); % SStotala=(length(ta)-1)*var(pa2cum); % rsqa21=1-SSresida2/SStotala % timea=linspace(1,1000, 1000); % fita21=myfita11(timea); % timeaP=linspace(736436, 737425,1000); % %rsq12=1-SSresid2/SStotal % %figure, % %plot(myfit11, 'r', t, y1);hold on;xlim([5 1000]); % %plot(myfit12, 'b', t, y1); % % %plot(myfit12,'b', t, y1); % % % figure, % errorbar(timeA131, pa1cum, pa1_err, 'ko', 'LineWidth', 2); % ylim([-9 1]);grid on;hold on; % plot(timeaP, fita11, 'r');hold on; % scatter(timeA131, pa1cum, 54, 'filled','r', 'MarkerEdgeColor', 'k', 'LineWidth', 2); % plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% %Point at north flank %p3x=826; p3y=901; pa3data=[ A131int20160426_20160510a(p3y, p3x); A131int20160510_20160719a(p3y, p3x); A131int20160719_20161206a(p3y, p3x); A131int20161206_20170314a(p3y, p3x); A131int20170314_20180213a(p3y, p3x); A131int20180213_20180703a(p3y, p3x); A131int20180703_20180828a(p3y, p3x)]; pa3cum=cumsum(pa3data); pa3_err=[ std(std((A131int20160426_20160510a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((A131int20160510_20160719a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((A131int20160719_20161206a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((A131int20161206_20170314a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((A131int20170314_20180213a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((A131int20180213_20180703a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran)))); std(std((A131int20180703_20180828a(p3y-err_ran:p3y+err_ran, p3x-err_ran:p3x+err_ran))))]; ta=timeA131dif;ya3=pa3cum; pa3_fittype1=fittype('a*log(1+(ta/tau))', 'dependent', {'ya3'}, 'independent', {'ta'}, ... 'coefficients', {'a', 'tau'}); pa3_fittype2=fittype('b*ta+c', 'dependent', {'ya3'}, 'independent', {'ta'}, ... 'coefficients', { 'b', 'c'}); %startPoints=[-3 0.1 -3 1 0.1];ExcludeP=[100 1000 100 1000]; myfita31=fit(ta, ya3, pa3_fittype1, 'Lower',[-100 0.000001] , 'Upper', [100 1000] ) myfita32=fit(ta, ya3, pa3_fittype2, 'Lower',[-30 -100] , 'Upper', [1000 1000 ] ) SSresida31=sum((myfita31(ta)-pa3cum).^2); SSresida32=sum((myfita32(ta)-pa3cum).^2); SStotala=(length(ta)-1)*var(pa3cum); rsqa31=1-SSresida31/SStotala rsqa32=1-SSresida32/SStotala timea=linspace(1,1000, 1000); fita31=myfita31(timea); fita32=myfita32(timea); timeaP=linspace(736436, 737425,1000); %rsq12=1-SSresid2/SStotal %figure, %plot(myfit11, 'r', t, y1);hold on;xlim([5 1000]); %plot(myfit12, 'b', t, y1); %plot(myfit12,'b', t, y1); % figure, % errorbar(timeA131, pa3cum, pa3_err, 'ko', 'LineWidth', 2); % ylim([-9 1]);grid on;hold on; % plot(timeaP, fita31, 'r');hold on; % scatter(timeA131, pa3cum, 54, 'filled','r', 'MarkerEdgeColor', 'k', 'LineWidth', 2); % plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% %Point at south flank pa4data=[ A131int20160426_20160510a(p4y, p4x); A131int20160510_20160719a(p4y, p4x); A131int20160719_20161206a(p4y, p4x); A131int20161206_20170314a(p4y, p4x); A131int20170314_20180213a(p4y, p4x); A131int20180213_20180703a(p4y, p4x); A131int20180703_20180828a(p4y, p4x)]; pa4cum=cumsum(pa4data); pa4_err=[ std(std((A131int20160426_20160510a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((A131int20160510_20160719a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((A131int20160719_20161206a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((A131int20161206_20170314a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((A131int20170314_20180213a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((A131int20180213_20180703a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran)))); std(std((A131int20180703_20180828a(p4y-err_ran:p4y+err_ran, p4x-err_ran:p4x+err_ran))))]; ta=timeA131dif;ya4=pa4cum; pa4_fittype1=fittype('a*log(1+(ta/tau))', 'dependent', {'ya4'}, 'independent', {'ta'}, ... 'coefficients', {'a', 'tau'}); pa4_fittype2=fittype('b*ta+c', 'dependent', {'ya4'}, 'independent', {'ta'}, ... 'coefficients', {'b', 'c'}); %startPoints=[-3 0.1 -3 1 0.1];ExcludeP=[100 1000 100 1000]; myfita41=fit(ta, ya4, pa4_fittype1, 'Lower',[-100 0.000001] , 'Upper', [100 1000] ) myfita42=fit(ta, ya4, pa4_fittype2, 'Lower',[-30 -100] , 'Upper', [1000 1000] ) SSresida41=sum((myfita41(ta)-pa4cum).^2); SSresida42=sum((myfita42(ta)-pa4cum).^2); SStotala=(length(ta)-1)*var(pa4cum); rsqa41=1-SSresida41/SStotala rsqa42=1-SSresida42/SStotala timea=linspace(1,1000, 1000); fita41=myfita41(timea); fita42=myfita42(timea); timeaP=linspace(736436, 737425,1000); %rsq12=1-SSresid2/SStotal %figure, %plot(myfit11, 'r', t, y1);hold on;xlim([5 1000]); %plot(myfit12, 'b', t, y1); %plot(myfit12,'b', t, y1); % figure, % errorbar(timeA131, pa4cum, pa4_err, 'ko', 'LineWidth', 2); % ylim([-9 1]);grid on;hold on; % plot(timeaP, fita41, 'r');hold on; % scatter(timeA131, pa4cum, 54, 'filled','r', 'MarkerEdgeColor', 'k', 'LineWidth', 2); % plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% %Point at fault junction pa5data=[ A131int20160426_20160510a(p5y, p5x)+2; A131int20160510_20160719a(p5y, p5x); A131int20160719_20161206a(p5y, p5x); A131int20161206_20170314a(p5y, p5x); A131int20170314_20180213a(p5y, p5x); A131int20180213_20180703a(p5y, p5x); A131int20180703_20180828a(p5y, p5x)]; pa5cum=cumsum(pa5data); pa5_err=[ std(std((A131int20160426_20160510a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((A131int20160510_20160719a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((A131int20160719_20161206a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((A131int20161206_20170314a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((A131int20170314_20180213a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((A131int20180213_20180703a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran)))); std(std((A131int20180703_20180828a(p5y-err_ran:p5y+err_ran, p5x-err_ran:p5x+err_ran))))]; ta=timeA131dif;ya5=pa5cum; pa5_fittype1=fittype('a*log(1+(ta/tau))', 'dependent', {'ya5'}, 'independent', {'ta'}, ... 'coefficients', {'a', 'tau'}); %p5_fittype2=fittype('a*log(1+(t/tau))+b*exp(-t/taug)+c', 'dependent', {'y1'}, 'independent', {'t'}, ... % 'coefficients', {'a', 'tau', 'b','taug', 'c'}); %startPoints=[-3 0.1 -3 1 0.1];ExcludeP=[100 1000 100 1000]; myfita51=fit(ta, ya5, pa5_fittype1, 'Lower',[-100 0.000001] , 'Upper', [100 1000] ) %myfit12=fit(t, y1, p5_fittype2, 'Lower',[-3 0.1 -3 0.1 -100] , 'Upper', [100 1000 100 1000 100] ) SSresida5=sum((myfita51(ta)-pa5cum).^2); %SSresid2=sum((myfit12(ta)-p5cum).^2); SStotala=(length(ta)-1)*var(pa5cum); rsqa51=1-SSresida5/SStotala timea=linspace(1,1000, 1000); fita51=myfita51(timea); timeaP=linspace(736436, 737425,1000); %rsq12=1-SSresid2/SStotal %figure, %plot(myfit11, 'r', t, y1);hold on;xlim([5 1000]); %plot(myfit12, 'b', t, y1); %plot(myfit12,'b', t, y1); % figure, % errorbar(timeA131, pa5cum, pa5_err, 'ko', 'LineWidth', 2); % ylim([0 10]);grid on;hold on; % plot(timeaP, fita51, 'r');hold on; % scatter(timeA131, pa5cum, 54, 'filled','r', 'MarkerEdgeColor', 'k', 'LineWidth', 2); % plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% %Point at Kumamoto delta pa6data=[ A131int20160426_20160510a(p6y, p6x); A131int20160510_20160719a(p6y, p6x); A131int20160719_20161206a(p6y, p6x); A131int20161206_20170314a(p6y, p6x); A131int20170314_20180213a(p6y, p6x); A131int20180213_20180703a(p6y, p6x); A131int20180703_20180828a(p6y, p6x)]; pa6cum=cumsum(pa6data); pa6_err=[ std(std((A131int20160426_20160510a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((A131int20160510_20160719a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((A131int20160719_20161206a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((A131int20161206_20170314a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((A131int20170314_20180213a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((A131int20180213_20180703a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran)))); std(std((A131int20180703_20180828a(p6y-err_ran:p6y+err_ran, p6x-err_ran:p6x+err_ran))))]; ta=timeA131dif;ya6=pa6cum; pa6_fittype1=fittype('a*log(1+(ta/tau))', 'dependent', {'ya6'}, 'independent', {'ta'}, ... 'coefficients', {'a', 'tau'}); pa6_fittype2=fittype('b*ta+c', 'dependent', {'ya6'}, 'independent', {'ta'}, ... 'coefficients', {'b', 'c'}); %startPoints=[-3 0.1 -3 1 0.1];ExcludeP=[100 1000 100 1000]; myfita61=fit(ta, ya6, pa6_fittype1, 'Lower',[-100 0.000001] , 'Upper', [100 1000] ) myfita62=fit(ta, ya6, pa6_fittype2, 'Lower',[-100 -100] , 'Upper', [100 100] ) SSresida61=sum((myfita61(ta)-pa6cum).^2); SSresida62=sum((myfita62(ta)-pa6cum).^2); SStotala=(length(ta)-1)*var(pa6cum); rsqa61=1-SSresida61/SStotala rsqa62=1-SSresida62/SStotala timea=linspace(1,1000, 1000); fita61=myfita61(timea); fita62=myfita62(timea); timeaP=linspace(736436, 737425,1000); %rsq12=1-SSresid2/SStotal %figure, %plot(myfit11, 'r', t, y1);hold on;xlim([5 1000]); %plot(myfit12, 'b', t, y1); %plot(myfit12,'b', t, y1); % figure, % errorbar(timeA131, pa6cum, pa6_err, 'ko', 'LineWidth', 2); % ylim([0 10]);grid on;hold on; % plot(timeaP, fita61, 'r');hold on; % plot(timeaP, fita62, 'r');hold on; % scatter(timeA131, pa6cum, 54, 'filled','r', 'MarkerEdgeColor', 'k', 'LineWidth', 2); % plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% %Point at Kumamoto city pa7data=[ A131int20160426_20160510a(p7y, p7x); A131int20160510_20160719a(p7y, p7x); A131int20160719_20161206a(p7y, p7x); A131int20161206_20170314a(p7y, p7x); A131int20170314_20180213a(p7y, p7x); A131int20180213_20180703a(p7y, p7x); A131int20180703_20180828a(p7y, p7x)]; pa7cum=cumsum(pa7data); pa7_err=[ std(std((A131int20160426_20160510a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((A131int20160510_20160719a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((A131int20160719_20161206a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((A131int20161206_20170314a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((A131int20170314_20180213a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((A131int20180213_20180703a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran)))); std(std((A131int20180703_20180828a(p7y-err_ran:p7y+err_ran, p7x-err_ran:p7x+err_ran))))]; ta=timeA131dif;ya7=pa7cum; pa7_fittype1=fittype('a*log(1+(ta/tau))', 'dependent', {'ya7'}, 'independent', {'ta'}, ... 'coefficients', {'a', 'tau'}); pa7_fittype2=fittype('b*ta+c', 'dependent', {'ya7'}, 'independent', {'ta'}, ... 'coefficients', {'b', 'c'}); myfita71=fit(ta, ya7, pa7_fittype1, 'Lower',[-100 0.000001] , 'Upper', [100 1000] ) myfita72=fit(ta, ya7, pa7_fittype2, 'Lower',[-10 -10 ] , 'Upper', [10 10 ] ) SSresida71=sum((myfita71(ta)-pa7cum).^2); SSresida72=sum((myfita72(ta)-pa7cum).^2); SStotala=(length(ta)-1)*var(pa7cum); rsqa71=1-SSresida71/SStotala rsqa72=1-SSresida72/SStotala timea=linspace(1,1000, 1000); fita71=myfita71(timea); fita72=myfita72(timea); timeaP=linspace(736436, 737425,1000); %rsq12=1-SSresid2/SStotal %figure, %plot(myfit11, 'r', t, y1);hold on;xlim([5 1000]); %plot(myfit12, 'b', t, y1); %plot(myfit12,'b', t, y1); % figure, % errorbar(timeA131, pa7cum, pa7_err, 'ko', 'LineWidth', 2); % ylim([-1 10]);grid on;hold on; % plot(timeaP, fita71, 'r');hold on; % plot(timeaP, fita72, 'r');hold on; % scatter(timeA131, pa7cum, 54, 'filled','r', 'MarkerEdgeColor', 'k', 'LineWidth', 2); % plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% pa8data=[ A131int20160426_20160510a(p8y, p8x); A131int20160510_20160719a(p8y, p8x); A131int20160719_20161206a(p8y, p8x); A131int20161206_20170314a(p8y, p8x); A131int20170314_20180213a(p8y, p8x); A131int20180213_20180703a(p8y, p8x); A131int20180703_20180828a(p8y, p8x)]; pa8cum=cumsum(pa8data); pa8_err=[ std(std((A131int20160426_20160510a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((A131int20160510_20160719a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((A131int20160719_20161206a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((A131int20161206_20170314a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((A131int20170314_20180213a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((A131int20180213_20180703a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran)))); std(std((A131int20180703_20180828a(p8y-err_ran:p8y+err_ran, p8x-err_ran:p8x+err_ran))))]; ta=timeA131dif;ya8=pa8cum; pa8_fittype1=fittype('a*log(1+(ta/tau))', 'dependent', {'ya8'}, 'independent', {'ta'}, ... 'coefficients', {'a', 'tau'}); pa8_fittype2=fittype('b*ta+c', 'dependent', {'ya8'}, 'independent', {'ta'}, ... 'coefficients', {'b', 'c'}); myfita81=fit(ta, ya8, pa8_fittype1, 'Lower',[-100 0.000001] , 'Upper', [100 1000] ) myfita82=fit(ta, ya8, pa8_fittype2, 'Lower',[-10 -10 ] , 'Upper', [10 10 ] ) SSresida81=sum((myfita81(ta)-pa8cum).^2); SSresida82=sum((myfita82(ta)-pa8cum).^2); SStotala=(length(ta)-1)*var(pa8cum); rsqa81=1-SSresida81/SStotala rsqa82=1-SSresida82/SStotala timea=linspace(1,1000, 1000); fita81=myfita81(timea); fita82=myfita82(timea); timeaP=linspace(736436, 737425,1000); %rsq12=1-SSresid2/SStotal %figure, %plot(myfit11, 'r', t, y1);hold on;xlim([5 1000]); %plot(myfit12, 'b', t, y1); %plot(myfit12,'b', t, y1); % figure, % errorbar(timeA131, pa8cum, pa8_err, 'ko', 'LineWidth', 2); % ylim([-1 10]);grid on;hold on; % plot(timeaP, fita81, 'r');hold on; % plot(timeaP, fita82, 'r');hold on; % scatter(timeA131, pa8cum, 54, 'filled','r', 'MarkerEdgeColor', 'k', 'LineWidth', 2); % plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 14); %% figure, subplot(4,2,1); errorbar(timeA131, pa1cum, pa1_err, 'ko', 'LineWidth', 2);ylim([-11 3]);grid on;hold on; plot(timeaP, fita11, 'r', 'LineWidth', 3);hold on; plot(timeaP, fita12, 'b', 'LineWidth', 1);hold on; scatter(timeA131, pa1cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); % subplot(3,1,2); % errorbar(timeA131, p2cum, p2_err, 'ko', 'LineWidth', 2);ylim([0 10]);grid on;hold on; % plot(timeP, fit21, 'k');hold on; % scatter(timeA131, p2cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); % plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); % xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); % xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); % ylabel('LOS change [cm]'); % set(gca, 'FontSize', 12); subplot(4,2,5); errorbar(timeA131, pa3cum, pa3_err, 'ko', 'LineWidth', 2);ylim([-11 3]);grid on;hold on; plot(timeaP, fita31, 'r', 'LineWidth', 3);hold on; plot(timeaP, fita32, 'b', 'LineWidth', 1);hold on; scatter(timeA131, pa3cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); subplot(4,2,7); errorbar(timeA131, pa4cum, pa4_err, 'ko', 'LineWidth', 2);ylim([-11 3]);grid on;hold on; plot(timeaP, fita41, 'r', 'LineWidth', 3);hold on; plot(timeaP, fita42, 'b', 'LineWidth', 1);hold on; scatter(timeA131, pa4cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); subplot(4,2,2); errorbar(timeA131, pa5cum, pa5_err, 'ko', 'LineWidth', 2);ylim([-3 11]);grid on;hold on; plot(timeaP, fita51, 'r', 'LineWidth', 3);hold on; scatter(timeA131, pa5cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); subplot(4,2,4); errorbar(timeA131, pa6cum, pa6_err, 'ko', 'LineWidth', 2);ylim([-3 11]);grid on;hold on; plot(timeaP, fita61, 'r', 'LineWidth', 3);hold on; plot(timeaP, fita62, 'b', 'LineWidth', 1);hold on; scatter(timeA131, pa6cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); subplot(4,2,6); errorbar(timeA131, pa7cum, pa7_err, 'ko', 'LineWidth', 2);ylim([-3 11]);grid on;hold on; plot(timeaP, fita71, 'r', 'LineWidth', 3);hold on; plot(timeaP, fita72, 'b', 'LineWidth', 1);hold on; scatter(timeA131, pa7cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); subplot(4,2,8); errorbar(timeA131, pa8cum, pa8_err, 'ko', 'LineWidth', 2);ylim([-3 11]);grid on;hold on; plot(timeaP, fita81, 'r', 'LineWidth', 3);hold on; plot(timeaP, fita82, 'b', 'LineWidth', 1);hold on; scatter(timeA131, pa8cum, 54, 'k', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 2); plot([736436 736436], [-20 20], 'k', 'LineWidth', 3); xticks([736512, 736696, 736877, 737061, 737242]);xlim([736390 737425]); xticklabels({'Jul. 2016', 'Jan. 2017', 'Jul. 2017', 'Jan. 2018', 'Jul. 2018'}); ylabel('LOS change [cm]'); set(gca, 'FontSize', 12); %% rsqa=[rsqa11; NaN; rsqa31; rsqa41; rsqa51; rsqa61; rsqa62; rsqa71; rsqa72]; rsqatag=['rsqa11='; 'rsqa21='; 'rsqa31='; 'rsqa41='; 'rsqa51='; 'rsqa61='; 'rsqa62='; 'rsqa71='; 'rsqa72=']; disp([rsqatag num2str(rsqa)]); disp(myfita11); %disp(myfita21); disp(myfita31); disp(myfita41); disp(myfita51); disp(myfita61); disp(myfita62); disp(myfita71); disp(myfita72); [tv,pva11]=ttest(myfita11(ta)-pa1cum) [tv,pva12]=ttest(myfita12(ta)-pa1cum) %[tv,pva21]=ttest(myfita21(ta)-pa2cum) [tv,pva31]=ttest(myfita31(ta)-pa3cum) [tv,pva41]=ttest(myfita41(ta)-pa4cum) [tv,pva51]=ttest(myfita51(ta)-pa5cum) [tv,pva61]=ttest(myfita61(ta)-pa6cum) [tv,pva62]=ttest(myfita62(ta)-pa6cum) [tv,pva71]=ttest(myfita71(ta)-pa7cum) [tv,pva72]=ttest(myfita72(ta)-pa7cum) [tv,pva81]=ttest(myfita81(ta)-pa8cum) [tv,pva82]=ttest(myfita82(ta)-pa8cum) disp(tv); %%EOF