34.2 Solution

Let’s start with the container and its functionality.

const N_COLS = 80
const N_ROWS = 40

function addBorders!(container::Matrix{Char})::Nothing
    container[:, 1] .= '|'
    container[:, N_COLS] .= '|'
    container[1, :] .= '-'
    container[N_ROWS, :] .= '-'
    return nothing
end

function getEmptyContainer()::Matrix{Char}
    container::Matrix{Char} = fill(' ', N_ROWS, N_COLS);
    addBorders!(container)
    return container
end

function printContainer(container::Matrix{Char})::Nothing
    for r in 1:N_ROWS
        println(container[r, :] |> join)
    end
    return nothing
end

function clearDisplay(nLinesUp::Int)::Nothing
    @assert 0 < nLinesUp "nLinesUp must be a positive integer"
    # "\033[xxxA" - xxx moves cursor up xxx LINES
    print("\033[" * string(nLinesUp) * "A")
    # "\033[0J" - clears from cursor position till the end of the screen
    print("\033[0J")
    return nothing
end

The container is just a Matrix (a table) of characters (Chars) that’s constrained by the borders (| and -) and initially contains nothing inside (fill(' ', N_ROWS, N_COLS)).

Now we need our molecules. The above will be defined as a vector of positions denoting their locations (row and column) within the matrix (our container).

const Pos = Tuple{Int, Int} # position (row, col) in 2D container

function isWithinContainer(molecule::Pos)::Bool
    row, col = molecule
    # accounts for borders
    return (1 < row < N_ROWS) && (1 < col < N_COLS)
end

At first molecules will occupy random positions.

# assumption: molecules may pass through each other
# (or occupy the same pixel in 2D) since they move
# past each other in the third (not drawn) dimension
function placeMoleculesRandomly!(molecules::Vec{Pos},
                                 rowMin::Int, rowMax::Int,
                                 colMin::Int, colMax::Int)::Nothing
    @assert(isWithinContainer((rowMin, colMin)),
            "(rowMin, colMin) outside of container")
    @assert(isWithinContainer((rowMax, colMax)),
            "(rowMax, colMax) outside of container")
    r::Int, c::Int = 0, 0
    for i in eachindex(molecules)
        r = rand(rowMin:rowMax)
        c = rand(colMin:colMax)
        molecules[i] = (r, c)
    end
    return nothing
end

const MOLECULE = '.'

function addMolecules2container!(molecules::Vec{Pos},
                                 container!::Matrix{Char})::Nothing
    for molecule in molecules
        if isWithinContainer(molecule)
            container![molecule...] = MOLECULE
        end
    end
    return nothing
end

Using rowMin/rowMax, colMin/colMax allows to place the molecules only in some part of the matrix (per task specification it will be the left side of container). Notice ! character in addMolecules2container!. Per Julia’s convention it was added to the name of the function that modifies its contents. However, both molecules (Vec{Pos}) and container (Matrix{Char}) are passed by reference. So a question may arise which one of the two (or maybe both) will get modified. To help with the answer the second parameter was named container! to emphasize that only it will be modified by the function. On the other hand, placeMoleculesRandomly! may modify at most one of its arguments (molecules) so there is no need for an extra ! which might be confusing at first glance. Anyway, the key thing is that based on the positions (molecules) we placed a marker MOLECULE = '.' in our container which later on will be displayed to the user.

Time to implement Brownian motion or:

[…] a normal distribution with the mean \(\mu = 0\) and variance \(\sigma^{2} = 2Dt\) usually called Brownian motion \(B_{t}\) […]

(quote from the Wiki link provided above)

Here, we are OK with all the molecules being identical (and sharing identical properties). Moreover, for simplicity we’ll just assume that D = 0.5 and t = 1, hence 2Dt = 1. This last action will give us a normal distribution with the mean \(\mu = 0\) and variance \(\sigma^{2} = 1\) (standard deviation sd is also 1, since \(sd = \sqrt{variance}\)). Luckily, that is what the built-in randn function provides. Hence, we will calculate the new position of a molecule to be new_position = old_position + shift (randn() provides the shift) and implement it like so:

round2int(f::Flt)::Int = round(Int, f)

function getNewPosition(molecule::Pos)::Pos
    rowShift::Int = randn() |> round2int
    colShift::Int = randn() |> round2int
    return molecule .+ (rowShift, colShift)
end

