% ### EXkoch1.m ###     10.14.14
% Plot Koch's 'Snowflake' (see http://mathworld.wolfram.com/KochSnowflake.html)
% [source - http://www.math.columbia.edu/~fusco/Hmwk01_APAM4300.pdf]

clear; clf;
% --------------------------------
% User Inputs
n = 4;      % 'order' of the snowflake
% --------------------------------

% ---
x = [-1/2 0 1/2 -1/2];  % starting matricies (for n=0)
y = [0 1 0 0];

% ---
% loop to create 'snowflake' - each time this outer loop executes, the
% inner loop updates to the next iteration
for zz = 1:n
    k = length(x);
    v = zeros(4*k-3,1);   % dummies (v-->x, w-->y) to fill up in inner loop;
    w = zeros(4*k-3,1);   % size stems from increase in resolution
    % this loop determines the 'iterated' version, each time updating based
    % upon the required rotations/reductions in step-size
    for j = 1:k-1,
        v(4*j-3) = x(j);
        w(4*j-3) = y(j);
        dirx = x(j+1) - x(j);
        diry = y(j+1) - y(j);
        v(4*j-2) = x(j) + 1/3*dirx;
        w(4*j-2) = y(j) + 1/3*diry;
        orthox = -diry;
        orthoy = dirx;
        v(4*j-1) = x(j) + 1/2*dirx + 1/3*1/2*sqrt(3)*orthox;
        w(4*j-1) = y(j) + 1/2*diry + 1/3*1/2*sqrt(3)*orthoy;
        v(4*j) = x(j) + 2/3*dirx;   % change fraction here to manipulate snowflake shape
        w(4*j) = y(j) + 2/3*diry;
    end
    v(4*k-3) = x(k);    % bring back to starting point
    w(4*k-3) = y(k);
    x = v; y = w;       % update x and v (re inner loop)
    
end
% ---
% visualize
plot(x,y,'LineWidth',2)
axis([-0.75 0.75 -sqrt(3)/6 1])