# List comprehension and integration of a function

I’m currently attempting to integrate a function in Julia using QuadGK() of the following form. It’s computing the potential of a ring of charge in space.
\phi=\int_0^{2\pi} \frac{1}{\sqrt{(x-cos(t))^2+(y-sin(t))^2+z^2)}}dt
If x=cos(t),y=sin(t) and z=0 then obviously my function blows up. To prevent this, while passing a 3d set of data points to the function I used a conditional list comprehension. It looks like so.

inflateunknot(f,N,p,A,t)=[f(x,y,z,A) for x in -1*N:p:N, y in -1*N:p:N, z in -1*N:p:N if z ≠ 0 && y  ≠ sin.(t) && x  ≠ cos.(t) ]


The function integrates just fine without error, however when I go to graph my data I note a problem. This 2d slice of the 3d output matrix A shows a series of points through the origin where no value has been calculated. I assume I’ve got an error in my list comprehension but I’m not too sure why. Should I be using a different tool to disallow these values in my function? I know in python one can use a lambda keyword for a similar issue. All of this is running in a pluto notebook, I include the complete code below.

### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ 6eeb2f90-43c4-11eb-2fb0-23fb34f03217
md"## Introduction"

# ╔═╡ 0848c260-43c0-11eb-218b-59a047c6233f
md"This notebook implements [Max Lipton's](https://e.math.cornell.edu/people/ml2437/) Python code on knot electric potentials in [Julia](https://julialang.org/). The original code it is based on can be viewed at the following [link](https://github.com/ml2437/knot-potential-surfaces/blob/master/Electrostatic%20Knot%20Theory.ipynb). Max has also written two papers on this topic, while a lot of the technical details are advanced and I don't understand it all, there's a lot of useful information in them. They deal with the lower bound on the critical points in a [charged knot] (https://e.math.cornell.edu/people/ml2437/tunnel_number_bound.pdf) and the relationship between the [critical sets and Morse code](https://e.math.cornell.edu/people/ml2437/Morse_Code.pdf) of a charged knot. Future plans for this notebook include restating the math in a geometric algebra formalism, implementing posits instead of floats, and optimizing the code for parallel computation. In general I interleave the translated code with markdown comments that explain the math for my own understanding."

# ╔═╡ a2d3f430-43c5-11eb-0f42-0f689f8fb15e
#setup pkg environment, why? I dunno - I think I can add all sorts of pkgs meow
# or is this necessary
begin
import Pkg
Pkg.activate(mktempdir())
Pkg.Registry.update()
end

# ╔═╡ 7bd300a0-43cb-11eb-153a-1fe35dc585bb
#think I want to use PlutoUI - this gives a bunch of user interface options like
#sliders and stuff that let you interact with your plots and data.
begin
using PlutoUI
end

# ╔═╡ 004c0740-44a4-11eb-2ec2-15a6bd0685a3
begin
using Plots
plotly()
end

# ╔═╡ 59570b70-4b91-11eb-252a-db35d8b6e956
begin
using GeometryBasics
using LinearAlgebra
using Images
end

# ╔═╡ b5a0612e-4419-11eb-118b-9db02b2086ea
begin
end

# ╔═╡ 2fc65542-4ae7-11eb-1e6e-5f53acc8d0df
# this might require some additional packages, in particular geometrybasics, geometrytypes and staticarrays
begin
using Meshing
end

md"## Knot Parametrizations"

# ╔═╡ 90400640-43c7-11eb-10a4-831d39da8bde
#Order of accuracy for Gaussian Quadrature - don't know why it shows sometimes
quadorder = 20 #like python Julia doesn't require ';'

# ╔═╡ 897f19d0-43d7-11eb-392c-6fee25b1eb92
#Discretization of the domain interval for knot parametrizations
t = range(0, stop = 2*pi, length = 1000)# an array of a thousand components

