%% Plot Generator for "Seismological Stress Drops for Confined Ruptures are Invariant To Normal Stress
% Will Steinhardt | UCSC | 2022
% For more information contact: wsteinh@gmail.com

close all;
clear all;

%% Load directtoy of event data
fileList = dir('AllEventData*.mat');
nsets = max(size(fileList));

% % % Plot Table of Experiment Info
% for kk = 1:nsets
%     load(fileList(kk).name);
%     tabl(kk).num = kk;
%     tabl(kk).desc = desc;
% end
% struct2table(tabl)

nsets = max(size(fileList));
lastFile = max(size(fileList));

powerfit3 = @(b,x) ( b(1).*x.^3 );
beta0 = [1e3];
nameStruct = [];


%% Clear figures and set positions
figure(1); clf; hold on; set(gcf,'Position',[200 200 1300 1050]); box on;
figure(3); clf; box on; set(gcf,'Position', [1505 325 560 420])
figure(4); clf; box on; set(gcf,'Position', [1950 575 540 410])
figure(5); clf; box on; set(gcf,'Position', [1505 830 560 420])

% Grouped lists of experiments under common experimental conditions, sorted in order of normal stress
noAsp = [2 4 3 5 6 7 32 33 34]; n_noAsp = max(size(noAsp));
cAspL = [10 11 12 13 14]; n_cAspL = max(size(cAspL));
thickSand = [40 41 42 43]; n_thickSand = max(size(thickSand));
noSand = [30 31 19 20 21]; n_noSand = max(size(noSand));
fric = [thickSand noSand]; n_fric = max(size(fric));
sample5x5 = [16 44 45 17]; n_5x5 = max(size(sample5x5));
sample3x3 = [47 46 23 15 24 18 26]; n_3x3 = max(size(sample3x3));
smallCirc = [27 28 29 8]; n_circ = max(size(smallCirc));
smallSamps = [sample5x5 sample3x3 smallCirc]; n_smallSamps = max(size(smallSamps));

exptNumMat = [noAsp cAspL fric smallSamps];
vertLines = 0.5 + cumsum([n_noAsp n_cAspL n_thickSand n_noSand n_5x5 n_3x3 n_circ]);

showEvents = 1:max(size(exptNumMat));
nsets = max(showEvents);

binEdges = linspace(0,3e4,25);
meanSDmat = nan(length(showEvents),2);

expFitMat = [];
avgAreaMat = [];
maxAreaMat = [];

