Simulation of a Scanning Tunneling Microscope's control algorithm (Pluto notebook)

STM.jl

### A Pluto.jl notebook ###
# v0.12.21
# (c) 2021 Lewin Bormann

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

# ╔═╡ 3ca61268-8a5d-11eb-34ab-9b0c7ce9c19a
begin
	using Plots
	using Statistics
	import Random
	
	gr()
end

# ╔═╡ 662c9a8a-8a5d-11eb-1ca0-972e75cce9de
md"""
## PID-Regelung Rastertunnelmikroskop

Wir haben
* ein Stromsignal $I(t)$,
* ein Fehlersignal $E(t) = I(t) - S$ (Abweichung von Sollstrom $S$), und
* die Stellgröße $Y(t)$, die bestimmt, wie die Stellung der Spitze verändert wird.

"""

# ╔═╡ 81b5f22c-8a6e-11eb-3c82-7d8a54b309cb
random_seed = 12349

# ╔═╡ 81616ba0-8a71-11eb-1f0f-d3bdd3e0fe00
md"""
### Oberfläche generieren.

Hier generieren wir durch Random Walking eine Oberfläche, die wir später abschreiten wollen.
"""

# ╔═╡ b9c9a60c-8a6e-11eb-2977-d7c1e06bd7ef
function generate_step_surface(
		; len::Int64=250, stepwidth::Int64=25, 
		stepheight::Float64=3.)::Vector{Float64}
	rng = Random.MersenneTwister(random_seed)
	current = 0
	thresh = 1.1
	factor = 1.

	function mapstep(x::Float64)
		if x > thresh
			current += stepheight * (x-thresh) * factor
		elseif x < -thresh
			current += stepheight * (x+thresh) * factor
		end
		return current
	end
	
	dst = Random.randn(rng, len)
	for i = 1:length(dst)
		dst[i] = mapstep(dst[i])
	end
	return dst
end

# ╔═╡ b1960dd8-8a6e-11eb-3abc-99286e3e8f1c
plot(generate_step_surface(len=500000))

# ╔═╡ ac96c664-8a6f-11eb-360c-27599ba4d5ce
md"""
### Hauptalgorithmus: Oberflächenverfolgung.

Unabhängig vom Reglermechanismus wird hier eine Oberfläche verfolgt.
"""

# ╔═╡ 03ac3e04-8a73-11eb-3d45-15c63bf21540
# Returns the current as function of distance.

function current_func(dist::Real; current_factor=1)::Float64
	return exp(-current_factor*dist)
end

# ╔═╡ 8e6f403a-8a74-11eb-0045-5318bd996f47
function antilinear_func(dist::Real; current_factor=1)::Real
	return -dist*current_factor
end

# ╔═╡ 94773e40-8a76-11eb-0ff7-250b0d6c5200
md"""
### Regleralgorithmen.

Hier sind die eigentlichen Regelalgorithmen definiert, die die Einstellung der Spitze kontrollieren.
"""

# ╔═╡ ba835830-8a76-11eb-14c9-e58c34a1227f
md"""
`proportional` gibt einen Wert zurück, der mit größerem Fehlersignal `E` stärker abfällt:
"""

# ╔═╡ 76f0ff18-8a72-11eb-3393-8585070a479f
function proportional(E::Real, state=nothing; g=0., h=nothing, j=nothing)
	return g*E, nothing
end

# ╔═╡ 4cb63808-8a72-11eb-1766-23047364c971
# surface: the trackable surface
# pos0: starting position
# target_S: Target current
# current: Function calculating the current from distance
# algo, algoargs: algorithm for tracking

function track_surface(surface::Vector{Float64}, target_S=.4;
		pos0=1, g::Real=.5, h::Real=.1, j::Real=.1,
		current=current_func,
		algo=proportional, algoargs=[])::Vector{Float64}
	pos = pos0
	state = nothing
	
	function maptrack(x)
		# Almost perfect tracking! -- use relative difference 
		# E = (current(pos - x) / current(target_S)) -1 
		
		# Imperfect tracking -- use absolute difference, not good for very small currents
		E = current(pos - x) - current(target_S)

		Y, state = algo(E, state; g=g, h=h, j=j)
		pos += Y
	end

	return collect(map(maptrack, surface))
end

# ╔═╡ 905587a2-8bb0-11eb-17e9-590c36dabc97
@bind rtimespan html"Ti: <input type=range>"

# ╔═╡ f263dcac-8a76-11eb-3116-4396a196a240
begin
	timespan = 4*rtimespan
	
	function integrated(E::Real, state::Union{Nothing, Vector{Float64}};
			g=0., h=0., j=nothing, timespan=timespan)
		if state == nothing
			state = ones(timespan) .* E
		end
		integrated = sum(state)/timespan
		state = [state[2:end]; E]
		return proportional(E, nothing; g=g)[1] + h * integrated, state
	end
end

# ╔═╡ 9eedb29a-8a86-11eb-1a1e-496d57b9ffdb
function differentiated(E::Real,
		state::Union{Nothing, Tuple{Float64, Vector{Float64}}};
		g=0., h=0., j=0., timespan=rtimespan)
	if state == nothing
		last, hist = 0, nothing
	else
		(last, hist) = state
	end
	int, hist_ = integrated(E, hist; g=g, h=h, timespan=timespan)
	dEdx = (E - last)
	
	return int+j*dEdx, (hist_[end], hist_)