# ╔═╡ 29152b00-43ca-11eb-0221-cd0ccf8aa735
md"We begin by defining the paramaterizations of a number of different knots. In doing so we use the discretization variable t we defined above which is a linear array of 1000 elements from 0 to 2$\pi$. This means that we've swept out an entire circle in 1000 discrete points. The quadorder variable is used for integration later on. We are going to define the following knot types, the Unknot, the Trefoil, the Figure 8, the (3,1) Torus Knot, the (5,1) Torus Knot, the Cinquefoil Knot, the 3-Twist Knot, the Granny Knot, the Square Knot and finally the Endless Knot. Each knot has 3 functions for returning x,y and z coordinates and a numerator function that returns the  derivative. "

# ╔═╡ dd692ae2-440e-11eb-1c97-dd1f4c7e6178
md"#### The Unknot"

md"The formula for a circle centered at the origin is of course $r=1$ in polar coordinates. Here to output $x$ and $y$ values we use $x = cos(t)$ and $y = sin(t)$. Our value for $z$ is $0$ since the unknot sits on an axis and has no depth."

# ╔═╡ 20492eb0-440e-11eb-06de-c147e3f93660
function unknotnum(t) #the numerator in the integral below
1
end

# ╔═╡ 5fe42430-440e-11eb-228b-53a81b657786
function unknotx(t)
cos.(t) #need the dot for elementwise
end

# ╔═╡ 68cb0aa0-440e-11eb-1565-c94c7e099a0a
function unknoty(t)
sin.(t)
end

# ╔═╡ 6eba5bf0-440e-11eb-03d6-4d2f8e8fc9f1
function unknotz(t)
0
end

# ╔═╡ 2b1c2810-449a-11eb-1ea3-2bed99095cff
md"The Unknot is a simple circle, we can visualise the structure using our functions above. "

# ╔═╡ bdd1fb60-493d-11eb-302d-aba1d7e3d74f
plot(unknotx(t),unknoty(t))
# want to improve the presentation, ideally it'd be in 3d and scaled symetrically

# ╔═╡ ba86cf00-440e-11eb-033c-8fd03e12333b
md"Having established a parameterization we can now move to integration to compute the electric potential at a specific point (a,b,c). The core task here is how to integrate the electric potential. Recall that the formula for computing an electric potential $\phi_{E}= -\int_{C} E \cdot dl$, which is a line integral. We can ignore the negative for convenience purposes. The electric potential is the voltage at a point in space. As in an electric circuit where we have to measure voltage relative to a reference point (there called ground), here we measure voltage relative to a point at $\infty$ where $\phi=0$. Of course to calculate the electric potential or voltage at a point, we have to know what the electric field is at that point. Having parametrized our knot shapes, we can calculate the electric field by measuring the distance of a point in space from the knot, and then assigning a fixed charge to the knot itself. The Unknot is simply a ring of charge, a special case of which I solved during my Electrical Engineering degree. Now we need to solve the general case in 3 dimensions for the electric field all round the ring of charge rather than at a specific point. "

md"We will simply state that a point $x$ at some distance $R$ from an electric charge $k$ has an electric potential proportional to $R^{-1}$. This gives us the integral below, where $r(t)$ is the function parametrizing our knot and $x$ is some point. This is the general formula for any knot parametrization, we will need to implement a specific version for the Unknot. "

# ╔═╡ 43aac0b0-4492-11eb-2b05-c51da990da76
md"$\phi=\int_0^{2\pi} \frac{|r'(t)|}{|x-r(t)|}dt$"

# ╔═╡ 6f51d500-4492-11eb-01a6-d9925465c853
md" In our case we're working in $\mathbb{R}^3$ so that means we need to implement a distance including $x,y$ and $z$ coordinates. Our formula for the Unknot will thus end up as below. With the implementation of the function to be integrated written in Julia immediately below it."

md"$\phi=\int_0^{2\pi} \frac{1}{\sqrt{(x-cos(t))^2+(y-sin(t))^2+z^2)}}dt$"

# ╔═╡ a33a0e30-4895-11eb-18c8-d3cf673935d9
function unknotintegralfunction(a,b,c,t)
1/(sqrt(((a.-cos.(t))^2).+((b.-sin.(t))^2).+c^2)) #using .math means element wise
end

md"We've skipped some details in the derivation, but via cursory examination it should be apparent that we're calculating a distance in $\mathbb{R}^3$. ie: see that we're calculating our distance in the $x$ coordinate by subtracting $cos(t)$, the position of the knot, from our position in $x$. So the denominator of our fraction is the tried and true pythagorean formula. Since the Unknot is centered at the origin and has no height in $z$ we don't subtract anything from it since one's location in $z$ is always the distance from $0$."

