% Create a random polygon at the bottom half to act as a bay/water inlet.
% The land mass is the top half.
% The line separating these two polygons is the shore line.
% Generate some-what predictable moving particles (points) bouncing
% around the 'bay'. When the particles hit the shore line, high light
% an exaggerated (more than the difference between the left and right
% most particles) affected area. Notes: the left/right/bottom of bay
% is considered as open water. They are used just for the purpose of
% turning the particles back into the bay.

close all;
clear all;
clc;

%  Simulation time
startIteration = 0;
iterationIncrement = 0.1;
endIteration = 100;
f=0;

% Simulation boundaries
xBoundary = 100;
yBoundary = 100;

% Derive change/movement from boundary above
xVelocity = xBoundary / (endIteration - startIteration) / iterationIncrement * 3;  %  meters/second
yVelocity = yBoundary / (endIteration - startIteration) / iterationIncrement * 7;  %  meters/second

% Generate initial area
pSize = 5;         % particle size

nParticles = 3;     % number of particles for moving patch

xSpill = zeros(1, nParticles);
ySpill = zeros(1, nParticles);

A = zeros(nParticles, nParticles);

n = 0;
loopCount = 0;

while n < nParticles
    x = mod(round(rand * 100 / nParticles), nParticles) + 1;
    y = mod(round(rand * 100 / nParticles), nParticles) + 1;
    if A(x,y) == 0
       A(x,y) = 1;
       n = n+1;
    end
end

n = 0;
while n < nParticles
    for x=1:nParticles
        for y=1:nParticles
            if A(x,y) == 1
                n = n+1;
                xSpill(n) = x;
                ySpill(n) = y;
            end
        end
    end
end

% Generate water way/inlet/bay

nPoints = floor(xBoundary / 10) + 2;
xBay=zeros(1, nPoints);
yBay=zeros(1, nPoints);

%% Create a zigzag shoreline
lowest = round(yBoundary/3);
xBay(1) = 1;
yBay(1) = 1;
xBay(2) = 1;
yBay(2) = lowest;
n = 2;

while n < nPoints - 1
    % Some pseudo random curve/line
    yt = mod(round(yBoundary / (rand * yBoundary)), yBoundary) + 1;
    if yt >= yBay(n-1)
        n = n+1;
        % Spacing out the X´s evenly
        xBay(n) = xBay(n-1)+nPoints;
        if yt < lowest
            yBay(n) = lowest;
        else
            yBay(n) = yt;
        end
    end
end

% Making the top right boundary line up to right most of grid

xBay(nPoints-1) = xBoundary;

xBay(nPoints) = xBoundary;      % close the polygon
yBay(nPoints) = 1;

% Create shore line

xShore = xBay(2:end-2);
yShore = yBay(2:end-2);

% Create land mass polygon

for n = 1:nPoints-2;
    xLand(n) = xBay(n+1);
    yLand(n) = yBay(n+1);
end;

xLand(nPoints-1) = xBoundary;
yLand(nPoints-1) = yBoundary;
xLand(nPoints) = 1;
yLand(nPoints) = yBoundary;

% Create x/y matrices to store affected areas

xAffectedArea = [];
yAffectedArea = [];

% end of setup.

% Run the simulation:

nCount = 0;
nMaxCount = xBoundary * 10;

% Simulate up to the end of iteration or each maximum number of nMaxCount
while (startIteration < endIteration) && (nCount < nMaxCount)
    
    % tracking loop count
    nCount = nCount + 1;

    axis([-1,xBoundary+1,-1,yBoundary+1]);
    axis manual;
    axis equal;
    axis tight;
    hold on;

    % Plotting the land, the bay, spill, and affected area
    fill(xLand, yLand, 'g');

    fill(xBay, yBay, 'b');

    plot(xSpill, ySpill, 'ro', 'MarkerSize', pSize);

    plot(xAffectedArea, yAffectedArea, 'r', 'linewidth', 5);

    drawnow;
  
    % increase iteration counter
    startIteration = startIteration + iterationIncrement;

    % Moving along based on some random velocity
    % TODO - it should be some current movement algorithm to be implemented here.
    xSpill = xSpill + xVelocity * iterationIncrement;
    ySpill = ySpill + yVelocity * iterationIncrement;

    hittingBoundary = 0;

    [inBay onLand] = inpolygon(xSpill, ySpill, xBay, yBay);

    while (~inBay | onLand) & (hittingBoundary < nMaxCount)
        if hittingBoundary == 0
            % only if it hit the top of polygons
            if max(ySpill) >= lowest
                % Go to the left to create exaggerated effect
                x1 = min(xSpill);
                if (x1 < 1)
                    x1 = 1;
                end
                if (x1 > xBoundary)
                    x1 = xBoundary;
                end

                % Go to the right to create exaggerated effect
                x2 = max(xSpill);
                if (x2 < 1)
                    x2 = 1;
                end
                if (x2 > xBoundary)
                    x2 = xBoundary;
                end

                xx = [x1 x2];
                yy = [];

                % Get the x1 and x2's y locations
                yy = [yy interp1(xShore, yShore, x1)];
                yy = [yy interp1(xShore, yShore, x2)];

                plot(xx, yy, 'r', 'linewidth', 5);
                drawnow;

                % Adding to the overall affected area for drawing on next loop
                xAffectedArea = [xAffectedArea xx NaN];
                yAffectedArea = [yAffectedArea yy NaN];
            end
        end

        hittingBoundary = hittingBoundary + 1;

        % change direction
        y1 = interp1(xShore, yShore, min(xSpill));
        y2 = interp1(xShore, yShore, max(xSpill));
        minYSpill = min(ySpill);
        maxYSpill = max(ySpill);

        if (y1 == NaN) | (y2 == NaN) | (max(xSpill) >= xBoundary) | (min(xSpill) <= 1)
            xVelocity = -1 * xVelocity;
        end

        if (y1 == NaN) | (y2 == NaN) | (minYSpill >= y1) | (minYSpill >= y2) | (maxYSpill >= y1) | (maxYSpill >= y2) | (minYSpill <= 1)
            yVelocity = -1 * yVelocity;
        end

        xSpill = xSpill + (xVelocity * iterationIncrement * hittingBoundary);
        ySpill = ySpill + (yVelocity * iterationIncrement * hittingBoundary);
        [inBay onLand] = inpolygon(xSpill, ySpill, xBay, yBay);
    end

    if hittingBoundary >= nMaxCount
        break;
    end
    %hold off;
    f=f+1;
    saveas(gcf,['oil_spill_simulation1', sprintf('%03d',f),'.png'])
  
end

