% This script loads all experimental data and calculates a cumulative
% displacement curve of one point centered above the seismogenic zone.
%
% see Figure 2 in van Rijsingen et al. (2019). Rough subducting seafloor
% reduces interseismic coupling and mega-earthquake occurrence: insights
% from analogue models. Geophysical Research letters.

%% load Smooth A
load('smoothA_variables.mat') % load variables in m/s
% take value at centered point, convert to displacement in mm
cumdis = cumsum(Vx(25,64,:)*0.13333*-1000,3);

% loop over all frames to create array
for i = 1:12853
cumdis_smoothA(i)=cumdis(:,:,i);
end

clearvars -except cumdis_smoothA 

%% load Smooth B
load('smoothB_variables.mat') % load variables in m/s
% take value at centered point, convert to displacement in mm
cumdis = cumsum(Vx(29,68,:)*0.13333*-1000,3);

% loop over all frames to create array
for i = 1:10011
cumdis_smoothB(i)=cumdis(:,:,i);
end

clearvars -except cumdis_smoothA cumdis_smoothB

%% load smooth C
load('smoothC_variables.mat') % load variables in m/s
% take value at centered point, convert to displacement in mm
cumdis = cumsum(Vx(26,63,:)*0.13333*-1000,3);

% loop over all frames to create array
for i = 1:8490
cumdis_smoothC(i)=cumdis(:,:,i);
end

clearvars -except cumdis_smoothA cumdis_smoothB cumdis_smoothC

%% load smooth D
load('smoothD_variables.mat') % load variables in m/s
% take value at centered point, convert to displacement in mm
cumdis = cumsum(Vx(24,58,:)*0.13333*-1000,3); % take values at centered point

% loop over all frames to create array
for i = 1:9873 
cumdis_smoothD(i)=cumdis(:,:,i); 
end

% clear all variables except the cumdis arrays
clearvars -except cumdis_smoothA cumdis_smoothB cumdis_smoothC cumdis_smoothD


%% load Rough A
load('roughA_variables.mat') % load variables in m/s
% take value at centered point, convert to displacement in mm
cumdis = cumsum(Vx(21,52,:)*0.13333*-1000,3);

% loop over all frames to create array
for i = 1:10108
cumdis_roughA(i)=cumdis(:,:,i);
end

% clear all variables except the cumdis arrays
clearvars -except cumdis_smoothA cumdis_smoothB cumdis_smoothC cumdis_smoothD cumdis_roughA

%% load Rough B
load('roughB_variables.mat') % load variables in m/s
% take value at centered point, convert to displacement in mm
cumdis = cumsum(Vx(27,64,:)*0.13333*-1000,3);

% loop over all frames to create array
for i = 1:10244
cumdis_roughB(i)=cumdis(:,:,i);
end

clearvars -except cumdis_smoothA cumdis_smoothB cumdis_smoothC cumdis_smoothD cumdis_roughA cumdis_roughB

%% load Rough C
load('roughC_variables.mat') % load variables in m/s
% take value at centered point, convert to displacement in mm
cumdis = cumsum(Vx(25,59,:)*0.13333*-1000,3);

% loop over all frames to create array
for i = 1:9413
cumdis_roughC(i)=cumdis(:,:,i);
end

clearvars -except cumdis_smoothA cumdis_smoothB cumdis_smoothC cumdis_smoothD cumdis_roughA cumdis_roughB cumdis_roughC

%% load Rough D
load('roughD_variables.mat')% load variables in m/s
% take value at centered point, convert to displacement in mm
cumdis = cumsum(Vx(24,59,:)*0.13333*-1000,3); % select point in the middle of the seismogenic zone

% loop over all frames to create array
for i = 1:10040
cumdis_roughD(i)=cumdis(:,:,i);
end

clearvars -except cumdis_smoothA cumdis_smoothB cumdis_smoothC cumdis_smoothD cumdis_roughA cumdis_roughB cumdis_roughC cumdis_roughD

%% Plot data
minutes = [0:1:7]; % create time array for movement downgoing plate
dplate = 5.835*minutes; % create velocity array for downgoing plate
time = linspace(0,8490/450,8490); % create time array for cumdis

plot(time,cumdis_smoothA(1:8490)) % plot smooth A for 8490 frames
hold on
plot(time,cumdis_smoothB(1:8490)) % plot smooth B for 8490 frames
plot(time,cumdis_smoothC(1:8490)) % plot smooth C for 8490 frames
plot(time,cumdis_smoothD(1:8490)) % plot smooth D for 8490 frames
plot(time,cumdis_roughA(1:8490)) % plot rough A for 8490 frames
plot(time,cumdis_roughB(1:8490)) % plot rough B for 8490 frames
plot(time,cumdis_roughC(1:8490)) % plot rough C for 8490 frames
plot(time,cumdis_roughD(1:8490)) % plot rough D for 8490 frames
plot(minutes,dplate) % plot movement downgoing plate
legend('smooth A','smooth B','smooth C','smooth D','rough A','rough B','rough C','rough D','downgoing plate')
xlabel('time (min)')
ylabel('cumulative displacement point middle SZ (mm)')
xlim([0 22]);
ylim([0 50]);

