Plotting the Koch Snowflake

The Koch snowflake is a well known fractal. Beginning with an equilateral triangle, a smaller equilateral triangle is placed halfway along each edge of the shape. This process then repeats on each edge of the new shape. Here is my attempt to plot this in MATLAB.

```
clc, clearvars, close all;

length = 1; % Original side length of triangle
theta = 60*pi/180; % Angle of the equilateral triangle in radians

% Points of the original equilateral triangle
p1 = [0,0];
p2 = [1/2length,lengthsin(theta)];
p3 = [length,0];

flake_points = [p1; p2; p3]; % Vertices of the snowflake

for iteration = 1:7
n_points = size(flake_points,1); % No. points currently in data
new_triangle_points = nan(n_points*3,2); % New points to be added
for i = 1:n_points
p = flake_points(i,:);
if i == n_points
vector = flake_points(1,:)-flake_points(i,:);
else
% Vector between next point and current point
vector = flake_points(i+1,:)-flake_points(i,:);
end
b1 = p + 1/3*vector; % Base point 1 of the new triangle
b2 = p + 2/3*vector; % Base point 2 of the new triangle
top = p + 1/2vector + perp(vector)length/3*sin(theta);
new_triangle_points(3i-2:3i,:) = [b1; top; b2];
end
flake_points(end+1:end+9,:) = nan(9,2); % Add rows to data
new_flake_points = nan(n_points*3+n_points,2);
n = 1;
for i = 1:n_points
new_flake_points(n,:) = flake_points(i,:);
index = 3*i-2;
new_flake_points(n+1:n+3,:) = new_triangle_points(index:index+2,:);
n = n + 4; % Three new points have been inserted,
% so we now skip to the next empty spaces
end
flake_points = new_flake_points;
length = 1/3*length; % Revise length
end

plot(flake_points(:,1),flake_points(:,2),'w')
hold on; % Plot last data point with 1st (ie, close the loop)
plot([flake_points(1,1),flake_points(end,1)],...
[flake_points(1,2),flake_points(end,2)],'w')
set(gc

Very nice script!

I suggest you put this into a function, that takes n as an input argument instead of hardcoding it. It requires very little alterations, and will make the function a lot easier to use. You can skip all the clearing in the beginning and substitute it with:

function [] = koch_snowflake(iterations)  % or

and add an end at the bottom of your script.

length = 1; % Original side length of triangle

This is a no no. I strongly recommend you not to use length as a variable name as it's a very commonly used builtin function. This can cause a lot of errors that are hard to find. Note that you should always use the function numel instead of length, but it's still not a good idea to overload it.

theta = 60*pi/180; % Angle of the equilateral triangle in radians

I guess you think 60*pi/180 is easier to read than pi/3, but I disagree. If you want to use the number of degrees, and convert it to radians, you might want to use deg2rad.

I suggest substituting the iterator i to for instance ii if you ever work with complex numbers, to avoid potential errors and confusion. i is by default the imaginary unit in MATLAB (so is j, 1i, 1j). This is a matter of preference though.

Also, personal preference again, but I like it when comments are aligned (if possible):

length = 1;                   % Original side length of triangle
theta = 60*pi/180;            % Angle of the equilateral triangle in radians

% Points of the original equilateral triangle
p1 = [0,0];                   
p2 = [1/2*length,length*sin(theta)];
p3 = [length,0];              

flake_points = [p1; p2; p3];  % Vertices of the snowflake

for iteration = 1:7
    n_points = size(flake_points,1);          % No. points currently in data
    new_triangle_points = nan(n_points*3,2);  % New points to be added
...
...

StackExchange Code Review Q#144700, answer score: 3