const N_MOLECULES = 150

# assumption: molecules may pass through each other
# (or occupy the same pixel in 2D) since they move
# past each other in the third (not drawn) dimension
function make1BrownianCycleShift!(molecules::Vec{Pos})::Nothing
    i::Int = 1
    newPos::Pos = (0, 0)
    while i <= N_MOLECULES
        newPos = getNewPosition(molecules[i])
        if isWithinContainer(newPos)
            molecules[i] = newPos
            i += 1
        end
    end
    return nothing
end

The terminal display will require to give a shift in integers, hence round2int implemented as single expression function syntax. Moreover, we cannot allow a particle to go outside the walls of the container, thus the if isWithinContainer(newPos), etc. statement. Together with the surrounding while block it makes sure that the molecule ‘falls back’ to the container. Effectively, it kind of simulates a reflection of a molecule from the walls of the vessel in a random direction.

Now we are almost ready for running our simulation, but first two, rather self-explanatory, functions:

const N_CYCLES = 2_500

getCol((_, c)::Pos)::Int = c

function getRightLeftCountsInfo(molecules::Vec{Pos})::Str
    colsWithMolecules::Vec{Int} = map(getCol, molecules)
    midCol::Int = round2int(N_COLS/2)
    lCount::Int = sum(colsWithMolecules .<= midCol)
    rCount::Int = sum(colsWithMolecules .> midCol)
    return "Left count: $lCount | Right count: $rCount"
end

function redrawDisplay(container::Matrix{Char},
                       molecules::Vec{Pos},
                       nCycle::Int)::Nothing
    clearDisplay(N_ROWS+2) # container + 2 info lines below
    println("Cycle no: $nCycle/$N_CYCLES")
    println(getRightLeftCountsInfo(molecules))
    printContainer(container)
    return nothing
end

Finally, crème de la crème, the simulation itself.

const DELAY_SEC = 0.1

function simulateBrownianMotions(nCycles::Int=N_CYCLES)::Nothing
    @assert 500 <= nCycles <= 1e5 "nCycles must be in range [500, 1e5]"

    container::Matrix{Char} = getEmptyContainer()
    molecules::Vec{Pos} = fill((0, 0), N_MOLECULES)
    placeMoleculesRandomly!(molecules,
                            # adjusted for borders
                            2, N_ROWS-1,
                            2, round2int((N_COLS-2)/2)
                            )
    addMolecules2container!(molecules, container)
    redrawDisplay(container, molecules, 0)

    for cycleNumber in 1:nCycles
        make1BrownianCycleShift!(molecules)
        emptyContainer!(container)
        addMolecules2container!(molecules, container)
        sleep(DELAY_SEC)
        redrawDisplay(container, molecules, cycleNumber)
    end

    return nothing
end

Nothing special here, we just declare container/molecules and initialize them with the appropriate values. Then, for each cycle 1:nCycles we make1BrownianCycleShift, remove the old molecule symbols from the container (emptyContainer), add the symbols of new, shifted molecules (addMolecules2container), pause for a moment (sleep) and redraw everything (redrawDisplay).

All that’s left to do is to add the main function since the app is meant to be run from terminal.

const SECS_PER_MIN = 60
const DURATION_SEC = DELAY_SEC * N_CYCLES
const DURATION_MIN = DURATION_SEC / SECS_PER_MIN

rnd2(x::Flt)::Flt = round(x, digits=2)

function main()::Nothing
    println("\nThis is a toy program that models simplified diffusion.")
    println("Note: your terminal must support ANSI escape codes.\n")

    println("Estimated execution time of the program:")
    print("$(rnd2(DURATION_SEC)) seconds or $(rnd2(DURATION_MIN)) ")
    println("minutes.")
    print("WARNING: the screen may flicker ")
    println("(Ctrl-C should abort the program).")

    # y(es) - default choice (also with Enter), anything else: no
    println("\nContinue with the simulation? [Y/n]")
    choice::Str = readline()
    if lowercase(strip(choice)) in ["y", "yes", ""]
        simulateBrownianMotions()
    end

    println("\nThat's all. Goodbye!")

    return nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

The end result is to be seen below (Figure 20).

Amazing. So the diffusion works in the room you find yourself in and even in our oversimplified simulation. All it took was a mock-up of Brownian motion, a reflection from a wall of the container and a few thousand cycles of random shifts to properly mix the particles.