%% best fit position and slip- "look ahead" and "subevent" version
% Comparing observed pressure values with Pylith model to find best fit
% positions and slip amounts

% last updated 9/9/2024
% Joshua Edgington 



%widths=[5,10,20,40]; %assign all 4 patch widths
widths=[20]; %assign specific patch widths

event=2015; %event year based on loaded filenames
maxslip=2; %maximum amount of allowed slip per 2-day window (in cm)
time_delP=eq_timeC2; %assign single time vecto

for yyy=1:4
%% load pylith green's functions, and slip data
sizep=widths(yyy);

if sizep==5
py_C2=P5km.C2kPa;
py_C10=P5km.C10_1kPa;
py_C6=P5km.C6kPa;   
X_Center=P5km.X_Center;

elseif sizep==10
py_C2=P10km.C2kPa;
py_C10=P10km.C10_1kPa;
py_C6=P10km.C6kPa;
X_Center=P10km.X_Center;

elseif sizep==20
py_C2=P20km.C2kPa;
py_C10=P20km.C10_1kPa;
py_C6=P20km.C6kPa;
X_Center=P20km.X_Center;

elseif sizep==40
py_C2=P40km.C2kPa;
py_C10=P40km.C10_1kPa;
py_C6=P40km.C6kPa;
X_Center=P40km.X_Center;
end

% differences (for array length defining purposes)
diffs=cat(2,eq_dataC2, eq_dataC10, eq_dataC6);
diffs=diffs(1:end,:);

%% normalize green's functions 
diffs_py=cat(2,py_C2, py_C10, py_C6);

%% extract data values
p_dat=diffs;
time=eq_timeC2;
%% look ahead loop variables
for aa=1
% preallocation
x_pos=X_Center;

best_fit2=zeros(ww,1);

err_fit=zeros(ww,1);

% define model data
model_dat=zeros(ww,3);

model_dat(1,:)=p_dat(1,:); % set first point as data point

model_cumi=p_dat(1,:);

model_cum=zeros(ww,3);

model_cum(1,:)=p_dat(1,:);

grid_errors=zeros(ww,3);
% define slip amount and error
slip=zeros(ww,1);

err_iter=zeros(ww,1);

slipit=(0.00:0.1:maxslip); % slip amounts for grid search

slipit=slipit';

custompos=123;
end

