vthriller/plan.png

## plan.png

      
    Raw
  

              plan.png
            
          
## wifi-field.jl
### A Pluto.jl notebook ###
# v0.12.17

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 1cf76d74-450a-11eb-3900-5dbbe000047a
begin
	using Colors
	using ColorTypes
	using Images
	using SparseArrays
	using DSP: conv2
	using Printf
	using PlutoUI
end

# ╔═╡ 14033bfc-4506-11eb-345d-b968dbdd669f
md"# Where's my Wi-Fi signal?"

# ╔═╡ b2375ede-4509-11eb-02c1-fde861ed36ad
md"## Solving the Helmhotz equation to find the good spots in your home"

# ╔═╡ b76f6990-4509-11eb-1f2c-8933290c6292
md"[*Helmhurts* by Jason Cole](http://jasmcole.com/2014/08/25/helmhurts/) explains how to apply the Maxwell equations, and in particular their smaller sibling the Helmholtz equations, to calculate the propagation of Wi-Fi signals."

# ╔═╡ bc937a5e-4509-11eb-0cfa-c5634086f187
md"What follows is an attempt to use Julia to go through the whole process described in the blog entry."

# ╔═╡ 1718edba-450a-11eb-2377-113207096427
md"### Setting the up the data"

# ╔═╡ 18908fcc-450a-11eb-18fb-53d6a09bd211
md"We'll start first by importing a floor plan, in that case a bitmap image that I produced with an image editor, but you can of course create your own image of your own place (or an imaginary one)."

# ╔═╡ 22808b72-450a-11eb-25ae-15d190b9c876
img = load("plan.png")

# ╔═╡ 3432951a-450a-11eb-0a87-05625cdca2f4
size(img)

# ╔═╡ 87240306-450f-11eb-2b4b-93419692a783
typeof(img)

# ╔═╡ 3531c2c2-450a-11eb-3e62-e31c5ad04546
md"Here we have a 2D-array of gray pixels.

For the purpose of the simulation we only need to know if a given point is air or concrete. Let's extract that material information from the pixels:"

# ╔═╡ 3d210916-450a-11eb-0ab4-4bde681a8fa2
plan = gray.(img);

# ╔═╡ 434cbaa6-450a-11eb-1e22-d50acd4eaabc
md"And set the necessary variables:"

# ╔═╡ 476333fe-450a-11eb-1b11-5bf8767c8847
begin
	dimx, dimy = size(plan)  # spatial dimensions
	fsize = length(plan)     # full size of the problem
	δ = 0.03                 # spatial resolution of each pixel
end;

# ╔═╡ 4efa869e-450a-11eb-222f-79c99c8ee078
md"Note that the spatial resolution is given by the actual size of the appartment divided by the number of pixel making the plan picture. Two important considerations come into play here:
1. Beware of not making the problem size (the number of pixels) too important, I found on my setup that you can go quickly from a problem that can be solved in seconds (100 x 100) to a never ending calculation or one that throws out of memory exceptions (1000 x 1000 in my case).
2. You can't on the other hand keep the pixel count low by increasing the spatial step δ (i.e. decreasing the resolution): when it gets too close to the wavelength (about 10 cm for Wi-Fi signals) the mesh becomes too coarse to give meaningful results.