# ╔═╡ 587dcdd0-440f-11eb-2db6-19abd03ce8c2
md"To actually integrate the formula above we will need an integration function for a single variable in Julia. The best current path seems to be a quadrature routine using the [QuadGK package](https://juliahub.com/ui/Packages/QuadGK/hq5Ol/2.4.1). "

# ╔═╡ 68d75000-4619-11eb-33b2-094cac9629bc
md"There's all sorts of fancy terminology and I need to understand how this works. To restate the basics, I'm integrating a value at a point in 3 dimensions $x$, $y$, and $z$ over a range from $[0,2\pi)$. The single variable is the $t$ of my parametrization function. So in my function below I pass the position (a,b,c) to the function which provides the point in space $t$ over the range specified by the second and third variables. Note the $t \rightarrow$, this tells Julia that the range for quadgk $0$ to $2\pi$ is assigned to the variable $t$. "

# ╔═╡ ded6aae0-4be0-11eb-266c-93a8daeaf540
#setup our storage matrix
A = zeros(41,41,41)

# ╔═╡ 50965470-488d-11eb-3b0a-313cbaf71261
function unknotpotential(a,b,c,A)
# we'll use the default errors to start
A[floor(Int,(10*a+21)),floor(Int,(10*b+21)),floor(Int,(10*c+21))]=ans# need a way to store a,b,c with the ans
return [a,b,c,ans,err]
end

# ╔═╡ 02920c30-4be2-11eb-244f-4de20ca35dd6
md"We have to generate a 3 dimensional array to store all of the data for our marching cubes algorithm later on. This array will be passed to the functions in question and they'll write the potential value to the specific location in the array. "

# ╔═╡ ecf19c60-4984-11eb-2919-771158df41e2
md"The integration function currently returns the correct value for the potential at the origin. The origin for the unknot is a critical point, so its natural for us to ask why did the integral return the value of $2\pi$ rather than $0$? Recall the definition of a [critical point](https://en.wikipedia.org/wiki/Critical_point_(mathematics)): for a real valued function a critical point is a point where the function is either not differentiable or its derivative is equal to $0$. Above I've calculated the critical value, which is the value of the function at the critical point. So there is a difference between the electric potential, $\phi$ and the electric field, $E = -\nabla\phi$. The electric field is the gradient of the electric potential, so the critical points of the potential are where the **electric field** is $0$, not where the electric potential is $0$. Since the electric potential is a voltage, which by definition is a relative value, the function is saying that there are $2\pi$ volts at the origin relative to a location at $\infty$. "

md"We can now try and calculate the electric potential for a whole region of space. Unlike older languages such as Matlab or Python you don't need to use a meshgrid type function in Julia. We can instead create our grid in $\mathbb{R}^3$ using list comprehension. This is a fun little programming technique that aims to copy set builder notation. So here I pass a function and a dimension to my inflate function, which then goes and gives all of these values to my passed function. A little bit of Julia syntax to take note of, for the range for $x$,$y$ and $z$ below we have $-1*N:0.1:N$. Which tells us to generate a range from $-1*N$ to $N$ in steps of $0.1$. So if I choose $N=4$ we would have a range $-4:0.1:4$ which would actually be $-4,-3.9,-3.8,...3.8,3.9,4$. We have to be careful with the range we feed this function, since we have a fraction of the form  $\frac{1}{a}$ where $a$ is our distance from the ring, for very small values of $a$ or when $a=0$ we're going to have a problem of blowup. Therefore we need to disallow some values from being included in the set. To do that we can apply some conditions to our list comprehension. In set builder notation we'll write: $\{x,y,z | x,y,z \in \mathbb{R} \wedge x \neq cos(t) \wedge y \neq sin(t) \wedge z \neq 0 \}$. Which translates quite nicely and directly into the Julia code one sees below. "

# ╔═╡ 7030e220-4a6b-11eb-1db0-6399f284f4bc
md"Having successfully calculated the electric potential at all of these points we now have an enormous tuple. The tuple is itself made up of tuples. The first three elements of each element in the tuple are the position in space and the last two are the electric potential and the error in the integration calculation.  "