end

# ╔═╡ afc29c12-8a76-11eb-3b46-b504dff087de
md"""
### Visualisierung
"""

# ╔═╡ f20ae626-8baf-11eb-3ed1-9901416098b2
@bind rdist html"Abstand: <input type=range>"

# ╔═╡ 02355084-8bb0-11eb-3305-8746fa5e5fe5
@bind rg html"g: <input type=range>"

# ╔═╡ 1d6b897e-8bb0-11eb-06fb-8dd988c72eb4
@bind rh html"h: <input type=range>"

# ╔═╡ ea82dbfc-8bb2-11eb-2252-61479b144fd1
@bind rj html"j: <input type=range>"

# ╔═╡ ade0d464-8a73-11eb-3677-7b5eaf669d81
begin
	
	# Distance between tip and surface
	dist = rdist

	# Weighting of correction factor (pos' = pos - g)
	g = round(rg/3, digits=1)
	
	# Weighting of current by distance (exp(-cf * dist))
	cf = .1
	
	# Which algorithm to use.
	algo = differentiated
	h = rh/5
	j = rj/3

	surf = generate_step_surface(len=1000)
	track = track_surface(surf, rdist; pos0=dist, g=g, h=h, j=j,
		current=(x) -> current_func(x; current_factor=cf),
		algo=algo)
	xs = range(1, length(surf), step=1)
	
	plot(xs, surf, label="Surface", title="z = $dist, g = $g, h = $h, j = $j")
	plot!(xs, track, label="Tip")
	p = plot!(xs, track-surf, label="Distance")
	p
end

# ╔═╡ fe045072-8eee-11eb-3751-53072b09f657
#savefig(p, "../img/simulation_d$(dist)_g$(g)_h$(h)_j$(j).pdf")

# ╔═╡ 13f27f32-9061-11eb-2fc4-c5965a95ba6e
rtimespan

# ╔═╡ 29534788-8a8c-11eb-1d2d-c51484c76caf
md"""
### Konstanthöhenmodus
"""

# ╔═╡ a1066da2-8af3-11eb-2209-b170613f9d8b
begin
	current_factor = .09
	ch_dist = rdist
	ch_surf = generate_step_surface(len=1500)
	ch_scale = maximum(ch_surf)
	scalefunc = identity
	
	currents = map(x -> ch_scale*scalefunc(current_func(ch_dist-x, current_factor=current_factor)),
		ch_surf)
	ch_xs = range(1, length(ch_surf), step=1)
	
	plot(ch_xs, ch_surf, label="Surface", title="z = $ch_dist")
	chp = plot!(ch_xs, currents, label="I", scale=:identity)
end

# ╔═╡ 93a44ae8-9062-11eb-0e77-8b28863fb6e1
#savefig(chp, "../img/simulation_current_d$(dist).pdf")

# ╔═╡ Cell order:
# ╠═3ca61268-8a5d-11eb-34ab-9b0c7ce9c19a
# ╠═662c9a8a-8a5d-11eb-1ca0-972e75cce9de
# ╠═81b5f22c-8a6e-11eb-3c82-7d8a54b309cb
# ╠═81616ba0-8a71-11eb-1f0f-d3bdd3e0fe00
# ╠═b9c9a60c-8a6e-11eb-2977-d7c1e06bd7ef
# ╠═b1960dd8-8a6e-11eb-3abc-99286e3e8f1c
# ╠═ac96c664-8a6f-11eb-360c-27599ba4d5ce
# ╠═03ac3e04-8a73-11eb-3d45-15c63bf21540
# ╠═8e6f403a-8a74-11eb-0045-5318bd996f47
# ╠═4cb63808-8a72-11eb-1766-23047364c971
# ╠═94773e40-8a76-11eb-0ff7-250b0d6c5200
# ╠═ba835830-8a76-11eb-14c9-e58c34a1227f
# ╠═76f0ff18-8a72-11eb-3393-8585070a479f
# ╠═905587a2-8bb0-11eb-17e9-590c36dabc97
# ╠═f263dcac-8a76-11eb-3116-4396a196a240
# ╠═9eedb29a-8a86-11eb-1a1e-496d57b9ffdb
# ╠═afc29c12-8a76-11eb-3b46-b504dff087de
# ╠═f20ae626-8baf-11eb-3ed1-9901416098b2
# ╠═02355084-8bb0-11eb-3305-8746fa5e5fe5
# ╠═1d6b897e-8bb0-11eb-06fb-8dd988c72eb4
# ╠═ea82dbfc-8bb2-11eb-2252-61479b144fd1
# ╠═ade0d464-8a73-11eb-3677-7b5eaf669d81
# ╠═fe045072-8eee-11eb-3751-53072b09f657
# ╠═13f27f32-9061-11eb-2fc4-c5965a95ba6e
# ╠═29534788-8a8c-11eb-1d2d-c51484c76caf
# ╠═a1066da2-8af3-11eb-2209-b170613f9d8b
# ╠═93a44ae8-9062-11eb-0e77-8b28863fb6e1