% if piecewise
%ww=length(eq_timeC2);
%% 
%%%%%%%%% THE LOOP %%%%%%%%%
for ii=1:ww-1
    
    %%%%%%%%%%%%%%%%%%
    %%% set deltaP %%%
    %%%%%%%%%%%%%%%%%%

    delP=diffs(ii,:); % select sub event
    
    % preallocate arrays
    err_py=zeros(length(diffs_py(:,3)),1);

    slip_py=zeros(length(diffs_py(:,3)),1);

    %%%%%%%%%%%%%%%
    %%% No delP %%%
    %%%%%%%%%%%%%%%

     if sum(delP)==0 

        slip(ii)=nan;
        
        best_fit2(ii)=NaN;
        
        model_dat(ii+1,:)=model_dat(ii,:); 
       
        model_cumi=model_cumi+[0,0,0];

        model_cum(ii+1,:)=model_cumi;

     else

    %%%%%%%%%%%%%%%%%%%%
    %%% thresholding %%%
    %%%%%%%%%%%%%%%%%%%%
    
   mdp=0.01; % minimum delP before setting to zero

   % C2
   if abs(delP(:,1))<mdp
       delP(:,1)=0;
   end

   % C10
   if abs(delP(:,2))<mdp
    delP(:,2)=0;
   end

   % C6
   if abs(delP(:,3))<mdp
    delP(:,3)=0;
   end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% best fit position loop %%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        err_all=zeros(length(slipit),length(diffs_py(:,3))); % store error in grid
        err_mins=zeros(55,1); % store lowest errors for each position/GF

        all_grids=zeros(21,3,55);

        for yy=1:length(diffs_py(:,3)) %for each slip position, one set of GFs
            
           % look across all positions' GF extract GFs
           diffa=diffs_py(yy,:);
            
           % build grid
            p_mod=zeros(length(slipit),3);

                for xx=1:length(slipit) 

                    p_mod(xx,:)=slipit(xx)*diffa; %generate model delP

                end

            if event==2016 %turn off C10 and C6 if not present
                p_mod(:,3)=nan;
                p_mod(:,2)=nan;
                delP(:,2)=nan;
                delP(:,3)=nan;
            end   

            if event==2015 %turn off C6 if not present
                p_mod(:,3)=nan;
                delP(:,3)=nan;
            end
            
            %%% RMSE error scaled by 2 sigma %%%

            errpre=p_mod-delP; %calculate error between observed delP and model del P
            all_grids(:,:,yy)=errpre; %store error grid for lookup


            % scale errors by inverse of noise (by 2-sigma)
            if event==2015
                errpre2=errpre.*[(1/C2_2sig),(1/C10_2sig),0];
            elseif event==2016
                errpre2=errpre.*[(1/C2_2sig),0,0];
            else
                errpre2=errpre.*[1/C2_2sig,(1/C10_2sig),1/C6_2sig];  
            end

            %all_grids(:,:,yy)=errpre2; %store error grid for lookup

            % RMSE
            err = sqrt(mean((errpre2).^2,2,'omitnan'));
           
            err_all(:,yy)=err;

            [~,idx]=min(err);

            err_py(yy)=min(err);

            slip_py(yy)=slipit(idx); %store in array 

            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
        end 

        %%%%%%%%%%%%%%%%%%%%%
        %%% assign values %%%
        %%%%%%%%%%%%%%%%%%%%%

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% MODEL POSITION AND SLIP %%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            % extract exact errors
            
            [~,idx]=min(err_py);

            [~,gridpage]=min(err_py); %get page of lowest error

            [row_err,~]=find(slipit==slip_py(idx)); %row of errors from page
            
            grid_errors(ii,:)=all_grids(row_err,:,gridpage); %extract exact errors and store

            
            % define best fits

            best_fit2(ii)=x_pos(idx);
            
            err_iter(ii)=err_py(idx);
    
            slip(ii)=slip_py(idx); %store in array 
    
            model_dat(ii,:)=slip(ii)*diffs_py(idx,:);  
    
            model_cumi=model_cumi+(slip(ii)*diffs_py(idx,:));
    
            model_cum(ii+1,:)=model_cumi;

     
     end
%generate error plot

        if ii==30
                figure
                h=surf(X_Center./1000,slipit,err_all)
                hold on

                mini = min(min(err_all));
                [x,y]=find(err_all==mini)


                plot3(X_Center(y)/1000,slipit(x),min(min(err_all)),'r*')
                set(h,'LineStyle','none')
                xlabel('Centroid Location');
                ylabel('Slip');
                zlabel('Error')
                title('19 Mar 2020 [5km patch]');

        end
end



best_fit2(end)=nan;
slip(end)=nan;
err_iter(end)=nan;

%% compare delta P
for aa=1
firstin = find(~isnan(best_fit2), 1);
lastin = find(~isnan(best_fit2), 1, 'last');

mod_shift_C2=model_dat(:,1);
mod_shift_C10=model_dat(:,2);
mod_shift_C6=model_dat(:,3);

mod_shift_C2(1:firstin-1)=mod_shift_C2(1:firstin-1)-model_dat(1,1);
mod_shift_C10(1:firstin-1)=mod_shift_C10(1:firstin-1)-model_dat(1,2);
mod_shift_C6(1:firstin-1)=mod_shift_C6(1:firstin-1)-model_dat(1,3);

diffs_errplot=diffs;