# ╔═╡ 379a5fb0-4be5-11eb-1338-f360ec868a99
A #Just confirming data was written to the array

# ╔═╡ 5e06cf02-4a6c-11eb-24dd-6b1bceb4f832
md"The next step is to visualize the data. Max calculated a level potential surface. Recall the definition of a [level surface](https://web.ma.utexas.edu/users/m408m/Display12-6-4.shtml): For a function $\phi = f(x,y,z): U \subseteq \mathbb{R}^3 \rightarrow \mathbb{R}$ the level surface of value $c$ is the surface $s$ in $U \subseteq \mathbb{R}^3$ on which $\phi \rvert_{s}=c$. So in our case we're collecting all the points that have the same voltage."

# ╔═╡ 89db1b90-4a71-11eb-0280-e9986755c8ef
md"The trick of course is to generate the surface to plot. How do we collect those points? To do that Max used the [marching cubes algorithm](https://en.wikipedia.org/wiki/Marching_cubes). There is of course a package that gives us this capability in Julia called [meshing](https://juliageometry.github.io/Meshing.jl/stable/). There's a lovely introduction to marching cubes and the isosurface problem at the following [link](https://0fps.net/2012/07/12/smooth-voxel-terrain-part-2/), and if you'd like a very good exposition of the technical details involved you can investigate the classic explanation by [Paul Bourke](http://paulbourke.net/geometry/polygonise/). Very basically the marching cubes algorithm takes some sort of input, either an equation or an array, and outputs a mesh. The technical name for this problem is isosurface extraction."

md"The first bit of important information for us is to note that the scalar values that we're going to use to generate our manifold are stored at each vertex. Recall from graph theory that a vertex is a node, and an edge is just an edge. So for a generic cube, we'd have 8 vertices (corners), and 12 edges connecting all the vertices. What the marching cubes algorithm does is generate a surface based on a boundary the user specifies. So I should be able to specify a value for the electric potential, and find a boundary around the ring where that value holds. To do that we need to grab the matrix A we defined earlier and feed it into the marching cubes algorithm. "

# ╔═╡ a4cac4a0-4be8-11eb-04e6-a7328783b6aa
md"We'll just copy what Max did for the level value of our isosurface, which is to take the potential value at the origin and arbitrarily add 0.5 to it."

# ╔═╡ 5b0caa10-4be6-11eb-2047-6deab2f3b7fd
c=A[20,20,20]+0.5

# ╔═╡ ba45dac0-4b27-11eb-1675-ef9d2ae9f9ca
#generate a mesh
points,surfaces = isosurface(A,MarchingCubes(iso=c))

# ╔═╡ d4f6ebe0-4be8-11eb-1239-3320defd23ba
md"Now that the points and surfaces of the mesh have been generated we need to plot them! To do that we can use the Mesh3d function. And to access the values in my points tuple I'm going to have to use a list comprehension. Having provided the points mesh3d will autogenerate the mesh for me."

# ╔═╡ f215b4d0-4bfd-11eb-0a33-d5711fbef747
mesh3d([a for (a,_,_) in points], [b for (_,b,_) in points], [c for (_,_,c) in points])

# ╔═╡ b2bc6ea0-4d2e-11eb-32da-bb24b71c18c4
md"This show method along with the BWImage function were copied from the following [spot](https://gist.github.com/pbouffard/3d48d3c47d9bd70e7c9f52f984d14245). It apparently comes from the [2020 JuliaCon talk](https://www.youtube.com/watch?v=IAF8DjrQSSk). Either way now that it's included I've got a very quick and dirty way to view a 2D matrix as an image."

# ╔═╡ 522123b0-4d2e-11eb-3fc4-c946ba08cbaa
begin
struct BWImage
data::Array{UInt8, 2}
zoom::Int
end
function BWImage(data::Array{T, 2}; zoom::Int=1) where T <: Real
BWImage(floor.(UInt8, clamp.(((data .- minimum(data)) / (maximum(data) .- minimum(data))) * 255, 0, 255)), zoom)
end

import Base: show