for kk = showEvents

    %% Load data file, set plotting info, check if there are any found events for the datafile
    jj = exptNumMat(kk);
    load(fileList(jj).name);

    tag = fileList(jj).name(25:end-4);
    nameStruct(kk).name = desc;

    mks = 10;
    filledYN = 'y';

    if isempty(allEvents)
        continue;
    end

    %% Initailize variables for cutoff and checking if there any any points near the edge of the imaging window
    slipAreas = vertcat(allEvents.SlipArea2);
    slipLengths = sqrt(slipAreas./pi);
    moments = vertcat(allEvents.Moment);

    %% Remove events below defined cutoff radius
    cutoff = 1e-3;
    ind_low = find(slipLengths < cutoff);
    allEvents(ind_low) = [];

    %% Check if any of the perimeter points in the event reach the edge of the window, and if so, remove them from the catalog
    nevents = size(allEvents,2);
    datDims = allEvents(1).DataDimensions;
    datDimsX = datDims(1); datDimsY = datDims(2);
    edge_ind = [];

    for nn = 1:nevents
        shp = allEvents(nn).Perim2;
        verts = shp.Vertices;

        maxX = max(verts(:,1));
        maxY = max(verts(:,2));

        xedgepts = find(verts(:,1) >= datDimsX-1);
        yedgepts = find(verts(:,2) >= datDimsY-1);

        uniqueEdgePts = unique([xedgepts; yedgepts]);

        if isempty(maxX) || isempty(maxY)
            continue;
        end

        nEdgePts = size(uniqueEdgePts,1);
        pctEdgePts = nEdgePts / size(verts,1);

        if nEdgePts > 0
            edge_ind = [edge_ind; nn];
        end

    end

    allEvents(edge_ind) = [];

    if isempty(allEvents)
        continue;
    end

    %% define variables for calculation and plotting
    durations = vertcat(allEvents.Duration);
    slipAreas = vertcat(allEvents.SlipArea2);
    slipLengths = sqrt(slipAreas./pi);
    moments = vertcat(allEvents.Moment);
    durations = vertcat(allEvents.Duration);

    %% Plot Moment as a function of rupture length
    figure(5);
    plot(slipLengths,moments,symb,'MarkerSize',mks-5,'MarkerFaceColor',colN,'MarkerEdgeColor',colN); hold on;
    xline(cutoff,'--k');
    ylim([0.6e-5 4e-1]); xlim([(cutoff-cutoff/4), 2.1e-2]);
    set(gca,'XScale','log','YScale','log')

    %% Calculate stress drop for each event
    stressDrops = 7/16 * moments ./ (slipLengths.^3);
    meanSD = mean(stressDrops);
    stderrSD = std(stressDrops) ./ sqrt(length(stressDrops));
    meanSDmat(kk,:) = [kk meanSD];

    %% Plot average stress drop
    figure(1); subplot(5,1,[3 4])
    errorbar(kk, meanSD, stderrSD, symb, 'MarkerSize',8,'Color',colN,'LineWidth',2,'MarkerFaceColor',colN); hold on;

    %% Generate violin plot of stress drops
    [counts binEdges] = histcounts(stressDrops,50,'Normalization','pdf');
    binCenters = (binEdges(1:end-1)+binEdges(2:end))./2;

    figure(1); subplot(5,1,[1 2]); hold on;
    countsW=counts.*5700;
    xpts = [kk-countsW'; flipud(kk+countsW')];
    ypts = [binCenters'; flipud(binCenters')];
    p=fill(xpts,ypts,'r'); hold on;
    p.FaceColor = colN; p.FaceAlpha = 0.75; p.EdgeAlpha = 0;

    meanSD = mean(stressDrops);
    stderrSD = std(stressDrops) ./ sqrt(length(stressDrops));
    meanSDmat(kk,:) = [kk meanSD];

    plot(kk,meanSD,'o','MarkerSize',7,'MarkerEdgeColor',[0.8 0.8 0.8],'MarkerFaceColor',[0 0 0],'LineWidth',2)

    %% Plot the number of events for each experiment
    figure(4); hold on;
    if strcmp(filledYN,'y')
        plot(kk, size(allEvents,2), symb, 'MarkerSize',8,'Color',colN,'LineWidth',2,'MarkerFaceColor',colN);
    else
        plot(kk, size(allEvents,2), symb, 'MarkerSize',8,'Color',colN,'LineWidth',2);
    end

    %% Make b-value plots for each dataset
    eqMags = (log10(moments)-9.1)./1.5;
    minMag = min(eqMags);
    nbins = 15;
    bins = linspace(minMag,-7.5,nbins);

    if jj == 8
        bins = linspace(minMag,max(eqMags),nbins);
    end

    [counts binEdges] = histcounts(eqMags, bins);
    counts = cumsum(counts,'reverse');
    minBinSize = 1;

    if jj == 40
        minBinSize = 1;
    end

    minbin = find(counts < minBinSize,1);

    if ~isempty(minbin)
        minEdge = binEdges(minbin+1);
        eqMagsThresh = eqMags(eqMags < minEdge);

        bins = linspace(minMag,minEdge,nbins);
        [counts binEdges] = histcounts(eqMagsThresh, bins);

        counts = cumsum(counts,'reverse');

        minbin = find(counts < minBinSize,1);
        minEdge = binEdges(minbin+1);
        eqMagsThresh = eqMags(eqMags < minEdge);
    end

    binCenters = (binEdges(1:end-1)+binEdges(2:end))./2;
    yvals = log10(counts./size(eqMags,1));

    figure(3); hold on;
    plot(binCenters,yvals, symb, 'MarkerSize',5,'Color',colN,'LineWidth',2,'MarkerFaceColor',colN); hold on;

    %% Plot normal force and normal stress at each experiment
    normStress = str2num(desc(end-12:end-9)) * 10; % factor of 10 converts to kN / m^2

    figure(1); subplot(5,1,5)
    plot(kk, normStress,'.k','MarkerSize',20); hold on;

end

vertLines = vertLines(vertLines < nsets);

if isempty(vertLines)
    vertLines = [0];
end

for ii=1:max(size(vertLines))
    nameStruct(nsets+ii).name = '';
end

figure(1); subplot(5,1,[3 4]);
xlim([0 nsets+1]);
ylim([0.4e4 3e4])
ylabel('Average Stress Drop (Pa)'); xlabel('Experiment');
set(gca,'FontName','Arial','FontSize',12,'FontWeight','Bold','LineWidth',1.5)
xline(vertLines,'-k','LineWidth',2)
xticks([])
subplot(5,1,[1 2]); box on;
xlim([0 nsets+1]);
ylabel('Stress Drop (Pa)'); xlabel('Experiment #');
set(gca,'FontName','Arial','FontSize',12,'FontWeight','Bold','LineWidth',1.5)
xline(vertLines,'-k','LineWidth',2)
xticks([])
set(gca,'YScale','log')
subplot(5,1,5)
xticks([]);
xlim([0 nsets+1]);
ylabel('\sigma_n (kN/m^2)');
set(gca,'FontName','Arial','FontSize',12,'FontWeight','Bold','LineWidth',1.5)
xline(vertLines,'-k','LineWidth',2)

figure(3);
xlabel('Magnitude'); ylabel('log(N/N_{tot})');
set(gca,'FontName','Arial','FontSize',12,'FontWeight','Bold','LineWidth',1.5)
xvals = linspace(-11,-7,100);
plot(xvals,-xvals-11,'--k','LineWidth',5)
plot(xvals,-1.5*xvals-15.5,'--r','LineWidth',5)
plot(xvals,-2*xvals-20.5,'--b','LineWidth',5)
ylim([-3 0.5]); xlim([-9.5 -6.5])

figure(5)
xlabel('Rupture Radius (m)'); ylabel('Moment (N m)');
set(gca,'FontName','Arial','FontSize',12,'FontWeight','Bold','LineWidth',1.5)
set(gca,'XScale','log','YScale','log')

figure(4)
xlabel('Experiment #'); ylabel('n_{events}');
set(gca,'FontName','Arial','FontSize',12,'FontWeight','Bold','LineWidth',1.5)