[question] Finding number of closed potential contours and area enclosed by them

Specific Domains Visualization
, , , ,
1

Hello, I’m trying to calculate the number of closed contours in a potential as well the area enclosed by closed,

i.e. In this image below at certain energy level there are 4 closed contours. I’m trying to figure a way to calculate the number of closed contours (4 here) and the area enclosed by them, given the energy levels.

con_1
con_1680×400 63.6 KB

sample code for above (adopted from Pluto.jl notebook example),

using  GLMakie

len = 50
x = range(-1, 1, length=len)
y = range(-1, 1, length=len)
four_wells(x,y) = x^4 -  x^2 + y^4 - y^2
z = [four_wells(x[i], y[j]) for i in 1:len, j in 1:len]

# energy
height = -0.4

f = Figure(;resolution=(680,400))
ax1 = Axis3(f[1, 1])
ax2 = Axis(f[1, 2])

# Draw the contour in the 3D graph
plane = [height for i in 1:length(x), j in 1:length(y)];

# Drawing the function and the contour on the 3G graph
surface!(ax1, x, y, z, colormap=:viridis)
surface!(ax1, x, y, plane; transparency=true, colormap=[:white,:white])

# Drawing the function and the countour on the heatmap
hm = contourf!(ax2, x, y, z)
Colorbar(f[:, 3], hm; vertical=true)
contour!(ax2, x, y, z; linewidth=3, levels=[height], color=:white)
f
Pluto notebook file for the same 
### A Pluto.jl notebook ###
# v0.19.35

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 iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ ca18be6d-8a9b-4c53-bd2e-553ea01ebb13
using PlutoUI, GLMakie

# ╔═╡ 503aa0a0-ece8-4acc-b437-3119055d22bc
k = @bind height PlutoUI.Slider(-1:0.1:1; default=0, show_value=true)

# ╔═╡ 2dab2fcd-77bf-4651-8223-6854863f425a
md"""
# Appendix
"""

# ╔═╡ 58d43864-9a74-42da-96d5-afb59f22cd0f
begin
	len = 50
	x = range(-1, 1, length=len)
	y = range(-1, 1, length=len)
	four_wells(x,y) = x^4 -  x^2 + y^4 - y^2
	z = [four_wells(x[i], y[j]) for i in 1:len, j in 1:len]
end

# ╔═╡ e39373bd-8e12-469f-bd66-2d005ff6cff2
begin
	# Setup of the figure
    f = Figure(;resolution=(680,400))
    ax1 = Axis3(f[1, 1])
    ax2 = Axis(f[1, 2])

	# Draw the contour in the 3D graph
	plane = [height for i in 1:length(x), j in 1:length(y)];
	
	# Drawing the function and the contour on the 3G graph
	surface!(ax1, x, y, z, colormap=:viridis)
	surface!(ax1, x, y, plane; transparency=true, colormap=[:white,:white])

	# Drawing the function and the countour on the heatmap
	hm = contourf!(ax2, x, y, z)
	Colorbar(f[:, 3], hm; vertical=true)
    contour!(ax2, x, y, z; linewidth=3, levels=[height], color=:white)
end;

# ╔═╡ 9c17c0d8-e4cb-41c0-b180-b18a3c87e265
# Displaying the function
f

# ╔═╡ Cell order:
# ╠═503aa0a0-ece8-4acc-b437-3119055d22bc
# ╟─9c17c0d8-e4cb-41c0-b180-b18a3c87e265
# ╠═e39373bd-8e12-469f-bd66-2d005ff6cff2
# ╟─2dab2fcd-77bf-4651-8223-6854863f425a
# ╠═ca18be6d-8a9b-4c53-bd2e-553ea01ebb13
# ╠═58d43864-9a74-42da-96d5-afb59f22cd0f

There’s a python implementation on stackoverflow, any pointers on which packages to look on Julia will be more helpful.

2

I’m not familiar with the internals of Makie’s contour code, so I can’t offer a solution as elegant as the python equivalent. Nevertheless, here is an idea for how to do a the calculation using other Julia tools:

  • To obtain a representation of the contour, MDBM.jl will be able to sample the iso-surface of the potential
  • The list of points given by MDBM.jl will not distinguish the number of closed contours, so you will have to cluster the points into disconnected contours. It should be fairly easy to write a function for this if there aren’t any packages available.
  • With the list of points in each contour, you can sort these points in a counterclockwise fashion and calculate their area, e.g. using the shoelace formula as in the python example

I hope these ideas help

1 Like
3

Funny, this seems the fill pits steps in DEMs. Found this package WhereTheWaterFlows that has to solve it as an intermediate step.

1 Like