figure
cx(1)=subplot(3,1,1);
plot(time_delP,diffs_errplot(:,1),'k_','linewidth',3);
hold on

if event==2015 || event==2016
    err5=C2_2sig+zeros(size(mod_shift_C2));
    errorbar(time_delP,diffs(:,1),err5,'ks','LineWidth', 1)
%plot(time,diffs_errplot(:,1)+0.039067044194372,'r--','linewidth',2);
%plot(time,diffs_errplot(:,1)-0.039067044194372,'r--','linewidth',2);
else
    err5=C2_2sig+zeros(size(mod_shift_C2));
    errorbar(time_delP,diffs(:,1),err5,'ks','LineWidth', 1)
%plot(time,diffs_errplot(:,1)+0.035050569836934,'r--','linewidth',2);
%plot(time,diffs_errplot(:,1)-0.035050569836934,'r--','linewidth',2); 
end

plot(time_delP,mod_shift_C2,'b','linewidth',2); 

ytickformat('%.2f');
datetick('x','keeplimits');
axis 'tight';
xlabel('Time');
ylabel('Pressure (KPa)');
title('delta P C2');
legend('data','2sigma','model')
%xlim([time(1) time(end)])

cx(2)=subplot(3,1,2);
plot(time_delP,diffs_errplot(:,2),'k_','linewidth',3);
hold on

if event==2015 || event==2016

    err5=C10_2sig+zeros(size(mod_shift_C10));
    errorbar(time_delP,diffs(:,2),err5,'ks','LineWidth', 1)
%plot(time,diffs_errplot(:,2)+0.045691413022237,'r_','linewidth',2);
%plot(time,diffs_errplot(:,2)-0.045691413022237,'r_','linewidth',2);
else
    err5=C10_2sig+zeros(size(mod_shift_C10));
    errorbar(time_delP,diffs_errplot(:,2),err5,'ks','LineWidth', 1)
%plot(time,diffs_errplot(:,2)+0.174477359714062,'r_','linewidth',2);
%plot(time,diffs_errplot(:,2)-0.174477359714062,'r_','linewidth',2);
end

plot(time_delP,mod_shift_C10,'b','linewidth',2); 

ytickformat('%.2f');
datetick('x','keeplimits');
axis 'tight';
xlabel('Time');
ylabel('Pressure (KPa)');
title('delta P C10');
legend('data','2sigma','model')
%xlim([time(1) time(end)])

cx(3)=subplot(3,1,3);
plot(time_delP,diffs_errplot(:,3),'k_','linewidth',3);
hold on


if event==2015 || event==2016
%'meme'
else
    err5=C6_2sig+zeros(size(mod_shift_C6));
    errorbar(time_delP,diffs_errplot(:,3),err5,'ks','LineWidth', 1)

%plot(time,diffs_errplot(:,3)+0.045063135745060,'r--','linewidth',2);
%plot(time,diffs_errplot(:,3)-0.045063135745060,'r--','linewidth',2);
end

plot(time_delP,mod_shift_C6,'b','linewidth',2); 
ytickformat('%.2f');
datetick('x','keeplimits');
axis 'tight';
xlabel('Time');
ylabel('Pressure (KPa)');
title('delta P C6');
legend('data','2sigma','model')
%xlim([time(1) time(end)])
end
%% variable assignment
for aa=1
%err_iter(err_iter==0)=NaN;

if sizep==5
 cumpC2_5km=model_cum(:,1);
 cumpC10_5km=model_cum(:,2);
 cumpC6_5km=model_cum(:,3);
 err_5km=err_iter;
 slip_5km=slip;
 best_fit5km=best_fit2;
end

if sizep==10
 cumpC2_10km=model_cum(:,1);
 cumpC10_10km=model_cum(:,2);
 cumpC6_10km=model_cum(:,3);
 err_10km=err_iter;
 slip_10km=slip;
 best_fit10km=best_fit2;
