% this script loads all data and then plots one or multiple earthquakes
% from experiment rough B. 

%% load and set parameters
load('roughB_variables.mat') % load variables previously defined in findearthquakes.m
load('roughness.mat') % load pattern seamounts for plotting
Earthquakes = importdata('Earthquake_table_roughB.txt'); % load earthquake data
firstframes = Earthquakes(:,1); % define first frames all eq's
lastframes = Earthquakes(:,3); % define last frames all eq's 
duration = Earthquakes(:,5); % define duration all eq's

  ROIx = 516.1; % ROI length in mm, taken from the PIV analysis
  ROIy = 307.0; % ROI width in mm, taken from the PIV analysis
  
% numbers of columns and rows (based on output PIV analysis)
  Nx=90;
  Ny=53;
  
  Dx = ROIx/Nx; % x distance one unit
  Dy = ROIy/Ny; % y distance one unit
  
  % distance backstop and top setup
  ROIback = 16.6; % distance ROI from backstop experiment in mm
  ROItop = 20.4; % distance ROI from top flank experiment in mm
  
  % Paramters seismogenic zone
  SZd = (305-ROIback); % horizontal distance downdip limit SZ
  SZu = (462.8-ROIback); % horizontal distance updip limit SZ
  Roughstart = round((295-ROIback)); % define starting point roughness pattern
  
  x = linspace(0,ROIx,Nx); 
  y = linspace(0,ROIy,Ny);
  [Xc,Yc] = meshgrid(x,y);
  
   slipthreshold = 0.05*0.13333; % eq threshold in cm
  
Dsum=zeros(Ny,Nx,length(firstframes));
Dsumall=zeros(Ny,Nx);
%% loop over eq's for plotting
%for i = 1:length(firstframes) %loop over all eq's
for i = 21 % or choose a single event from the earthquake table (example 21)
        clf
        disp('earthquake = ')
        disp(i)
     for j = firstframes(i):lastframes(i) %loop over eq frames
           disp(j)
           D=Vx(:,:,j)*13.333; %take displacement matrix in cm (*100*0.1333)
           Dsum(:,:,i)=Dsum(:,:,i)+D; % sum the displacement of all frames to Dsum
           Dsum(isnan(Dsum))=0; % convert NaN values to zero
     end
     % calculate source parameters
     Meanslip(i) = mean(mean(Dsum(:,:,i)))*10; % calculate max slip in mm
     Maxslip(i) = max(max(Dsum(:,:,i)))*10; % calculate max slip in mm
     Dsumfind = find (Dsum(:,:,i) > slipthreshold); % select nr of cells that exceed slip threshold
     Area(i)=length(Dsumfind)*(Dx/10)*(Dy/10); % multiply nr of cells with Dx and Dy to calculate area
     DsumfindSA = find (Dsum(:,:,i) > 0.05*Maxslip(i)); % find all cells that meet the seismic asperity criteria
     SA(i)=length(DsumfindSA)*(Dx/10)*(Dy/10); % multiply nr of cells with Dx and Dy to get seismic asperity area
         
     % plot cumulative slip distribution earthquake i
     contours=(slipthreshold:0.01:0.3); % start from 0.05 threshold
     contourf(Xc,Yc,Dsum(:,:,i),contours); % plot slip distribution
     hold on
     c = colorbar;
     caxis([0 0.3]); % set limits colorbar
     set(c,'fontsize',10); % set fontsize colorbar
     c.Label.String = 'Cumulated displacement (in cm)'; % set axis label colorbar
     seismicasperity=[0.05*Maxslip(i) 0.05*Maxslip(i)]; % determine seismic asperity contour
     contour(Xc,Yc,Dsum(:,:,i),seismicasperity,'w','LineWidth',1.5); % plot seismic asperity contour

     % plot Seismogenic zone
     line([SZd SZd],[0 ROIy], 'Color','k', 'LineWidth',1,'LineStyle','--')
     line([SZu SZu],[0 ROIy], 'Color','k', 'LineWidth',1,'LineStyle','--')
    
     axis equal
     % reverse y axis, since pivlab has the origin at the top left corner
     set(gca,'Ydir','reverse')
     xlim([0 ROIx]);
     ylim([0 ROIy]);
         
     % plot projection roughness pattern
     rate = 35.55/2700; % 35.55 cm/h = 9.875E-3 cm/s = 1.31666E-3 cm/f
     RoughExtra = zeros(Roughstart-(round(firstframes(i)*rate)),340);
     Zroughnew = [RoughExtra;Zrough];
     xrough = linspace(-ROItop,339-ROItop,340);
     yrough = linspace(0,(length(Zroughnew)-1),length(Zroughnew));
     [Xrough,Yrough] = meshgrid(xrough,yrough);
     contour(Yrough,Xrough,Zroughnew,[0.05 0.05],'LineColor','k');
       
               
     % save image for each rupture
     %    save_path = 'yourpath';
     %    filename = [save_path,num2str(i)];
     %    res = '-r300';
     %    format = '-djpeg';
     %    print(format,res,filename)
    
end
%%