% Earthquake simulation

% 1. create a grid 100x100m long, sampling 1 m
time       = 80; % seconds
SouthNorth = 100;
EastWest   = 100;
dx         = 1;
canvas = zeros(SouthNorth,EastWest);
time_sampling = 1;
measurement_time = 0:1/time_sampling:time-1/time_sampling;

% 2. create the positions for the optical fiber (L shaped)
fiber_position_SouthNorth = [21,80];
fiber_position_WestEast   = [21,80];
fiber_plot = canvas;
fiber_plot(fiber_position_SouthNorth(1):fiber_position_SouthNorth(2),fiber_position_WestEast(1)) = 1;
fiber_plot(fiber_position_SouthNorth(2),fiber_position_WestEast(1):fiber_position_WestEast(2))  = 1;

% figure, imagesc(fiber_plot);
nVect = [canvas(fiber_position_SouthNorth(1):fiber_position_SouthNorth(2),fiber_position_WestEast(1))' canvas(fiber_position_SouthNorth(2),fiber_position_WestEast(1):fiber_position_WestEast(2))];
nVect = repmat(nVect,[length(measurement_time),1]);

wavenumber     = 1; %rad/m
phase_velocity = 6;
omeg = wavenumber*phase_velocity;

angle_deg=85;

main_dir    = [cosd(angle_deg), sind(angle_deg)];
rotMatrix   = [cosd(90) -sind(90); sind(90) cosd(90)];
ortogVector = main_dir*rotMatrix;

idx = 1;
h = figure(3);
set(h, 'WindowStyle', 'Docked');

max_pert_length = ceil(sqrt(size(canvas,1)^2+size(canvas,2)^2));
ps = 0:max_pert_length;
for t = measurement_time

    perturbation = sin(omeg*t - wavenumber*dx*ps).*exp(-((ps-phase_velocity*t).^2)/32);

    pert_points = 1:0.5:(length(perturbation)-1);

    for idd = 1:(length(pert_points)-1)

        main_pos = floor(pert_points(idd)*main_dir);
        idd2 = -100;

        while idd2 < 100
            p = main_pos+floor(idd2*ortogVector)+1;
            if all(p    > 0)          && ...
               all(p(2) < EastWest)   && ...
               all(p(1) < SouthNorth)
                canvas(p(1),p(2)) = perturbation(floor(pert_points(idd)));
            end
            idd2 = idd2 + 1;
        end
    end
    subplot(1,2,1), imagesc(canvas+fiber_plot), drawnow;
    idx = idx+1;
    nVect(idx,:) = [canvas(fiber_position_SouthNorth(1):fiber_position_SouthNorth(2),fiber_position_WestEast(1))' canvas(fiber_position_SouthNorth(2),fiber_position_WestEast(1):fiber_position_WestEast(2))];
    subplot(1,2,2), plot(nVect(idx,:)),   axis([0 size(nVect,2) -1 1]),    drawnow;
end
%
% y = mx+b
% function [fiber_canvas] = generateFiberPositionMatrix(canvas, nodes)
%
% y_vec = 1:size(canvas,1)
% x_vec = 1:size(canvas,2)
%
% % 2. create the positions for the optical fiber (L shaped)
% fiber_position_SouthNorth = [21,80];
% fiber_position_WestEast   = [21,80];
% fiber_plot = canvas;
% fiber_plot(fiber_position_SouthNorth(1):fiber_position_SouthNorth(2),fiber_position_WestEast(1)) = 1;
% fiber_plot(fiber_position_SouthNorth(2),fiber_position_WestEast(1):fiber_position_WestEast(2))  = 1;
%
% % figure, imagesc(fiber_plot);
% nVect = [canvas(fiber_position_SouthNorth(1):fiber_position_SouthNorth(2),fiber_position_WestEast(1))' canvas(fiber_position_SouthNorth(2),fiber_position_WestEast(1):fiber_position_WestEast(2))];
% nVect = repmat(nVect,[length(measurement_time),1]);
%
% prev_node = nodes(1,:)
%     for node_idx = 2:size(nodes,1)
%         node = nodes(node_idx,:)
%         m = (node(2)-prev_node(2))/(node(1)-prev_node(1));
%         nVect = [nVect canvas(
%     end
% end