end

if sizep==20
 cumpC2_20km=model_cum(:,1);
 cumpC10_20km=model_cum(:,2);
 cumpC6_20km=model_cum(:,3);
 err_20km=err_iter;
 slip_20km=slip;
 best_fit20km=best_fit2;
end

if sizep==40
 cumpC2_40km=model_cum(:,1);
 cumpC10_40km=model_cum(:,2);
 cumpC6_40km=model_cum(:,3);
 err_40km=err_iter;
 slip_40km=slip;
 best_fit40km=best_fit2;
end


end
%% plotting results %%
for aa=1

%%%%%%%%%%%%%%%%%%%%% ONLY C2/10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if event==2015 || event==2016
    for ss=1
   
figure
cx(1)=subplot(4,1,1);
%plot(C2time,C2d,'k','linewidth',1);
hold on
plot(C10time,fits_C2,'b','linewidth',2);

if sizep==5
    plot(time,cumpC2_5km,'r:','linewidth',2);
elseif sizep==10
    plot(time,cumpC2_10km,'r:','linewidth',2);
elseif sizep==20
    plot(time,cumpC2_20km,'r:','linewidth',2);
elseif sizep==40
    plot(time,cumpC2_40km,'r:','linewidth',2);
end
    
ytickformat('%.2f');
datetick('x','keeplimits');
%axis 'tight';
%ylim([-2 2])
xlabel('Time');
ylabel('Pressure (KPa)');
title('C2');
legend('Data','Smooth','Result')
xlim([time(1) time(end)])
%xlim([time(1) time(end)])

cx(2)=subplot(4,1,2);
%plot(C10time,C10d,'k','linewidth',1);
hold on
plot(C10time,fits_C10,'b','linewidth',2);


if sizep==5
    plot(time,cumpC10_5km,'r:','linewidth',2);
elseif sizep==10
    plot(time,cumpC10_10km,'r:','linewidth',2);
elseif sizep==20
    plot(time,cumpC10_20km,'r:','linewidth',2);
elseif sizep==40
    plot(time,cumpC10_40km,'r:','linewidth',2);
end

ytickformat('%.2f');
datetick('x','keeplimits');
%axis 'tight';
%ylim([-2 2])
xlabel('Time');
ylabel('Pressure (KPa)');
title('C10');
xlim([time(1) time(end)])


% cx(3)=subplot(3,1,3);
% scatter(time(firstin:lastin),best_fit2(firstin:lastin)./1000,300,'sk','filled');
% hold on
% hline(-2, 'r-', 'C6');
%  hline(-34, 'r-', 'C2');
%  hline(-23,'r-', 'C10');
% hold off
% ylabel('Position (km)');
% set(gca,'YDir','reverse')
% hold on
% yyaxis right
% scatter(time(firstin:lastin),slip(firstin:lastin),300,'dr','filled');
% dx=gca;
% dx.YColor = [0, 0, 0]; % X labels are red.
% ylabel('Slip (cm)');
% ytickformat('%.2f');
% datetick('x','keeplimits');
% 
% xlabel('Time');
% 
% title('Position & Slip');
% 
% 
% linkaxes(cx,'x');
%legend('Input','Output')

cx(3)=subplot(4,1,3);
pointsize=100;

if sizep==5
    scatter(time_delP, best_fit5km, pointsize, slip_5km,'s','filled');
elseif sizep==10
    scatter(time_delP, best_fit10km, pointsize, slip_10km,'s','filled');
elseif sizep==20
    scatter(time_delP, best_fit20km, pointsize, slip_20km,'s','filled');
elseif sizep==40
    scatter(time_delP, best_fit40km, pointsize, slip_40km,'s','filled');
end

     hline(-34000, 'r-', 'C2');
     hline(-23000,'r-', 'C10');
     hline(-2000, 'r-', 'C6');
     
