Building fractals with the Julia programming language
Tags: julialang, fractal, plotting
Published: 2022-03-22 22:26:11

Hello everybody,

let me introduce myself first:
- Highschool student.
- Pretty much always had a passion to math.
- Started programming around a year ago.
- French (sorry if my English sounds sketchy)

As I always liked maths, I wondered how I could have fun while programming, something that I would find fun, and that I could share to others even if they don't know the (technical) details.
That leaves plenty of rooms for choice, but I end up doing lots of fractals, as the visual aspect makes it especially rewarding.

I do it in Julia, high-level, fast programming language. Be aware that this isn't a tutorial to Julia per se, a little bit on its libraries, it's just me showing off stuffs I did and how.

Packages used:
- Luxor.jl
- GLMakie.jl
- CliffordAlgebras.jl
- Lindenmayer.jl

Let's start really simple, and 2D, as a warm-up, let's just do a spiral. For that, we will plot every points on this spiral and join them. And let the user choose how many circle they want (it would be infinity to make a real fractal).

```julia
using GLMakie


function spiral_point(n)
	rot = (ℯ^(im*n)) # set up the rotation of the point
	point = rot/n
	(point.re, point.im)  # this point converge to (0, 0) as n grows
end

function spiral(n) # where n is the number of spirals
	V = Tuple{Float64, Float64}[]
	for i in 1:0.01:2π*n
		push!(V, spiral_point(i))
	end
	V
end

println("How many spirals do you want? (number in N)")
const num = parse(Int, readline()) 
spiral(num) |> lines |> display
```
![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/g50rvlxhaysx5jfoe4zx.png)

I personally don't like that it doesn't make enough detail changes as n grows, would be cooler if it was slower first, then faster, so that it shows more details. Let's do a gradient function for that.

```julia
using GLMakie
function gradient(n)
	f(x) =  √(x/n)
end

function spiral_point(n, f)
	rot = (ℯ^(im*n))
	point = rot/f(n)
	(point.re, point.im)
end 

function spiral(n, f)
	V = Tuple{Float64, Float64}[]
	for i in 1:0.01:2π*n
		push!(V, spiral_point(i, f))
	end
	V
end

println("How many spirals do you want? (number in N)")
const num = parse(Int, readline())
spiral(num, gradient(num)) |> lines |> display
```
![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/kb9xfcpr0h5b8qjwutic.png)

Looks much better. Useful technique for customization.

Recursion is a common technic to create fractals, let me show an example which looks really cool when n is small, and really familiar when n is big :)
```julia
using Luxor

const N = 9 # change that to change the deps
const SIZE = 2700 # change that to resize the image
const NAME = "MyLovelyFractal.png" # change that to change the name of your .png file you'll receive

function points(x, y, size, l=Tuple{Int, Int}[])
	0 == size && return x, y, [l; (x, y)]
	x, y, l = points(x, y, size-1, l)
	y -= 2^(size-1)
	x -= 2^(size-1) - 1
	x, y, l = points(x, y, size-1, l)
	x += 1
	y += 2^(size-1)
	points(x, y, size-1, l)
end
p(s) = points(0, 0, s)[end]

function main()
	v = p(NUM)
	Drawing(SIZE, SIZE, NAME)
	origin(0, SIZE)
	setcolor("blue")
	for i in 1:length(v)-1
		line(Point(v[i])*5, Point(v[i+1])*5, :stroke)
	end
	finish()
	preview()
end

main()
```
![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/s94wav5t1zq4pvraws0m.png)

I'm quite proud of this little one :)
As you can see, as n grows larger, it looks like the Sierpinski triangle.

Speaking of Sierpinski triangle, and as we are in the *line* section, what about using Lindenmayer-System to do it?
```julia
using Lindenmayer
Sierpinski = LSystem(Dict("G" => "F+G+F", "F" => "G-F-G"), "F")

drawLSystem(Sierpinski,
    forward = 6,
    startingorientation = 0,
    turn = 60,
    iterations = 6,
    filename = "Sierpinsky.png")
```
![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/ch8eukesww848mssn0uz.png)

For Hilbert:
```julia
Hilbert = LSystem(Dict("X" => "-YF+XFX+FY-", "Y" => "+XF-YFY-FX+"), "X")
# rotation unit = -90
```
![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/0eggebpatqza35ner1j4.png)

and for some plants:
```julia
Plant = LSystem(Dict("F" => "FF+[+F-F-F]-[-F+F+F]"), "F")
# rotation unit = 22.5
```
![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/v6o7jdmdo18fyqpuk0l3.png)

Think about the possibilities!✨✨✨

You'll see lots of ways to draw plants when looking for fractals, so let's make another one, but this time, with random number and finally starting to scatter point-wise instead of lines (it allows for more flexibility).