function show(io::IO, ::MIME"image/bmp", i::BWImage)

orig_height, orig_width = size(i.data)
height, width = (orig_height, orig_width) .* i.zoom
datawidth = Integer(ceil(width / 4)) * 4

palette_size = 256 * 4
data_size = datawidth * height * 1

write(io, 0x42, 0x4d)
write(io, 0x00, 0x00)
write(io, 0x00, 0x00)

write(io, Int32(width))
write(io, Int32(-height))
write(io, UInt16(1))
write(io, UInt16(8))
write(io, UInt32(0))
write(io, UInt32(0))
write(io, 0x12, 0x0b, 0x00, 0x00)
write(io, 0x12, 0x0b, 0x00, 0x00)
write(io, UInt32(0))
write(io, UInt32(0))

# color palette
write(io, [[x, x, x, 0x00] for x in UInt8.(0:255)]...)

# data
padding = fill(0x00, datawidth - width)
for y in 1:orig_height
for z in 1:i.zoom
line = vcat(fill.(i.data[y,:], (i.zoom,))...)
end
end
end
end

# A black and white image of the potential
BWImage(A[:,20,:], zoom=10)

# ╔═╡ 6b5d40c0-4d2e-11eb-3ffc-3b2857b90f7b
#A different view of the potential
BWImage(A[:,:,20], zoom=10)

# ╔═╡ 2fa8fa52-4d39-11eb-1014-4f20761386d3
BWImage(A[18,:,:], zoom=10)

# ╔═╡ e158b030-48d2-11eb-3a88-5f96538a8c9c
md"Following the calculation of the actual electric potential, I want to analyze the orbits of the vector field. To do that I need to restate the equation for the electric potential as 3 separate equations for $\frac{dx}{dt}$, $\frac{dy}{dt}$ and $\frac{dz}{dt}$. Recall that $t \in (0,2\pi]$, which is sort of unnecessary to state since $sin(t)$ and $cos(t)$ are periodic, but forewarned is forearmed."

# ╔═╡ 0345b880-4934-11eb-326d-79484684c8e9
md"$\frac{dx}{dt}=\frac{1}{x-cos(t)}$"

md"$\frac{dy}{dt}=\frac{1}{y-sin(t)}$"

md"$\frac{dz}{dt}=\frac{1}{z}$"

# ╔═╡ 4459aab0-493a-11eb-1982-1965fec7b47e
md"Now that we've got these 3 equations setup we can start searching for some zeros! Since Max has already done the work we know that the only critical point for this potential is at the origin. Is it possible to show that algebraically using the techniques I learned in my ODE's course?"

# ╔═╡ 86323420-493a-11eb-1bec-b5f36fdf424d
md"#### The Trefoil"

# ╔═╡ fe126c90-4b96-11eb-315a-f31f194d9e1b
points

# ╔═╡ 20df66d0-4cb7-11eb-2b0f-cdd79ddbfe53
x = A[20,:,:]

# ╔═╡ 4a56c450-4a0e-11eb-2432-fde548b9f590
#we're going to plagiarize this function from Tomas' answer in google groups
#N should be a positive integer
inflateunknot(f,N,p,A,t)=[f(x,y,z,A) for x in -1*N:p:N, y in -1*N:p:N, z in -1*N:p:N if z ≠ 0 && y  ≠ sin.(t) && x  ≠ cos.(t) ]

# ╔═╡ e7a1b8be-4a1b-11eb-16cc-0b68c71c2418
#now we're going to call the function itself
ans = inflateunknot(unknotpotential,2,0.1,A,t)