datetick('x','keeplimits');
ax = gca;
ax.YRuler.Exponent = 0;
set(gca,'YDir','reverse')
colorbar
%axis 'tight';
xlabel('Time');
ylabel('KM');
title('Center of Slip Patch (km from trench)');
xlim([time(1) time(end)])
% 
cx(4)=subplot(4,1,4);
hold on
plot(time_delP,slip,'m_','linewidth',2)
datetick('x','keeplimits');
ax = gca;
ax.YRuler.Exponent = 0;
xlabel('Time');
ylabel('cm');
title('slip');
xlim([time(1) time(end)])

linkaxes(cx,'x');

end
end


%%%%%%%%%%%%%%%%%%%%% ALL THREE SITES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if event==2020
    for ss=1
  
figure
cx(1)=subplot(5,1,1);
%plot(C2time,C2d,'k','linewidth',1);
hold on
plot(C10time,fits_C2,'b','linewidth',2);

if sizep==5
    plot(time,cumpC2_5km,'r:','linewidth',2);
elseif sizep==10
    plot(time,cumpC2_10km,'r:','linewidth',2);
elseif sizep==20
    plot(time,cumpC2_20km,'r:','linewidth',2);
elseif sizep==40
    plot(time,cumpC2_40km,'r:','linewidth',2);
end
    
ytickformat('%.2f');
datetick('x','keeplimits');
%axis 'tight';
%ylim([-3 0.5])
xlabel('Time');
ylabel('Pressure (KPa)');
title('C2');
legend('Data','Smooth','Result')
xlim([time(1) time(end)])
%xlim([time(1) time(end)])

cx(2)=subplot(5,1,2);
%plot(C10time,C10d,'k','linewidth',1);
hold on
plot(C10time,fits_C10,'b','linewidth',2);


if sizep==5
    plot(time,cumpC10_5km,'r:','linewidth',2);
elseif sizep==10
    plot(time,cumpC10_10km,'r:','linewidth',2);
elseif sizep==20
    plot(time,cumpC10_20km,'r:','linewidth',2);
elseif sizep==40
    plot(time,cumpC10_40km,'r:','linewidth',2);
end

ytickformat('%.2f');
datetick('x','keeplimits');
%axis 'tight';
%ylim([-3 0.5])
xlabel('Time');
ylabel('Pressure (KPa)');
title('C10');
xlim([time(1) time(end)])
%ylim([-0.5 0.5])
%legend('Input','Output')

cx(3)=subplot(5,1,3);
%plot(C6time,C6d,'k','linewidth',1);
hold on
plot(C10time,fits_C6,'b','linewidth',2);


if sizep==5
    plot(time,cumpC6_5km,'r:','linewidth',2);
elseif sizep==10
    plot(time,cumpC6_10km,'r:','linewidth',2);
elseif sizep==20
    plot(time,cumpC6_20km,'r:','linewidth',2);
elseif sizep==40
    plot(time,cumpC6_40km,'r:','linewidth',2);
end

ytickformat('%.2f');
datetick('x','keeplimits');
%axis 'tight';
%ylim([-3 0.5])
xlabel('Time');
ylabel('Pressure (KPa)');
title('C6');
xlim([time(1) time(end)])
%ylim([-0.15 0.15])
%legend('Input','Output')

% cx(4)=subplot(4,1,4);
% scatter(time(firstin:lastin),best_fit2(firstin:lastin)./1000,300,'sb','filled');
% hold on
% hline(-2, 'r-', 'C6');
%  hline(-34, 'r-', 'C2');
%  hline(-23,'r-', 'C10');
% hold off
% ylabel('Position (km)');
% set(gca,'YDir','reverse')
% hold on
% yyaxis right
% scatter(time(firstin:lastin),slip(firstin:lastin),300,'dk','filled');
% dx=gca;
% dx.YColor = [0, 0, 0]; % X labels are red.
% ylabel('Slip (cm)');
% ytickformat('%.2f');
% datetick('x','keeplimits');
% 
% xlabel('Time');
% 
% title('Position & Slip');
% 
% linkaxes(cx,'x');

