function DiffusionSimulation(data)
%% Show window with diffusion sumulation the real data
%% repsesnted as a set on polygons and calculated cylinders.
% After simullation graphs show up.
global TimePeriod GTimeStep np charge interCoef diffCoef simComplexity plotTimeStep
global figures stopSimulation
stopSimulation = false;
% Find geometry params
[cylinderRadii, convexHulls, realLevelHeight, cylinderHeight] = CalculateCylinders(data);
hulls3d = Calculate3DLevelHulls(convexHulls);
hullParams = FlatHullsTo3DHullParams(convexHulls);
cylinderRadiiSquares = cylinderRadii.^2;
cylinderGeometry = 1;
realGeometry = 2;
epsilonDistance = 1e-6;
% Number of levels
levelsSize = zeros(2, 1);
levelsSize(cylinderGeometry, 1) = size(cylinderRadii, 1);
levelsSize(realGeometry, 1) = size(hullParams, 1);
% Heights of levels
levelHeight = zeros(2, 1);
levelHeight(cylinderGeometry, 1) = cylinderHeight;
levelHeight(realGeometry, 1) = realLevelHeight;
% Plot colors
plotColor = cell(2, 1);
plotColor{cylinderGeometry} = 'blue';
plotColor{realGeometry} = 'red';
particlesLeftPlot = cell(2, 1);
% Simulation window
figure(figures.simulationWindow);
clf;
% Calculate plot boundaries
xyAxis = max(cylinderRadii);
for hull = 1:size(convexHulls, 1)
maxHullXY = max(max(abs(convexHulls{hull}(:, 1:2))));
if maxHullXY > xyAxis
xyAxis = maxHullXY;
end
end
xyAxis = xyAxis * 1.1;
% Setup plots and plot geometry
for geometry = [cylinderGeometry, realGeometry]
subplot(3, 2, [0, 2] + geometry);
view([0, 0]);
daspect([1, 1, 1]);
axis([-xyAxis, xyAxis, ...
-xyAxis, xyAxis, ...
0, levelsSize(geometry, 1) * levelHeight(geometry, 1) * 1.1]);
grid on;
hold on;
if (geometry == realGeometry)
PlotHulls(convexHulls);
title('Polygonal model');
elseif (geometry == cylinderGeometry)
PlotCylinders(cylinderRadii, cylinderHeight);
title('Cylinder model');
end
end
% Add controls
interCoefEdit = uicontrol('Style', 'edit', ...
'Position', [20, 20, 70, 20], ...
'String', interCoef);
uicontrol('Style', 'pushbutton', ...
'Position', [95, 20, 150, 20], ...
'String', 'Set interCoef and restart', ...
'UserData', {data, interCoefEdit}, ...
'Callback', @button_RestartSimulation_Callback);
uicontrol('Style', 'pushbutton', ...
'Position', [255, 20, 150, 20], ...
'String', 'Stop simulation', ...
'Callback', @button_StopSimmulation_Callback);
% Simulation params
StepRad = sqrt(6*diffCoef*GTimeStep); % Space step for the given diffusion coefficient
Time = 0:GTimeStep:TimePeriod; % Vector Time
NIteretionTime = size(Time, 2) - 1; % Number of Impulses/change in track
% Initialization of the position of particles..i.e origin
DifCoefTotal = cell(2, 1);
DifCoefTotal{1, 1} = zeros(1);
DifCoefTotal{2, 1} = zeros(1);
ParticleatFinish = zeros(2, 1);
TotalNumberinteration = 1;
for NumberIteretion = 1:TotalNumberinteration
startX = zeros(2, np);
startY = zeros(2, np);
startZ = zeros(2, np);
endX = zeros(2, np);
endY = zeros(2, np);
endZ = zeros(2, np);
% Statistics
particleLeftStatus = zeros(2, np);
particlesLeftThisStep = zeros(2, 1);
particlesLeftInTime = zeros(2, size(Time, 2));
particlesLeftInTime(:, 2:end) = NaN;
Distance = zeros(2, NIteretionTime+1);
TempDistance = zeros(2, 1);
TimeofOutControl = zeros(2, 1);
TimeofOut = zeros(2, NIteretionTime+1);
DifCoef = zeros(2, NIteretionTime+1);
% Plot initial point positions in subplots
pointPlots = zeros(2, 1);
for geometry = [cylinderGeometry, realGeometry]
subplot(3, 2, [0, 2] + geometry);
pointPlots(geometry, 1) = scatter3(...
startX(geometry, 1, :), ...
startY(geometry, 1, :), ...
startZ(geometry, 1, :), 12, 'filled', plotColor{geometry});
end
% Add particles-left plot
subplot(3, 2, [5, 6]);
ylabel('Particles passed');
xlabel('Time (ms)');
hold on;
for geometry = [cylinderGeometry, realGeometry]
particlesLeftPlot{geometry} = plot(Time, particlesLeftInTime(geometry, :), plotColor{geometry});
end
grid on;
hold off;
% Iterate over time
for i = 1:NIteretionTime % Time repeated
% Iterate over geometries
for geometry = [cylinderGeometry realGeometry]
% Iterate over particles
for j = 1:np % particles repeated
if stopSimulation
disp('Simulation stopped.');
return;
end
% Check if particle already left the geometry
if particleLeftStatus(geometry, j)
% Do nothing if particle is already left
endX(geometry, j) = startX(geometry, j);
endY(geometry, j) = startY(geometry, j);
endZ(geometry, j) = startZ(geometry, j);
else
% Calculate new position for the particle
% Simulate random movement
endX(geometry, j) = startX(geometry, j) - StepRad *randn();
endY(geometry, j) = startY(geometry, j) - StepRad *randn();
endZ(geometry, j) = startZ(geometry, j) - StepRad *randn();
% Add charge
if (startX(geometry, j)^2+startY(geometry, j)^2+startZ(geometry, j)^2) > 0.001^2 && charge ~= 0
endZ(geometry, j) = endZ(geometry, j) + charge * 0.1 * StepRad * (StepRad/0.2121) *startZ(geometry, j) / sqrt(startX(geometry, j)^2+startY(geometry, j)^2+startZ(geometry, j)^2);
end
% Choose between simulation complexity
if simComplexity == 1
% Simple simulation
if endZ(geometry, j) < 0
% Dont move particle if goes under the floor
endX(geometry, j) = startX(geometry, j);
endY(geometry, j) = startY(geometry, j);
endZ(geometry, j) = startZ(geometry, j);
else
particleLevel = floor(endZ(geometry, j)/levelHeight(geometry, 1)) + 1;
% Check if particle left the structure
if (particleLevel <= levelsSize(geometry, 1))
% Check if particle is inside the structure
if (geometry == cylinderGeometry)
particleInside = ((endX(geometry, j))^2+(endY(geometry, j))^2) < cylinderRadiiSquares(particleLevel);
elseif (geometry == realGeometry)
particleInside = IsPointInsideHull(hullParams{particleLevel, 1}, hullParams{particleLevel, 2}, endX(geometry, j), endY(geometry, j), endZ(geometry, j));
end
if ~particleInside
% Don't move particle if it's new
% cooridnate out of the structure
endX(geometry, j) = startX(geometry, j);
endY(geometry, j) = startY(geometry, j);
endZ(geometry, j) = startZ(geometry, j);
end
elseif particleLeftStatus(geometry, j) == 0
% Update statistics when particle left
particleLeftStatus(geometry, j) = 1;
TimeofOutControl(geometry, 1) = TimeofOutControl(geometry, 1) + i * GTimeStep;
ParticleatFinish(geometry, 1) = ParticleatFinish(geometry, 1) + 1;
particlesLeftThisStep(geometry, 1) = particlesLeftThisStep(geometry, 1) + 1;
end
end
else
% Complex simulation
if startZ(geometry, j) == 0 && endZ(geometry, j) < 0
% Dont move particle if starts and goes under the floor
endX(geometry, j) = startX(geometry, j);
endY(geometry, j) = startY(geometry, j);
endZ(geometry, j) = startZ(geometry, j);
else
particleLeft = false;
% Levels of the old and new position
particleLevel = floor(endZ(geometry, j)/levelHeight(geometry, 1))+1;
previousLevel = floor(startZ(geometry, j)/levelHeight(geometry, 1))+1;
% Current and new postition of the particle
preX = startX(geometry, j);
preY = startY(geometry, j);
preZ = startZ(geometry, j);
postX = endX(geometry, j);
postY = endY(geometry, j);
postZ = endZ(geometry, j);
% Coord change
xDiff = postX - preX;
yDiff = postY - preY;
zDiff = postZ - preZ;
% Check if particle is inside the geometry
if previousLevel > levelsSize(geometry, 1)
particleLeft = true;
else
% The level on which the particle
% should stay
leaveParticleAtLevel = previousLevel;
% Movement up/down
movesUp = particleLevel >= previousLevel;
% Find the level on which the particle
% should stay
if particleLevel ~= previousLevel
if movesUp
for level = previousLevel:min(particleLevel - 1, levelsSize(geometry, 1))
% Find intersection coords
pathPercentageBeforeFloorIntersection = (level * levelHeight(geometry, 1) - preZ) / zDiff;
floorIntersectionX = preX + xDiff * pathPercentageBeforeFloorIntersection;
floorIntersectionY = preY + yDiff * pathPercentageBeforeFloorIntersection;
% Check if particle passed
% through the gap between
% levels
if geometry == cylinderGeometry
if level < levelsSize(geometry, 1)
smallerRadiusSquare = min(cylinderRadiiSquares(level), cylinderRadiiSquares(level + 1));
else
smallerRadiusSquare = cylinderRadiiSquares(level);
end
intersectedSmallerRadius = (floorIntersectionX^2 + floorIntersectionY^2) <= smallerRadiusSquare;
else
inside = inpolygon(floorIntersectionX, floorIntersectionY, convexHulls{level + 1}(:, 1), convexHulls{level + 1}(:, 2));
intersectedSmallerRadius = inside > 0;
end
if ~intersectedSmallerRadius
leaveParticleAtLevel = level;
break;
else
leaveParticleAtLevel = level + 1;
end
end
else
for level = previousLevel:-1:max(particleLevel + 1, 1)
% Find intersection coords
pathPercentageBeforeFloorIntersection = ((level - 1) * levelHeight(geometry, 1) - preZ) / zDiff;
floorIntersectionX = preX + xDiff * pathPercentageBeforeFloorIntersection;
floorIntersectionY = preY + yDiff * pathPercentageBeforeFloorIntersection;
% Check if particle passed
% through the gap between
% levels
if geometry == cylinderGeometry
if level > 1
smallerRadiusSquare = min(cylinderRadiiSquares(level), cylinderRadiiSquares(level - 1));
else
smallerRadiusSquare = -1;
end
intersectedSmallerRadius = (floorIntersectionX^2 + floorIntersectionY^2) <= smallerRadiusSquare;
else
if level > 1
inside = inpolygon(floorIntersectionX, floorIntersectionY, convexHulls{level}(:, 1), convexHulls{level}(:, 2));
intersectedSmallerRadius = inside > 0;
else
intersectedSmallerRadius = false;
end
end
if ~intersectedSmallerRadius
leaveParticleAtLevel = level;
break;
else
leaveParticleAtLevel = level - 1;
end
end
end
end
% Found leaveParticleAtLevel
% Check if particle sould stay inside
% the structure
if leaveParticleAtLevel > levelsSize(geometry, 1)
particleLeft = true;
else
% Find the position of the particle
% inside its level
levelCeiling = leaveParticleAtLevel * levelHeight(geometry, 1);
levelFloor = (leaveParticleAtLevel - 1) * levelHeight(geometry, 1);
if movesUp
intersectedPlane = postZ >= levelCeiling;
intersectionHeight = levelCeiling;
else
intersectedPlane = postZ < levelFloor;
intersectionHeight = levelFloor;
end
particleHitWall = ~intersectedPlane;
% Check if particle left the floor or
% the ceiling.
if intersectedPlane
% Find intersection coords
pathPercentageBeforeFloorIntersection = (intersectionHeight - preZ) / zDiff;
floorIntersectionX = preX + xDiff * pathPercentageBeforeFloorIntersection;
floorIntersectionY = preY + yDiff * pathPercentageBeforeFloorIntersection;
% Find if particle intersected
% its level floor/ceiling area
if geometry == cylinderGeometry
intersectedRadius = (floorIntersectionX^2 + floorIntersectionY^2) <= cylinderRadiiSquares(leaveParticleAtLevel);
else
if movesUp
inside = inpolygon(floorIntersectionX, floorIntersectionY, convexHulls{leaveParticleAtLevel+1}(:, 1), convexHulls{leaveParticleAtLevel+1}(:, 2));
else
inside = inpolygon(floorIntersectionX, floorIntersectionY, convexHulls{leaveParticleAtLevel}(:, 1), convexHulls{leaveParticleAtLevel}(:, 2));
end
intersectedRadius = inside > 0;
end
% Check if particle hit the
% floor/ceiling
if intersectedRadius
% Left particle where it
% intersected the
% floor/ceiling
endX(geometry, j) = floorIntersectionX;
endY(geometry, j) = floorIntersectionY;
if movesUp
endZ(geometry, j) = intersectionHeight - epsilonDistance;
else
endZ(geometry, j) = intersectionHeight;
end
else
particleHitWall = true;
end
end
% Check if particle hit the wall.
% Otherwise its new coordinate is fine and
% shouldn't be changed.
if particleHitWall
if geometry == cylinderGeometry
% Cylinder geometry
if (postX^2 + postY^2) > cylinderRadiiSquares(leaveParticleAtLevel)
intersections = intersectLineCircle([preX preY 1 yDiff/xDiff], [0, 0, cylinderRadii(leaveParticleAtLevel)]);
interX = intersections(:, 1)';
interY = intersections(:, 2)';
if ((postX - interX(1))^2 + (postY - interY(1))^2) < ((postX - interX(2))^2 + (postY - interY(2))^2)
circleIntersectionX = interX(1);
circleIntersectionY = interY(1);
else
circleIntersectionX = interX(2);
circleIntersectionY = interY(2);
end
pathPercentageBeforeWallIntersection = sqrt(((circleIntersectionX - preX)^2 + (circleIntersectionY - preY)^2)/(xDiff^2 + yDiff^2));
circleIntersectionZ = preZ + zDiff * pathPercentageBeforeWallIntersection;
if circleIntersectionX < 0
circleIntersectionX = circleIntersectionX + epsilonDistance;
else
circleIntersectionX = circleIntersectionX - epsilonDistance;
end
if circleIntersectionY < 0
circleIntersectionY = circleIntersectionY + epsilonDistance;
else
circleIntersectionY = circleIntersectionY - epsilonDistance;
end
endX(geometry, j) = circleIntersectionX;
endY(geometry, j) = circleIntersectionY;
endZ(geometry, j) = circleIntersectionZ;
end
else
% Polygonal geometry
interPoints = zeros(0, 3);
polyCoords = hulls3d{leaveParticleAtLevel, 1};
polyIndices = hulls3d{leaveParticleAtLevel, 2};
% Find intersected polygons
for k = 1:size(polyIndices, 1)
vertex1 = polyCoords(polyIndices(k, 1), :);
vertex2 = polyCoords(polyIndices(k, 2), :);
vertex3 = polyCoords(polyIndices(k, 3), :);
maxDiff = max([xDiff, yDiff, zDiff]);
interPoint = intersectLineTriangle3d([preX, preY, preZ, xDiff/maxDiff, yDiff/maxDiff, zDiff/maxDiff], ...
[vertex1, vertex2, vertex3]);
if ~isnan(interPoint)
dotProduct = dot([xDiff yDiff zDiff], ([preX preY preZ] - interPoint));
if abs(dotProduct) > epsilonDistance && dotProduct < 0
interPoints = cat(1, interPoints, interPoint);
break;
end
end
end
numOfInters = size(interPoints, 1);
if numOfInters == 1
postDist = sum([xDiff yDiff zDiff].^2);
interDist = sum((interPoints(1, :) - [preX preY preZ]).^2);
% Chose closest point
if interDist < postDist
endX(geometry, j) = interPoints(1, 1);
endY(geometry, j) = interPoints(1, 2);
endZ(geometry, j) = interPoints(1, 3);
end
elseif numOfInters == 0
% Leave particle where
% it was
endX(geometry, j) = preX;
endY(geometry, j) = preY;
endZ(geometry, j) = preZ;
end
end
end
end
end
% Update statistics
if particleLeft && ~particleLeftStatus(geometry, j)
particleLeftStatus(geometry, j) = 1;
TimeofOutControl(geometry, 1) = TimeofOutControl(geometry, 1) + i * GTimeStep;
ParticleatFinish(geometry, 1) = ParticleatFinish(geometry, 1) + 1;
particlesLeftThisStep(geometry, 1) = particlesLeftThisStep(geometry, 1) + 1;
end
end
end
end
% Update statistics
TempDistance(geometry, 1) = TempDistance(geometry, 1) + (endX(geometry, j)^2 + endY(geometry, j)^2 + endZ(geometry, j)^2) / (6 * i * GTimeStep);
end % Repeate over particles
end
if (strcmp(get(figures.simulationWindow, 'Visible'), 'off'))
break;
end
% Update plots
for geometry = [cylinderGeometry, realGeometry]
% Update particles
set(pointPlots(geometry, 1),'XData',endX(geometry, :));
set(pointPlots(geometry, 1),'YData',endY(geometry, :));
set(pointPlots(geometry, 1),'ZData',endZ(geometry, :));
% Update plot
particlesLeftInTime(geometry, i+1) = particlesLeftThisStep(geometry, 1);
set(particlesLeftPlot{geometry}, 'XData', Time);
set(particlesLeftPlot{geometry}, 'YData', particlesLeftInTime(geometry, :));
% Update statistics
if ParticleatFinish(geometry, 1) > 0
TimeofOut(geometry, i+1) = ParticleatFinish(geometry, 1); %TimeofOutControl /ParticleatFinish;
else
TimeofOut(geometry, i+1) = 0;
end
ParticleatFinish(geometry, 1) = 0;
Distance(geometry, i+1) = (TempDistance(geometry, 1)^2) / np;
DifCoef(geometry, i+1) = (TempDistance(geometry, 1) / np / 3) / 10^6;
TempDistance(geometry, 1) = 0;
DifCoefTotal{geometry, 1} = DifCoefTotal{geometry, 1} + DifCoef(geometry, :);
end
pause(0.0000005);
startX = endX;
startY = endY;
startZ = endZ;
end
end
if (strcmp(get(figures.simulationWindow, 'Visible'), 'off'))
return;
end
% Plot resulting graphs
figure(figures.graphDistance);
clf;
grid on;
figure(figures.graphDiffCoef);
clf;
grid on;
figure(figures.graphTimeOfOut);
clf;
grid on;
plotTime = 0:plotTimeStep:TimePeriod;
plotTimeSize = size(plotTime, 2);
for geometry = [cylinderGeometry realGeometry]
figure(figures.graphDistance);
hold on;
plot(plotTime, ReduceVectorToSize(Distance(geometry, :), plotTimeSize), plotColor{geometry});
xlabel('Time');
ylabel('Distance');
figure(figures.graphDiffCoef);
hold on;
plot(plotTime, ReduceVectorToSize(DifCoefTotal{geometry, :}, plotTimeSize)/TotalNumberinteration, plotColor{geometry});
xlabel('Time');
ylabel('Diff coef.');
figure(figures.graphTimeOfOut);
hold on;
% Create plot
reduced = IntegrateVectorToSize(particlesLeftInTime(geometry, :), plotTimeSize, GTimeStep, plotTimeStep);
for i = 2:size(reduced, 2)
plot(plotTime, reduced, plotColor{geometry});
end
xlabel('Time, ms');
ylabel('Particles per ms');
end
end