% 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