cx(4)=subplot(5,1,4);
pointsize=100;

if sizep==5
    scatter(time_delP, best_fit5km, pointsize, slip_5km,'s','filled');
elseif sizep==10
    scatter(time_delP, best_fit10km, pointsize, slip_10km,'s','filled');
elseif sizep==20
    scatter(time_delP, best_fit20km, pointsize, slip_20km,'s','filled');
elseif sizep==40
    scatter(time_delP, best_fit40km, pointsize, slip_40km,'s','filled');
end

     hline(-34000, 'r-', 'C2');
     hline(-23000,'r-', 'C10');
     hline(-2000, 'r-', 'C6');
     
datetick('x','keeplimits');
ax = gca;
ax.YRuler.Exponent = 0;
set(gca,'YDir','reverse')
colorbar
%axis 'tight';
xlabel('Time');
ylabel('KM');
title('Center of Slip Patch (km from trench)');
xlim([time(1) time(end)])

cx(5)=subplot(5,1,5);
hold on
plot(time_delP,slip,'m_','linewidth',2)
datetick('x','keeplimits');
ax = gca;
ax.YRuler.Exponent = 0;
xlabel('Time');
ylabel('cm');
title('slip');
xlim([time(1) time(end)])
linkaxes(cx,'x');
end
end

end
%close all
clc
end
%% dual axis slip-position plot

% firstin = find(~isnan(best_fit2), 1);
% lastin = find(~isnan(best_fit2), 1, 'last');
% 
% figure
% scatter(time(firstin:lastin),best_fit2(firstin:lastin)./1000,300,'sk','filled');
% hold on
% hline(-2, 'r-', 'C6');
%  hline(-34, 'r-', 'C2');
%  hline(-23,'r-', 'C10');
% hold off
% ylabel('Position (km)');
% ax = gca;
% ax.YRuler.Exponent = 0;
% set(gca,'YDir','reverse')
% hold on
% yyaxis right
% scatter(time(firstin:lastin),slip(firstin:lastin),300,'dr','filled');
% ylabel('Slip (cm)');
% ytickformat('%.2f');
% datetick('x','keeplimits');
% 
% xlabel('Time');
% 
% title('Position & Slip');


%% all four inversion results

% ymin2=0;
% ymax2=1.4;
% 
% figure
% 
% cx(1)=subplot(4,1,1)
% firstin = find(~isnan(best_fit5km), 1);
% lastin = find(~isnan(best_fit5km), 1, 'last');
% scatter(time(firstin:lastin),best_fit5km(firstin:lastin)./1000,300,'sk','filled');
% ylabel('Position (km)');
% ax = gca;
% ax.YRuler.Exponent = 0;
% set(gca,'YDir','reverse')
% hold on
% yyaxis right
% scatter(time(firstin:lastin),slip_5km(firstin:lastin),300,'dr','filled');
% dx=gca;
% dx.YColor = [0, 0, 0]; % X labels are red.
% dx.FontSize=12
% ylim([ymin2,ymax2])
% ylabel('Slip (cm)');
% ytickformat('%.2f');
% datetick('x','keeplimits');
% 
% title('5km patch');
% 
% 
% linkaxes(cx,'x');