# ╔═╡ Cell order:
# ╟─6eeb2f90-43c4-11eb-2fb0-23fb34f03217
# ╟─0848c260-43c0-11eb-218b-59a047c6233f
# ╠═a2d3f430-43c5-11eb-0f42-0f689f8fb15e
# ╠═7bd300a0-43cb-11eb-153a-1fe35dc585bb
# ╠═004c0740-44a4-11eb-2ec2-15a6bd0685a3
# ╠═59570b70-4b91-11eb-252a-db35d8b6e956
# ╠═90400640-43c7-11eb-10a4-831d39da8bde
# ╠═897f19d0-43d7-11eb-392c-6fee25b1eb92
# ╟─29152b00-43ca-11eb-0221-cd0ccf8aa735
# ╟─dd692ae2-440e-11eb-1c97-dd1f4c7e6178
# ╠═20492eb0-440e-11eb-06de-c147e3f93660
# ╠═5fe42430-440e-11eb-228b-53a81b657786
# ╠═68cb0aa0-440e-11eb-1565-c94c7e099a0a
# ╠═6eba5bf0-440e-11eb-03d6-4d2f8e8fc9f1
# ╟─2b1c2810-449a-11eb-1ea3-2bed99095cff
# ╠═bdd1fb60-493d-11eb-302d-aba1d7e3d74f
# ╟─ba86cf00-440e-11eb-033c-8fd03e12333b
# ╟─43aac0b0-4492-11eb-2b05-c51da990da76
# ╟─6f51d500-4492-11eb-01a6-d9925465c853
# ╠═a33a0e30-4895-11eb-18c8-d3cf673935d9
# ╟─587dcdd0-440f-11eb-2db6-19abd03ce8c2
# ╠═b5a0612e-4419-11eb-118b-9db02b2086ea
# ╟─68d75000-4619-11eb-33b2-094cac9629bc
# ╠═ded6aae0-4be0-11eb-266c-93a8daeaf540
# ╠═50965470-488d-11eb-3b0a-313cbaf71261
# ╟─02920c30-4be2-11eb-244f-4de20ca35dd6
# ╟─ecf19c60-4984-11eb-2919-771158df41e2
# ╠═4a56c450-4a0e-11eb-2432-fde548b9f590
# ╠═e7a1b8be-4a1b-11eb-16cc-0b68c71c2418
# ╟─7030e220-4a6b-11eb-1db0-6399f284f4bc
# ╠═379a5fb0-4be5-11eb-1338-f360ec868a99
# ╟─5e06cf02-4a6c-11eb-24dd-6b1bceb4f832
# ╟─89db1b90-4a71-11eb-0280-e9986755c8ef
# ╠═2fc65542-4ae7-11eb-1e6e-5f53acc8d0df
# ╟─a4cac4a0-4be8-11eb-04e6-a7328783b6aa
# ╠═5b0caa10-4be6-11eb-2047-6deab2f3b7fd
# ╠═ba45dac0-4b27-11eb-1675-ef9d2ae9f9ca
# ╠═d4f6ebe0-4be8-11eb-1239-3320defd23ba
# ╠═f215b4d0-4bfd-11eb-0a33-d5711fbef747
# ╟─b2bc6ea0-4d2e-11eb-32da-bb24b71c18c4
# ╟─522123b0-4d2e-11eb-3fc4-c946ba08cbaa
# ╠═6b5d40c0-4d2e-11eb-3ffc-3b2857b90f7b
# ╠═2fa8fa52-4d39-11eb-1014-4f20761386d3
# ╟─e158b030-48d2-11eb-3a88-5f96538a8c9c
# ╟─0345b880-4934-11eb-326d-79484684c8e9
# ╠═4459aab0-493a-11eb-1982-1965fec7b47e
# ╟─86323420-493a-11eb-1bec-b5f36fdf424d
# ╠═fe126c90-4b96-11eb-315a-f31f194d9e1b
# ╠═20df66d0-4cb7-11eb-2b0f-cdd79ddbfe53



Should those &&'s be ||'s?

I don’t think so? I’m trying to prevent only the ring of charge being included. I think in that case a coordinate would have the form x,y,z = cos(t), sin(t), 0 I think I’d need &&'s for that to work no?

You should change the coordinate system to have radial coordinate r which will have singularity at ring radius r = 1, to be able to approach it more easily. Otherwise, you can just change the if in the comprehension to:

if z ≠ 0 || (x*x + y*y) ≠ 1.0


PS: reminding that the negation of “A and B” is “Not A or Not B.”. I.e., !(z == 0 && (x*x +y*y) == 1) === z ≠ 0 || (x*x +y*y) ≠ 1

1 Like

Thank you Rafael, the change in the if statement removed the problem.

I’ve been rereading the documentation on list comprehensions and I think I was getting mixed up with my use of \neq, apologies for the confusion!