```julia
using GLMakie
f1(x, y) = [0 0; 0 0.16] * [x, y]
f2(x, y) = [0.85 0.04; -0.04 0.85] * [x, y] + [0, 1.6]
f3(x, y) = [0.2 -0.26; 0.23 0.22] * [x, y] + [0, 1.6]
f4(x, y) = [-0.15 0.28; 0.26 0.24] * [x, y] + [0, 0.44]

function f(x, y)
	r = rand()
	r < 0.01 && return f1(x, y)
	r < 0.86 && return f2(x, y)
	r < 0.93 && return f3(x, y)
	f4(x, y)
end

function barnsley_fern(iter)
	X = [0.0]
	Y = [0.0]
	for i in 1:iter
		x, y = f(X[end], Y[end])
		push!(X, x)
		push!(Y, y)
	end
	scatter(X, Y, color=:green, markersize=1)
end

barnsley_fern(1_000_000)
```
![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/cxcci101krj6bvh20ygw.png)

Now that we are doing point-based fractal, we can do the mandelbrot and Julia's sets!

```julia
f(z, c) = z*z + c
function is_stable(iter, z, c)
	for _ in 1:iter
		if abs(z) > 2
			return false
		end
		z = f(z, c)
	end
	true
end

function mandel(precision, X, Y, f) # passing this f will be explained just after
    Points = Tuple{Float64, Float64}[]
	for x in X
		for y in Y
			z = f(x, y)
			if is_stable(precision, z, z)
				push!(Points, (x, y))
			end
		end
	end
	scatter(Points, markersize=1)
end


function julia(precision, X, Y, c, f)
	Points = Tuple{Float64, Float64}[]
	for x in X
		for y in Y
			z = f(x, y)
			if is_stable(precision, z, c)
				push!(Points, (x, y))
			end
		end
	end
	scatter(Points, markersize=4)
end

mandel(80, -2:0.0025:2, -2:0.0025:2, Complex) |> display
sleep(10)
julia(80, -2:0.0005:2, -2:0.0005:2, -0.8im, Complex) |> display
```

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/lvgkvh5i5vpl9faphrjk.png)
![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/p9qmlm5enclbfiiq4n59.png)

Julia set made in Julia, lovely.

Now... why did I pass the `f`? Welp, it's because you can decide to make your own to see how it will behave :)
This example is in 3D, so I didn't reuse some of the functions I made, but you will get the idea:
```julia
using CliffordAlgebras, GLMakie
const S = CliffordAlgebra(1, 1, 0)
const e1 = S.e1
const e2 = S.e2
# f definition
# is_stable definition
import Base.abs
abs(n::MultiVector) = vector(n) .^ 2 |> sum |> sqrt  # type piracy /(o_o)\

function mandel(precision, X, Y, Z, f)
	Points = Tuple{Float64, Float64, Float64}[]
	for x in X
		for y in Y
			for z in Z
				num = f(x, y, z)
				if is_stable(precision, num, num)
					push!(Points, (x, y, z))
				end
			end
		end
	end
	Points
end
const RANGE = -2:0.05:2
fun(a, b, c) = a + b*e1 + c*e2
mandel(50, RANGE, RANGE, RANGE, fun) |> (x -> scatter(x, markersize=4))
```
![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/suxmlx4wfo4flwur61jl.png)

Now you have a 3D slice of the mandelbrot set if defined in Cl(1, 1, 0).
And now, the last part, seeing this whole set thanks to an animation !

Let's suppose we have the same code until the `fun(a, b, c)`
```julia
n = Observable(0.0)
fn = lift(d -> (a, b, c) -> (a + b*e1 + c*e2 + d*e1*e2), n)
p = lift(f -> mandel(50, RANGE, RANGE, RANGE, f), fn)
fig, axe = scatter(p)
const frames = 1:28
record(fig, "Cl1-1-0_4DMandel.gif", frames; framerate=28) do frame
    n[] = -1.5+frame/10
end
```
![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/95blmc82td8ookv9ug75.gif)
See? with the help of passing `f`, you can do lots of things !
And so, here is the end, you're one of the lucky one who have seen the *whole* (except for all the missing points ...) mandelbrot set if defined in CliffordAlgebra(1, 1, 0) !
Hopes it was interesting.

So, why Julia?
You may not notice it because it's all compiling time, and TTFP is a real thing, but Julia is really fast, once I have my algo down, I can runs way more iterations than what I can in python.
It have some really great maths capacities, it's close to math in notation, and it feels quite natural.
It's extremely expressive, it just makes things easier to get something working well.
Awesome community, I really like some persons from the Human of Julia (discord) server.