%% display results
% 
% totslip=nansum(slip)
% peakslip=max(slip)
% 
% if event==2015
% toterr=nansum(err_iter(19:31))
% mederr=median(err_iter(19:31))
% meanerr=mean(err_iter(19:31))
% 
% C2RMSE=sqrt(sum((nonzeros(mod_shift_C2)-nonzeros(diffs(:,1))).^2)/length(nonzeros(diffs(:,1))));
% C10RMSE=sqrt(sum((nonzeros(mod_shift_C10)-nonzeros(diffs(:,2))).^2)/length(nonzeros(diffs(:,2))));
% 
% RMSE=C2RMSE+C10RMSE
% else
% toterr=nansum(err_iter(29:38))
% mederr=median(err_iter(29:38))
% meanerr=mean(err_iter(29:38))
% 
% C2RMSE=sqrt(sum((nonzeros(mod_shift_C2)-nonzeros(diffs(:,1))).^2)/length(nonzeros(diffs(:,1))));
% C10RMSE=sqrt(sum((nonzeros(mod_shift_C10)-nonzeros(diffs(:,2))).^2)/length(nonzeros(diffs(:,2))));
% C6RMSE=sqrt(sum((nonzeros(mod_shift_C6)-nonzeros(diffs(:,3))).^2)/length(nonzeros(diffs(:,3))));
% 
% RMSE=C2RMSE+C10RMSE+C6RMSE
% end


%% layered plots 2015
if event==2015

figure
cx(1)=subplot(4,1,1);
%plot(C10time,fits_C2-fits_C2(1),'b','linewidth',2);
plot(C2time,C2d,'k','linewidth',1);
hold on

    plot(time,cumpC2_5km,'linewidth',2);

    plot(time,cumpC2_10km,'linewidth',2);

    plot(time,cumpC2_20km,'linewidth',2);

    plot(time,cumpC2_40km,'linewidth',2);
    
ytickformat('%.2f');
datetick('x','keeplimits');
xlabel('Time');
ylabel('Pressure (KPa)');
title('C2');
xlim([time(1) time(end)])
ylim([-1.6 0.4])
legend('data','5km','10km','20km','40km')

cx(2)=subplot(4,1,2);
%plot(C10time,fits_C10-fits_C10(1),'b','linewidth',2);
plot(C10time,C10d,'k','linewidth',1);
hold on

    plot(time,cumpC10_5km,'linewidth',2);

    plot(time,cumpC10_10km,'linewidth',2);

    plot(time,cumpC10_20km,'linewidth',2);

    plot(time,cumpC10_40km,'linewidth',2);

ytickformat('%.2f');
datetick('x','keeplimits');
xlabel('Time');
ylabel('Pressure (KPa)');
title('C10');
xlim([time(1) time(end)])
ylim([-0.2 0.8])
%legend('data','5km','10km','20km','40km')

cx(3)=subplot(4,1,3);
pointsize=60;
pointsize2=5;
    scatter(time_delP, best_fit5km/1000, pointsize, 's','filled');
    hold on
    scatter(time_delP, best_fit10km/1000, pointsize, 's','filled');
    scatter(time_delP, best_fit20km/1000, pointsize, 's','filled');
    scatter(time_delP, best_fit40km/1000, pointsize, 's','filled');
    scatter(trem2015.time,trem2015.dist, pointsize2, 'ok', 'filled')

%      hline(-34000, 'r-', 'C2');
%      hline(-23000,'r-', 'C10');
%      hline(-2000, 'r-', 'C6');
     
datetick('x','keeplimits');
ax = gca;
ax.YRuler.Exponent = 0;
set(gca,'YDir','reverse')
xlabel('Time');
ylabel('KM');
title('Center of Slip Patch (km from trench)');
xlim([time(1) time(end)])
ylim([-45 5])
%legend('5km','10km','20km','40km')

cx(4)=subplot(4,1,4);
scatter(time_delP,slip_5km, pointsize, 's','filled')
hold on
scatter(time_delP,slip_10km, pointsize, 's','filled')
scatter(time_delP,slip_20km, pointsize, 's','filled')
scatter(time_delP,slip_40km, pointsize, 's','filled')
datetick('x','keeplimits');
ax = gca;
ax.YRuler.Exponent = 0;
xlabel('Time');
ylabel('cm');
title('slip');
xlim([time(1) time(end)])
ylim([-0.3 2.5])
linkaxes(cx,'x');
legend('5km','10km','20km','40km')

