description the koch snowflake is a fractal curve and one of the earliest fractals to have been described the koch snowflake can be built up iteratively in a sequence of stages the first stage is an equilateral triangle and each successive stage is formed by adding outward bends to each side of the previous stage making smaller equilateral triangles this can be achieved through the following steps for each line 1 divide the line segment into three segments of equal length 2 draw an equilateral triangle that has the middle segment from step 1 as its base and points outward 3 remove the line segment that is the base of the triangle from step 2 description adapted from https en wikipedia orgwikikochsnowflake for a more detailed explanation and an implementation in the processing language see https natureofcode combookchapter8fractals 84thekochcurveandthearraylisttechnique requirements pip matplotlib numpy initial triangle of koch snowflake uncomment for simple koch curve instead of koch snowflake initialvectors vector1 vector3 go through the number of iterations determined by the argument steps be careful with high values above 5 since the time to calculate increases exponentially iteratenumpy array0 0 numpy array1 0 1 array0 0 array0 33333333 0 array0 5 0 28867513 array0 66666667 0 array1 0 loops through each pair of adjacent vectors each line between two adjacent vectors is divided into 4 segments by adding 3 additional vectors inbetween the original two vectors the vector in the middle is constructed through a 60 degree rotation so it is bent outwards iterationstepnumpy array0 0 numpy array1 0 array0 0 array0 33333333 0 array0 5 0 28867513 array0 66666667 0 array1 0 standard rotation of a 2d vector with a rotation matrix see https en wikipedia orgwikirotationmatrix rotatenumpy array1 0 60 array0 5 0 8660254 rotatenumpy array1 0 90 array6 123234e17 1 000000e00 utility function to plot the vectors using matplotlib pyplot no doctest was implemented since this function does not have a return value avoid stretched display of graph matplotlib pyplot plot takes a list of all xcoordinates and a list of all ycoordinates as inputs which are constructed from the vectorlist using zip type ignore initial triangle of koch snowflake uncomment for simple koch curve instead of koch snowflake initial_vectors vector_1 vector_3 go through the number of iterations determined by the argument steps be careful with high values above 5 since the time to calculate increases exponentially iterate numpy array 0 0 numpy array 1 0 1 array 0 0 array 0 33333333 0 array 0 5 0 28867513 array 0 66666667 0 array 1 0 loops through each pair of adjacent vectors each line between two adjacent vectors is divided into 4 segments by adding 3 additional vectors in between the original two vectors the vector in the middle is constructed through a 60 degree rotation so it is bent outwards iteration_step numpy array 0 0 numpy array 1 0 array 0 0 array 0 33333333 0 array 0 5 0 28867513 array 0 66666667 0 array 1 0 standard rotation of a 2d vector with a rotation matrix see https en wikipedia org wiki rotation_matrix rotate numpy array 1 0 60 array 0 5 0 8660254 rotate numpy array 1 0 90 array 6 123234e 17 1 000000e 00 utility function to plot the vectors using matplotlib pyplot no doctest was implemented since this function does not have a return value avoid stretched display of graph matplotlib pyplot plot takes a list of all x coordinates and a list of all y coordinates as inputs which are constructed from the vector list using zip
from __future__ import annotations

import matplotlib.pyplot as plt  
import numpy


VECTOR_1 = numpy.array([0, 0])
VECTOR_2 = numpy.array([0.5, 0.8660254])
VECTOR_3 = numpy.array([1, 0])
INITIAL_VECTORS = [VECTOR_1, VECTOR_2, VECTOR_3, VECTOR_1]





def iterate(initial_vectors: list[numpy.ndarray], steps: int) -> list[numpy.ndarray]:
    
    vectors = initial_vectors
    for _ in range(steps):
        vectors = iteration_step(vectors)
    return vectors


def iteration_step(vectors: list[numpy.ndarray]) -> list[numpy.ndarray]:
    
    new_vectors = []
    for i, start_vector in enumerate(vectors[:-1]):
        end_vector = vectors[i + 1]
        new_vectors.append(start_vector)
        difference_vector = end_vector - start_vector
        new_vectors.append(start_vector + difference_vector / 3)
        new_vectors.append(
            start_vector + difference_vector / 3 + rotate(difference_vector / 3, 60)
        )
        new_vectors.append(start_vector + difference_vector * 2 / 3)
    new_vectors.append(vectors[-1])
    return new_vectors


def rotate(vector: numpy.ndarray, angle_in_degrees: float) -> numpy.ndarray:
    
    theta = numpy.radians(angle_in_degrees)
    c, s = numpy.cos(theta), numpy.sin(theta)
    rotation_matrix = numpy.array(((c, -s), (s, c)))
    return numpy.dot(rotation_matrix, vector)


def plot(vectors: list[numpy.ndarray]) -> None:
    
    
    axes = plt.gca()
    axes.set_aspect("equal")

    
    
    
    x_coordinates, y_coordinates = zip(*vectors)
    plt.plot(x_coordinates, y_coordinates)
    plt.show()


if __name__ == "__main__":
    import doctest

    doctest.testmod()

    processed_vectors = iterate(INITIAL_VECTORS, 5)
    plot(processed_vectors)