Here, with a spatial resolution of 3 cm we should be safe for the resolution while still simulating a normally sized appartment of about $(@sprintf(\"%.2fm x %.2fm\", dimx * δ, dimy * δ))"

# ╔═╡ 55198d04-450a-11eb-0529-8f14d12e67cf
begin
	η_air = 1.                 # refraction index for air

	# refraction index for concrete
	# the imaginary part conveys the absorption
	# https://www.itu.int/dms_pubrec/itu-r/rec/p/R-REC-P.2040-1-201507-I!!PDF-E.pdf
	# §3 Compilations of electrical properties of materials
	η_concrete = begin
		local f = 2.4 # GHz
		local a, b, c, d = 5.31, 0, 0.0326, 0.8095
		local η_re = a * (b != 0 ? f ^ b : 1)
		local σ    = c * (d != 0 ? f ^ d : 1)
		local η_im = 17.98im * σ / f
		η_re + η_im
	end

	λ = 0.12                  # for a 2.5 GHz signal, wavelength is ~ 12cm
	k = 2π / λ                # k is the wavenumber
end;

# ╔═╡ 6af3e020-450a-11eb-09ea-392a4fe9d03e
md"We'll now create the Complex Array containing the material information:"

# ╔═╡ 73c7354e-450a-11eb-121d-0b694833a20b
begin
	μ = similar(plan, Complex)
	μ[plan .!= 0] .= (k / η_air)^2
	μ[plan .== 0] .= (k / η_concrete)^2
end;

# ╔═╡ 77c4f118-450a-11eb-336c-1d76291a980c
md"Note that we directly calculate the $\left( \frac{k}{n} \right) ^2$ value appearing in the Helmholtz equation below."

# ╔═╡ 7c03f8aa-450a-11eb-2006-83dfdb3c3527
md"### Solving the equation"

# ╔═╡ 85ae4306-450a-11eb-02dd-7d43dc23b9e0
md"The Helmholtz equation for the amplitude of the electric field is: $$\Delta A + \left( \frac{k}{n} \right) ^2 A= f$$, where $\Delta$ is the Laplacian operator, $A$ the amplitude field and $f$ the wave source that will be equal to zero everywhere except where the antenna is."

# ╔═╡ 8dbcf3a8-450a-11eb-10a5-1f0f89515bbb
md"The first step is to rephrase that problem in one of the ways (there are several) a computer can tackle it: by discretizing the space and using finite differences.

In that framework the Laplacian can be approximated by the neighboring values, giving us:
$$\frac{A(x+δ,y) - 2 A(x,y) + A(x-δ,y)}{δ^2}+\frac{A(x,y+δ) - 2 A(x,y) + A(x,y-δ)}{δ^2}+ \left( \frac{k}{n} \right) ^2 A(x,y) = f(x,y)$$"

# ╔═╡ ab2a9cc4-450a-11eb-23bf-1f194db3142e
md"We have here a linear combination of elements of A equal to elements of f. In others words an ordinary system of linear equations that computers are good at solving.

There is however a roadblock: that system has as many equations as there are values in A: $dimx x $dimy = $fsize. The matrix describing the problem would occupy an huge amount of memory if stored as a dense array: $fsize x $fsize x 8 x 2 (don't forget these are complex numbers) = $(@sprintf(\"%.1f\", fsize ^ 2 * 16 / 2 ^ 30)) GiB!

This memory would be mostly wasted though as the matrix is almost empty (only 5 elements appear in each equation not the whole $fsize). Sparses matrices are the adapted structure in these situations as they only store the non zero elements in memory. Julia provides a sparse matrix type on which most dense matrix methods can apply."

# ╔═╡ b0e9f36a-450a-11eb-1bcf-af7ce83a9a0c
md"Let's now build the three vectors describing the non zero elements of the sparse matrix `S` (the `LinearIndices` is used to relate the (x,y) positions in the real space to an index in the matrix):"

# ╔═╡ cbf3f270-450a-11eb-0eba-0be7069aafc4
begin
	xs = Array{Int}(undef, 5*dimx*dimy)
	ys = Array{Int}(undef, 5*dimx*dimy)
	vs = Array{Complex}(undef, 5*dimx*dimy)
	local i = 1
	local ind = LinearIndices(plan)
	for x in 1:dimx, y in 1:dimy  #  x=1; y=1
	    xm = (x+dimx-2) % dimx + 1
	    xp =          x % dimx + 1
	    ym = (y+dimy-2) % dimy + 1
	    yp =          y % dimy + 1

	    xs[i] = ind[x, y]; ys[i] = ind[ x,  y]; vs[i] = μ[x,y] - 2*δ^-2
	    i += 1
	    xs[i] = ind[x, y]; ys[i] = ind[xp,  y]; vs[i] = δ^-2
	    i += 1
	    xs[i] = ind[x, y]; ys[i] = ind[xm,  y]; vs[i] = δ^-2
	    i += 1
	    xs[i] = ind[x, y]; ys[i] = ind[ x, yp]; vs[i] = δ^-2
	    i += 1
	    xs[i] = ind[x, y]; ys[i] = ind[ x, ym]; vs[i] = δ^-2
	    i += 1
	end
end

# ╔═╡ ce7c4ef2-450a-11eb-0120-4dd2b5b6ebe2
md"We can now declare `S` by calling the `sparse()` function:"

# ╔═╡ d466118e-450a-11eb-0260-e931d1fe5edb
S = sparse(xs, ys, vs, fsize, fsize);

# ╔═╡ d76ed16a-450a-11eb-1ff6-13a784896b1f
md"There is only the emiter field $f$ left to define:"

# ╔═╡ 8a9acc22-4582-11eb-10d1-af27fb9fa40e
# TODO re-render this cell with each update of `img`
@bind antenna_pos html"""
	<span class="hack-input">
		<!--
		the first element is the one that's being watched,
		even it it's a text node,
		hence this span
		-->
		Click to pick antenna position:	<br>
		<canvas width="200" height="200" style="position: relative"></canvas>
	</span>
	<script>
		const emitter = currentScript.closest('pluto-output').querySelector(".hack-input")
		const canvas = emitter.querySelector("canvas")
		const ctx = canvas.getContext("2d")

		canvas.clear = () => {
			// XXX UGLY HACK
			// I was not able to figure out how to interpolate `img` into HTML literals,
			// and HTML in Markdown strings is not supported,
			// so instead I clone already existing image from the cell that defined the `img`
			let img = document.getElementById("22808b72-450a-11eb-25ae-15d190b9c876").querySelector('pluto-output img')

			canvas.width = img.width
			canvas.height = img.height
			ctx.drawImage(img, 0, 0)
		}

		canvas.onclick = e => {
			// send the value to Pluto
			// N.B. coords are getting swapped:
			// JS: x→ y↓
			// Julia: x↓ y→
			emitter.value = [e.layerY, e.layerX]
			emitter.dispatchEvent(new CustomEvent("input"))

			canvas.clear();

			// draw the crosshair
			ctx.strokeStyle = '#c00'
			ctx.beginPath()
			ctx.moveTo(e.layerX - 3, e.layerY)
			ctx.lineTo(e.layerX + 3, e.layerY)
			ctx.stroke()
			ctx.beginPath()
			ctx.moveTo(e.layerX, e.layerY - 3)
			ctx.lineTo(e.layerX, e.layerY + 3)
			ctx.stroke()
		}

		// show the image for the first time
		canvas.clear()
	</script>
"""

# ╔═╡ a6105b1a-45a1-11eb-37bc-25e709d2aa61
antenna_pos

# ╔═╡ dd620024-450a-11eb-0d26-39444e86dbc7
begin
	f = fill(0 + 0im, (dimx, dimy))
	f[antenna_pos...] = 1.0 # our Wi-Fi emitter antenna will be here
end;

# ╔═╡ e13d4758-450a-11eb-21fa-21ef2b053b3c
md"Now, with the above notations, we have the linear system $S.A=f$ to solve for A:"

# ╔═╡ e53559e8-450a-11eb-10df-2f9bd44c829b
A = reshape(S \ vec(f), dimx, dimy);

# ╔═╡ ea215d76-450a-11eb-34c4-8f840a4c07e0
md"A contains the amplitude field we are looking for in the real part of the complex number. Since the relevant physical value is rather the Energy received not the amplitude, we'll calculate the square of the amplitude."

# ╔═╡ ed05c25c-450a-11eb-3ca4-fd98065261f4
E_raw = real(A) .* real(A);

# ╔═╡ f13d360c-450a-11eb-3ea5-a3c6738d3f7b
md"And since the receiving antennas of laptops and smartphones have a size of at least 10cm, we'll convolve the energy field to have a more accurate representation of the received signal strength"

# ╔═╡ 8faa5ac4-4574-11eb-113f-49565b459121
@bind E_convolve html"<input type='checkbox' checked> convolve energy field"

# ╔═╡ f4454b66-450a-11eb-1429-f3e5cd06ae3b
E = E_convolve ? conv2(E_raw, ones(5,5)/25 )[3:end-2, 3:end-2] : E_raw;

# ╔═╡ f7acc81a-450a-11eb-3c0c-61badd6b8835
md"### Plotting the result"

# ╔═╡ faf621ce-450a-11eb-0158-398569973267
md"We can now convert the E field back to an image to visualize where the good spots are for Wi-Fi reception in our appartment.

First we'll translate the energy field to 1-100 interger scale:"

# ╔═╡ ff2a96bc-450a-11eb-1d55-c1a1aab6dfc4
Ei = min.(100, max.(1, Int.(trunc.( 1 .+ 100 .* (E .- minimum(E)) / (maximum(E) - minimum(E)) ) )));

# ╔═╡ 039ef896-450b-11eb-231d-e1a15b86df0c
md"Then we'll create a color scale using the `colormap()` function of the `Colors` package:"

# ╔═╡ 07031f76-450b-11eb-1824-8367aa6aca36
cm = reverse(colormap("oranges"));

# ╔═╡ 14232eda-450b-11eb-23a1-9731174997ad
md"The last step is to create a color image by picking the RGB components of the color scale:"

# ╔═╡ 02251d32-4575-11eb-31de-cb4d694dd5b3
@bind overlay_walls html"<input type='checkbox' checked> show walls"

# ╔═╡ 176b18dc-450b-11eb-20bd-c3ab5982c5b6
begin
	field = similar(plan, RGB{Float64})
	local ind = CartesianIndices(plan)
	for i in 1:fsize
		field[ind[i]] = cm[Ei[i]]
	end

	if overlay_walls
		field[plan .== 0] .= RGB(0)
	end

	# show antenna position with a green crosshair
	# first, the outline to make it more visible on top of the bright background
	local ax, ay = antenna_pos
	for d in -3:3
		field[ax + d, ay - 1] = RGB(0, .2, 0)
		field[ax + d, ay + 1] = RGB(0, .2, 0)
		field[ax - 1, ay + d] = RGB(0, .2, 0)
		field[ax + 1, ay + d] = RGB(0, .2, 0)
	end
	# next, the crosshair itself
	for d in -3:3
		field[ax + d, ay] = RGB(0, 1, 0)
		field[ax, ay + d] = RGB(0, 1, 0)
	end
	field
end

# ╔═╡ e85b093c-4500-11eb-3441-c1d7431c634e
md"## et voilà!"

# ╔═╡ Cell order:
# ╟─14033bfc-4506-11eb-345d-b968dbdd669f
# ╟─b2375ede-4509-11eb-02c1-fde861ed36ad
# ╟─b76f6990-4509-11eb-1f2c-8933290c6292
# ╟─bc937a5e-4509-11eb-0cfa-c5634086f187
# ╟─1718edba-450a-11eb-2377-113207096427
# ╟─18908fcc-450a-11eb-18fb-53d6a09bd211
# ╠═1cf76d74-450a-11eb-3900-5dbbe000047a
# ╠═22808b72-450a-11eb-25ae-15d190b9c876
# ╠═3432951a-450a-11eb-0a87-05625cdca2f4
# ╠═87240306-450f-11eb-2b4b-93419692a783
# ╟─3531c2c2-450a-11eb-3e62-e31c5ad04546
# ╠═3d210916-450a-11eb-0ab4-4bde681a8fa2
# ╟─434cbaa6-450a-11eb-1e22-d50acd4eaabc
# ╠═476333fe-450a-11eb-1b11-5bf8767c8847
# ╟─4efa869e-450a-11eb-222f-79c99c8ee078
# ╠═55198d04-450a-11eb-0529-8f14d12e67cf
# ╟─6af3e020-450a-11eb-09ea-392a4fe9d03e
# ╠═73c7354e-450a-11eb-121d-0b694833a20b
# ╟─77c4f118-450a-11eb-336c-1d76291a980c
# ╟─7c03f8aa-450a-11eb-2006-83dfdb3c3527
# ╟─85ae4306-450a-11eb-02dd-7d43dc23b9e0
# ╟─8dbcf3a8-450a-11eb-10a5-1f0f89515bbb
# ╟─ab2a9cc4-450a-11eb-23bf-1f194db3142e
# ╟─b0e9f36a-450a-11eb-1bcf-af7ce83a9a0c
# ╠═cbf3f270-450a-11eb-0eba-0be7069aafc4
# ╟─ce7c4ef2-450a-11eb-0120-4dd2b5b6ebe2
# ╠═d466118e-450a-11eb-0260-e931d1fe5edb
# ╟─d76ed16a-450a-11eb-1ff6-13a784896b1f
# ╟─8a9acc22-4582-11eb-10d1-af27fb9fa40e
# ╠═a6105b1a-45a1-11eb-37bc-25e709d2aa61
# ╠═dd620024-450a-11eb-0d26-39444e86dbc7
# ╟─e13d4758-450a-11eb-21fa-21ef2b053b3c
# ╠═e53559e8-450a-11eb-10df-2f9bd44c829b
# ╟─ea215d76-450a-11eb-34c4-8f840a4c07e0
# ╠═ed05c25c-450a-11eb-3ca4-fd98065261f4
# ╟─f13d360c-450a-11eb-3ea5-a3c6738d3f7b
# ╟─8faa5ac4-4574-11eb-113f-49565b459121
# ╠═f4454b66-450a-11eb-1429-f3e5cd06ae3b
# ╟─f7acc81a-450a-11eb-3c0c-61badd6b8835
# ╟─faf621ce-450a-11eb-0158-398569973267
# ╠═ff2a96bc-450a-11eb-1d55-c1a1aab6dfc4
# ╟─039ef896-450b-11eb-231d-e1a15b86df0c
# ╠═07031f76-450b-11eb-1824-8367aa6aca36
# ╟─14232eda-450b-11eb-23a1-9731174997ad
# ╟─02251d32-4575-11eb-31de-cb4d694dd5b3
# ╠═176b18dc-450b-11eb-20bd-c3ab5982c5b6
# ╟─e85b093c-4500-11eb-3441-c1d7431c634e
	### A Pluto.jl notebook ###
	# v0.12.17

	using Markdown
	using InteractiveUtils

	# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
	macro bind(def, element)
	quote
	local el = $(esc(element))
	global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
	el
	end
	end

	# ╔═╡ 1cf76d74-450a-11eb-3900-5dbbe000047a
	begin
	using Colors
	using ColorTypes
	using Images
	using SparseArrays
	using DSP: conv2
	using Printf
	using PlutoUI
	end

	# ╔═╡ 14033bfc-4506-11eb-345d-b968dbdd669f
	md"# Where's my Wi-Fi signal?"

	# ╔═╡ b2375ede-4509-11eb-02c1-fde861ed36ad
	md"## Solving the Helmhotz equation to find the good spots in your home"

	# ╔═╡ b76f6990-4509-11eb-1f2c-8933290c6292
	md"[Helmhurts by Jason Cole](http://jasmcole.com/2014/08/25/helmhurts/) explains how to apply the Maxwell equations, and in particular their smaller sibling the Helmholtz equations, to calculate the propagation of Wi-Fi signals."

	# ╔═╡ bc937a5e-4509-11eb-0cfa-c5634086f187
	md"What follows is an attempt to use Julia to go through the whole process described in the blog entry."

	# ╔═╡ 1718edba-450a-11eb-2377-113207096427
	md"### Setting the up the data"

	# ╔═╡ 18908fcc-450a-11eb-18fb-53d6a09bd211
	md"We'll start first by importing a floor plan, in that case a bitmap image that I produced with an image editor, but you can of course create your own image of your own place (or an imaginary one)."

	# ╔═╡ 22808b72-450a-11eb-25ae-15d190b9c876
	img = load("plan.png")

	# ╔═╡ 3432951a-450a-11eb-0a87-05625cdca2f4
	size(img)

	# ╔═╡ 87240306-450f-11eb-2b4b-93419692a783
	typeof(img)

	# ╔═╡ 3531c2c2-450a-11eb-3e62-e31c5ad04546
	md"Here we have a 2D-array of gray pixels.

	For the purpose of the simulation we only need to know if a given point is air or concrete. Let's extract that material information from the pixels:"

	# ╔═╡ 3d210916-450a-11eb-0ab4-4bde681a8fa2
	plan = gray.(img);

	# ╔═╡ 434cbaa6-450a-11eb-1e22-d50acd4eaabc
	md"And set the necessary variables:"

	# ╔═╡ 476333fe-450a-11eb-1b11-5bf8767c8847
	begin
	dimx, dimy = size(plan) # spatial dimensions
	fsize = length(plan) # full size of the problem
	δ = 0.03 # spatial resolution of each pixel
	end;

	# ╔═╡ 4efa869e-450a-11eb-222f-79c99c8ee078
	md"Note that the spatial resolution is given by the actual size of the appartment divided by the number of pixel making the plan picture. Two important considerations come into play here:
	1. Beware of not making the problem size (the number of pixels) too important, I found on my setup that you can go quickly from a problem that can be solved in seconds (100 x 100) to a never ending calculation or one that throws out of memory exceptions (1000 x 1000 in my case).
	2. You can't on the other hand keep the pixel count low by increasing the spatial step δ (i.e. decreasing the resolution): when it gets too close to the wavelength (about 10 cm for Wi-Fi signals) the mesh becomes too coarse to give meaningful results.

	Here, with a spatial resolution of 3 cm we should be safe for the resolution while still simulating a normally sized appartment of about $(@sprintf(\"%.2fm x %.2fm\", dimx * δ, dimy * δ))"

	# ╔═╡ 55198d04-450a-11eb-0529-8f14d12e67cf
	begin
	η_air = 1. # refraction index for air

	# refraction index for concrete
	# the imaginary part conveys the absorption
	# https://www.itu.int/dms_pubrec/itu-r/rec/p/R-REC-P.2040-1-201507-I!!PDF-E.pdf
	# §3 Compilations of electrical properties of materials
	η_concrete = begin
	local f = 2.4 # GHz
	local a, b, c, d = 5.31, 0, 0.0326, 0.8095
	local η_re = a * (b != 0 ? f ^ b : 1)
	local σ = c * (d != 0 ? f ^ d : 1)
	local η_im = 17.98im * σ / f
	η_re + η_im
	end

	λ = 0.12 # for a 2.5 GHz signal, wavelength is ~ 12cm
	k = 2π / λ # k is the wavenumber
	end;

	# ╔═╡ 6af3e020-450a-11eb-09ea-392a4fe9d03e
	md"We'll now create the Complex Array containing the material information:"

	# ╔═╡ 73c7354e-450a-11eb-121d-0b694833a20b
	begin
	μ = similar(plan, Complex)
	μ[plan .!= 0] .= (k / η_air)^2
	μ[plan .== 0] .= (k / η_concrete)^2
	end;

	# ╔═╡ 77c4f118-450a-11eb-336c-1d76291a980c
	md"Note that we directly calculate the $\left( \frac{k}{n} \right) ^2$ value appearing in the Helmholtz equation below."

	# ╔═╡ 7c03f8aa-450a-11eb-2006-83dfdb3c3527
	md"### Solving the equation"

	# ╔═╡ 85ae4306-450a-11eb-02dd-7d43dc23b9e0
	md"The Helmholtz equation for the amplitude of the electric field is: $$\Delta A + \left( \frac{k}{n} \right) ^2 A= f$$, where $\Delta$ is the Laplacian operator, $A$ the amplitude field and $f$ the wave source that will be equal to zero everywhere except where the antenna is."

	# ╔═╡ 8dbcf3a8-450a-11eb-10a5-1f0f89515bbb
	md"The first step is to rephrase that problem in one of the ways (there are several) a computer can tackle it: by discretizing the space and using finite differences.

	In that framework the Laplacian can be approximated by the neighboring values, giving us:
	$$\frac{A(x+δ,y) - 2 A(x,y) + A(x-δ,y)}{δ^2}+\frac{A(x,y+δ) - 2 A(x,y) + A(x,y-δ)}{δ^2}+ \left( \frac{k}{n} \right) ^2 A(x,y) = f(x,y)$$"

	# ╔═╡ ab2a9cc4-450a-11eb-23bf-1f194db3142e
	md"We have here a linear combination of elements of A equal to elements of f. In others words an ordinary system of linear equations that computers are good at solving.

	There is however a roadblock: that system has as many equations as there are values in A: $dimx x $dimy = $fsize. The matrix describing the problem would occupy an huge amount of memory if stored as a dense array: $fsize x $fsize x 8 x 2 (don't forget these are complex numbers) = $(@sprintf(\"%.1f\", fsize ^ 2 * 16 / 2 ^ 30)) GiB!

	This memory would be mostly wasted though as the matrix is almost empty (only 5 elements appear in each equation not the whole $fsize). Sparses matrices are the adapted structure in these situations as they only store the non zero elements in memory. Julia provides a sparse matrix type on which most dense matrix methods can apply."

	# ╔═╡ b0e9f36a-450a-11eb-1bcf-af7ce83a9a0c
	md"Let's now build the three vectors describing the non zero elements of the sparse matrix `S` (the `LinearIndices` is used to relate the (x,y) positions in the real space to an index in the matrix):"

	# ╔═╡ cbf3f270-450a-11eb-0eba-0be7069aafc4
	begin
	xs = Array{Int}(undef, 5dimxdimy)
	ys = Array{Int}(undef, 5dimxdimy)
	vs = Array{Complex}(undef, 5dimxdimy)
	local i = 1
	local ind = LinearIndices(plan)
	for x in 1:dimx, y in 1:dimy # x=1; y=1
	xm = (x+dimx-2) % dimx + 1
	xp = x % dimx + 1
	ym = (y+dimy-2) % dimy + 1
	yp = y % dimy + 1

	xs[i] = ind[x, y]; ys[i] = ind[ x, y]; vs[i] = μ[x,y] - 2*δ^-2
	i += 1
	xs[i] = ind[x, y]; ys[i] = ind[xp, y]; vs[i] = δ^-2
	i += 1
	xs[i] = ind[x, y]; ys[i] = ind[xm, y]; vs[i] = δ^-2
	i += 1
	xs[i] = ind[x, y]; ys[i] = ind[ x, yp]; vs[i] = δ^-2
	i += 1
	xs[i] = ind[x, y]; ys[i] = ind[ x, ym]; vs[i] = δ^-2
	i += 1
	end
	end

	# ╔═╡ ce7c4ef2-450a-11eb-0120-4dd2b5b6ebe2
	md"We can now declare `S` by calling the `sparse()` function:"

	# ╔═╡ d466118e-450a-11eb-0260-e931d1fe5edb
	S = sparse(xs, ys, vs, fsize, fsize);

	# ╔═╡ d76ed16a-450a-11eb-1ff6-13a784896b1f
	md"There is only the emiter field $f$ left to define:"

	# ╔═╡ 8a9acc22-4582-11eb-10d1-af27fb9fa40e
	# TODO re-render this cell with each update of `img`
	@bind antenna_pos html"""
	<span class="hack-input">
	<!--
	the first element is the one that's being watched,
	even it it's a text node,
	hence this span
	-->
	Click to pick antenna position: <br>
	<canvas width="200" height="200" style="position: relative"></canvas>
	</span>
	<script>
	const emitter = currentScript.closest('pluto-output').querySelector(".hack-input")
	const canvas = emitter.querySelector("canvas")
	const ctx = canvas.getContext("2d")

	canvas.clear = () => {
	// XXX UGLY HACK
	// I was not able to figure out how to interpolate `img` into HTML literals,
	// and HTML in Markdown strings is not supported,
	// so instead I clone already existing image from the cell that defined the `img`
	let img = document.getElementById("22808b72-450a-11eb-25ae-15d190b9c876").querySelector('pluto-output img')

	canvas.width = img.width
	canvas.height = img.height
	ctx.drawImage(img, 0, 0)
	}

	canvas.onclick = e => {
	// send the value to Pluto
	// N.B. coords are getting swapped:
	// JS: x→ y↓
	// Julia: x↓ y→
	emitter.value = [e.layerY, e.layerX]
	emitter.dispatchEvent(new CustomEvent("input"))

	canvas.clear();

	// draw the crosshair
	ctx.strokeStyle = '#c00'
	ctx.beginPath()
	ctx.moveTo(e.layerX - 3, e.layerY)
	ctx.lineTo(e.layerX + 3, e.layerY)
	ctx.stroke()
	ctx.beginPath()
	ctx.moveTo(e.layerX, e.layerY - 3)
	ctx.lineTo(e.layerX, e.layerY + 3)
	ctx.stroke()
	}

	// show the image for the first time
	canvas.clear()
	</script>
	"""

	# ╔═╡ a6105b1a-45a1-11eb-37bc-25e709d2aa61
	antenna_pos

	# ╔═╡ dd620024-450a-11eb-0d26-39444e86dbc7
	begin
	f = fill(0 + 0im, (dimx, dimy))
	f[antenna_pos...] = 1.0 # our Wi-Fi emitter antenna will be here
	end;

	# ╔═╡ e13d4758-450a-11eb-21fa-21ef2b053b3c
	md"Now, with the above notations, we have the linear system $S.A=f$ to solve for A:"

	# ╔═╡ e53559e8-450a-11eb-10df-2f9bd44c829b
	A = reshape(S \ vec(f), dimx, dimy);

	# ╔═╡ ea215d76-450a-11eb-34c4-8f840a4c07e0
	md"A contains the amplitude field we are looking for in the real part of the complex number. Since the relevant physical value is rather the Energy received not the amplitude, we'll calculate the square of the amplitude."

	# ╔═╡ ed05c25c-450a-11eb-3ca4-fd98065261f4
	E_raw = real(A) .* real(A);

	# ╔═╡ f13d360c-450a-11eb-3ea5-a3c6738d3f7b
	md"And since the receiving antennas of laptops and smartphones have a size of at least 10cm, we'll convolve the energy field to have a more accurate representation of the received signal strength"

	# ╔═╡ 8faa5ac4-4574-11eb-113f-49565b459121
	@bind E_convolve html"<input type='checkbox' checked> convolve energy field"

	# ╔═╡ f4454b66-450a-11eb-1429-f3e5cd06ae3b
	E = E_convolve ? conv2(E_raw, ones(5,5)/25 )[3:end-2, 3:end-2] : E_raw;

	# ╔═╡ f7acc81a-450a-11eb-3c0c-61badd6b8835
	md"### Plotting the result"

	# ╔═╡ faf621ce-450a-11eb-0158-398569973267
	md"We can now convert the E field back to an image to visualize where the good spots are for Wi-Fi reception in our appartment.

	First we'll translate the energy field to 1-100 interger scale:"

	# ╔═╡ ff2a96bc-450a-11eb-1d55-c1a1aab6dfc4
	Ei = min.(100, max.(1, Int.(trunc.( 1 .+ 100 .* (E .- minimum(E)) / (maximum(E) - minimum(E)) ) )));

	# ╔═╡ 039ef896-450b-11eb-231d-e1a15b86df0c
	md"Then we'll create a color scale using the `colormap()` function of the `Colors` package:"

	# ╔═╡ 07031f76-450b-11eb-1824-8367aa6aca36
	cm = reverse(colormap("oranges"));

	# ╔═╡ 14232eda-450b-11eb-23a1-9731174997ad
	md"The last step is to create a color image by picking the RGB components of the color scale:"

	# ╔═╡ 02251d32-4575-11eb-31de-cb4d694dd5b3
	@bind overlay_walls html"<input type='checkbox' checked> show walls"

	# ╔═╡ 176b18dc-450b-11eb-20bd-c3ab5982c5b6
	begin
	field = similar(plan, RGB{Float64})
	local ind = CartesianIndices(plan)
	for i in 1:fsize
	field[ind[i]] = cm[Ei[i]]
	end

	if overlay_walls
	field[plan .== 0] .= RGB(0)
	end

	# show antenna position with a green crosshair
	# first, the outline to make it more visible on top of the bright background
	local ax, ay = antenna_pos
	for d in -3:3
	field[ax + d, ay - 1] = RGB(0, .2, 0)
	field[ax + d, ay + 1] = RGB(0, .2, 0)
	field[ax - 1, ay + d] = RGB(0, .2, 0)
	field[ax + 1, ay + d] = RGB(0, .2, 0)
	end
	# next, the crosshair itself
	for d in -3:3
	field[ax + d, ay] = RGB(0, 1, 0)
	field[ax, ay + d] = RGB(0, 1, 0)
	end
	field
	end

	# ╔═╡ e85b093c-4500-11eb-3441-c1d7431c634e
	md"## et voilà!"

	# ╔═╡ Cell order:
	# ╟─14033bfc-4506-11eb-345d-b968dbdd669f
	# ╟─b2375ede-4509-11eb-02c1-fde861ed36ad
	# ╟─b76f6990-4509-11eb-1f2c-8933290c6292
	# ╟─bc937a5e-4509-11eb-0cfa-c5634086f187
	# ╟─1718edba-450a-11eb-2377-113207096427
	# ╟─18908fcc-450a-11eb-18fb-53d6a09bd211
	# ╠═1cf76d74-450a-11eb-3900-5dbbe000047a
	# ╠═22808b72-450a-11eb-25ae-15d190b9c876
	# ╠═3432951a-450a-11eb-0a87-05625cdca2f4
	# ╠═87240306-450f-11eb-2b4b-93419692a783
	# ╟─3531c2c2-450a-11eb-3e62-e31c5ad04546
	# ╠═3d210916-450a-11eb-0ab4-4bde681a8fa2
	# ╟─434cbaa6-450a-11eb-1e22-d50acd4eaabc
	# ╠═476333fe-450a-11eb-1b11-5bf8767c8847
	# ╟─4efa869e-450a-11eb-222f-79c99c8ee078
	# ╠═55198d04-450a-11eb-0529-8f14d12e67cf
	# ╟─6af3e020-450a-11eb-09ea-392a4fe9d03e
	# ╠═73c7354e-450a-11eb-121d-0b694833a20b
	# ╟─77c4f118-450a-11eb-336c-1d76291a980c
	# ╟─7c03f8aa-450a-11eb-2006-83dfdb3c3527
	# ╟─85ae4306-450a-11eb-02dd-7d43dc23b9e0
	# ╟─8dbcf3a8-450a-11eb-10a5-1f0f89515bbb
	# ╟─ab2a9cc4-450a-11eb-23bf-1f194db3142e
	# ╟─b0e9f36a-450a-11eb-1bcf-af7ce83a9a0c
	# ╠═cbf3f270-450a-11eb-0eba-0be7069aafc4
	# ╟─ce7c4ef2-450a-11eb-0120-4dd2b5b6ebe2
	# ╠═d466118e-450a-11eb-0260-e931d1fe5edb
	# ╟─d76ed16a-450a-11eb-1ff6-13a784896b1f
	# ╟─8a9acc22-4582-11eb-10d1-af27fb9fa40e
	# ╠═a6105b1a-45a1-11eb-37bc-25e709d2aa61
	# ╠═dd620024-450a-11eb-0d26-39444e86dbc7
	# ╟─e13d4758-450a-11eb-21fa-21ef2b053b3c
	# ╠═e53559e8-450a-11eb-10df-2f9bd44c829b
	# ╟─ea215d76-450a-11eb-34c4-8f840a4c07e0
	# ╠═ed05c25c-450a-11eb-3ca4-fd98065261f4
	# ╟─f13d360c-450a-11eb-3ea5-a3c6738d3f7b
	# ╟─8faa5ac4-4574-11eb-113f-49565b459121
	# ╠═f4454b66-450a-11eb-1429-f3e5cd06ae3b
	# ╟─f7acc81a-450a-11eb-3c0c-61badd6b8835
	# ╟─faf621ce-450a-11eb-0158-398569973267
	# ╠═ff2a96bc-450a-11eb-1d55-c1a1aab6dfc4
	# ╟─039ef896-450b-11eb-231d-e1a15b86df0c
	# ╠═07031f76-450b-11eb-1824-8367aa6aca36
	# ╟─14232eda-450b-11eb-23a1-9731174997ad
	# ╟─02251d32-4575-11eb-31de-cb4d694dd5b3
	# ╠═176b18dc-450b-11eb-20bd-c3ab5982c5b6
	# ╟─e85b093c-4500-11eb-3441-c1d7431c634e