xlim([datenum('27-sep-2015') datenum('17-NOV-2015')])

elseif event==2020

figure
cx(1)=subplot(5,1,1);
plot(C2time,C2d,'k');
%plot(C10time,fits_C2-fits_C2(1),'b','linewidth',2);
hold on

    plot(time,cumpC2_5km,'linewidth',2);

    plot(time,cumpC2_10km,'linewidth',2);

    plot(time,cumpC2_20km,'linewidth',2);

    plot(time,cumpC2_40km,'linewidth',2);
    
ytickformat('%.2f');
datetick('x','keeplimits');
xlabel('Time');
ylabel('Pressure (KPa)');
title('C2');
%legend('Data','Smooth','Result')
xlim([time(1) time(end)])
ylim([-1.6 0.4])
legend('data','5km','10km','20km','40km')

cx(2)=subplot(5,1,2);
%plot(C10time,fits_C10-fits_C10(1),'b','linewidth',2);
plot(C10time,C10d,'k');
hold on

    plot(time,cumpC10_5km,'linewidth',2);

    plot(time,cumpC10_10km,'linewidth',2);

    plot(time,cumpC10_20km,'linewidth',2);

    plot(time,cumpC10_40km,'linewidth',2);

ytickformat('%.2f');
datetick('x','keeplimits');
xlabel('Time');
ylabel('Pressure (KPa)');
title('C10');
xlim([time(1) time(end)])
ylim([-0.2 0.8])
%legend('data','5km','10km','20km','40km')

cx(3)=subplot(5,1,3);
plot(C6time,C6d,'k');
%plot(C10time,fits_C6-fits_C6(1),'b','linewidth',2);
hold on

    plot(time,cumpC6_5km,'linewidth',2);

    plot(time,cumpC6_10km,'linewidth',2);

    plot(time,cumpC6_20km,'linewidth',2);

    plot(time,cumpC6_40km,'linewidth',2);

ytickformat('%.2f');
datetick('x','keeplimits');
xlabel('Time');
ylabel('Pressure (KPa)');
title('C6');
xlim([time(1) time(end)])
ylim([-0.1 0.25])
%legend('data','5km','10km','20km','40km')

cx(4)=subplot(5,1,4);
pointsize=60;

    scatter(time_delP, best_fit5km/1000, pointsize, 's','filled');
    hold on
    scatter(time_delP, best_fit10km/1000, pointsize, 's','filled');
    scatter(time_delP, best_fit20km/1000, pointsize, 's','filled');
    scatter(time_delP, best_fit40km/1000, pointsize, 's','filled');

%      hline(-34000, 'r-', 'C2');
%      hline(-23000,'r-', 'C10');
%      hline(-2000, 'r-', 'C6');
     
datetick('x','keeplimits');
ax = gca;
ax.YRuler.Exponent = 0;
set(gca,'YDir','reverse')
xlabel('Time');
ylabel('KM');
title('Center of Slip Patch (km from trench)');
xlim([time(1) time(end)])
ylim([-45 5])
%legend('5km','10km','20km','40km')


cx(5)=subplot(5,1,5);
scatter(time_delP,slip_5km, pointsize, 's','filled')
hold on
scatter(time_delP,slip_10km, pointsize, 's','filled')
scatter(time_delP,slip_20km, pointsize, 's','filled')
scatter(time_delP,slip_40km, pointsize, 's','filled')
datetick('x','keeplimits');
ax = gca;
ax.YRuler.Exponent = 0;
xlabel('Time');
ylabel('cm');
title('slip');
xlim([time(1) time(end)])
ylim([-0.3 2.5])
linkaxes(cx,'x');
legend('5km','10km','20km','40km')

xlim([datenum('05-MAR-2020') datenum('15-APR-2020')])
end