Julia Community 🟣: Gage Bonner

Firecrackers and Envelopes

Gage Bonner — Fri, 01 Sep 2023 21:59:21 +0000

Introduction

At the risk of going full recipe blog, I'm going to start with a boring but harmless personal anecdote.

When I was a wee lad, I thought that the secret to fame and fortune was solving Project Euler problems. If you aren't familiar, they describe it as,

Project Euler is a series of challenging mathematical/computer programming problems that will require more than just mathematical insights to solve. Although mathematics will help you arrive at elegant and efficient methods, the use of a computer and programming skills will be required to solve most problems.

It would still be a couple years before I discovered Julia, but armed with my rudimentary understanding of CALC101 and even more rudimentary understanding of Python 2.7, I got to work. I'm neither famous nor rich so the obvious conclusion is that I didn't solve enough problems, but I hacked my way through a few. Eventually it became clear that I simply didn't have the algorithmic chops (or patience) to make much more headway so I transitioned from solving problems to just sort of skimming them and seeing if I could think of anything in five minutes. You might be surprised to learn that this actually worked for one (1) problem.

I stumbled onto Problem 317: Firecracker. There were no gigantic prime numbers or arrays of obtusely defined integers, it was just good old projectile motion.

A firecracker explodes at a height of $m100\,\text{m}$ above level ground. It breaks into a large number of very small fragments, which move in every direction; all of them have the same initial velocity of $m/s20\,\text{m/s}$ .

We assume that the fragments move without air resistance, in a uniform gravitational field with $m/s29.81\,\text{m/s}^2$ .

Find the volume (in $m3\text{m}^3$ ) of the region through which the fragments move before reaching the ground. Give your answer rounded to four decimal places.

This problem will be today's topic. I'll show you the solution I came up with all those years ago, and atone for my past sins by also providing a Julia solution. For complete transparency, reading further will spoil the problem!

The Codeless Solution

First, clearly we need some physics - without knowing how a projectile behaves this problem is simply impossible (e.g. what if it had air drag?) A quick skim of the first chapter of any physics textbook will tell us the height of the firecracker $y$ as a function of its horizontal displacement $x$ is given by

y_0 + x \tan \theta - \frac{g \sec^2 \theta}{2 v^2} x^2

where

my_0 = 100\,\text{m}

is the initial height,

20\,\text{m/s}

is the initial velocity,

9.81\,\text{m/s}^2

is the acceleration due to gravity and

θ\theta

is the launch angle with respect to the horizontal. Let's get plotting; we'll assume that the firecracker is launched to the "right" and only look at the part of the trajectory above

y = 0

. I realize that the title of this section is "codeless solution" and I'm about to show some code, but that's showbiz baby.

using Makie, CairoMakie

const y0 = 100.0
const v = 20.0
const g = 9.81

# keep θ in degrees
y(x; θ) = y0 + x*tan(π*θ/180) - g*sec(π*θ/180)^2*x^2/(2*v^2) 

xs = range(0.0, 120.0, step = 0.1)
ys = map(x -> y(x, θ = 20), xs)

fig = Figure()
ax = Axis(
    fig[1,1], 
    limits = (0, 120, 0, 110),
    title = "A firecracker trajectory", 
    xlabel = "horizontal distance from start (m)", 
    ylabel = "distance above ground (m)"
)
lines!(ax, xs, ys, color = :blue, label = L"\theta = 20^\circ")
axislegend(ax)

fig

We see a parabola, which is reassuring. But what's all this about "volume" of the region the projectile moves through? The question advises us that "It breaks into a large number of very small fragments, which move in every direction...", so let's try plotting many trajectories for various launch angles.

fig = Figure()
ax = Axis(
    fig[1,1], 
    limits = (0, 120, 0, 130),
    title = "Many firecracker trajectories", 
    xlabel = "horizontal distance from start (m)", 
    ylabel = "distance above ground (m)"
)

for θ in range(-90.0, 90.0, step = 3.0)
    xs = range(0.0, 120.0, step = 0.1)
    ys = map(x -> y(x, θ = θ), xs)
    lines!(ax, xs, ys)
end

fig

Now the game is clear: the continuum of trajectories traces out a shape, so what we really need is somehow the "boundary" of all of these individual trajectories. If we knew that, we could determine the area traced out by all the trajectories, which is one step closer to the volume. But how to obtain such a thing?

To answer this, we have to think about what exactly is meant by "boundary" of the trajectories. Well, it's certainly some kind of function $B (x)$ . But, given some $x$ , what is the value of $B (x)$ ? If we examine the figure, the "boundary" of the trajectories at a given $x$ looks like the largest value attained by $y (x)$ for any $θ\theta$ . But that's exactly it:

\max_{\theta \in [-\pi/2, \pi/2]}{y(x; \theta) }

So for a given $x$ , we just need to compute the maximum of $y$ , treating it then as a function of $θ\theta$ . This is accomplished by routine calculus,

\frac{\partial y}{\partial \theta} = x \sec \theta^2 - \frac{g \tan \theta\sec^2 \theta}{v^2} x^2 \implies \tan \theta = \frac{v^2}{g x}

Substituting this into

\theta)

gives our boundary

B (x)

y_0 + \frac{v^2}{2g} - \frac{g}{2 v^2} x^2

So (fascinatingly, in my opinion) the boundary of a set of parabolic trajectories is itself a parabola! The boundary is more usually called the envelope and is obtained for a parametric function by simply extremizing the parameter in question (in our case,

θ\theta

) but I certainly didn't know that back in the day.

A neat observation: the quantity $v22g\tfrac{v^2}{2g}$ is actually the height attained by a projectile shot straight up in the air. Therefore, setting $ymax=v22gy_{\text{max}} = \tfrac{v^2}{2g}$ actually allows us to write

y_0 + y_\text{max}\left(1 - \frac{x^2}{4y_\text{max}^2}\right)

Anyway, let's check out $B (x)$ on our plot.

const ymax = v^2/(2*g)

B(x) = y0 + ymax*(1 - x^2/(4*ymax^2)) 

xs = range(0.0, 120.0, step = 0.1)
ys = B.(xs)
lines!(ax, xs, ys, 
    linewidth = 3.0, 
    color = :red, 
    label = L"\mathrm{Boundary: } B(x)")
axislegend(ax)

fig

Well there it is. So the area swept out by the trajectories is the area between the $x$ axis and the function $B (x)$ . But how to get the volume? This is, to date, the absolute only single time I have used a volume of revolution calculation. You may have realized that the calculations we've done so far only hold for a "slice" of the volume. After all, the firecrackers are being ejected in all directions in space. But then, since each direction looks the same as any other, we can get the total volume from the area under $B (x)$ by revolving $B (x)$ around the $y$ axis using the method of cylinders. In brief, this simply states that the volume of a little chunk of cylinder is its circumference $\pi r$ times its height $h$ times its width $dr\text{d} r$ . Noting that $B(2ymax1+y0/ymax)=0B(2y_\text{max} \sqrt{1 + y_0/y_\text{max}}) = 0$ , we therefore integrate

\text{Volume} = \int_{0}^{2y_\text{max} \sqrt{1 + y_0/y_\text{max}}} 2\pi x B(x) \,\text{d}x = 2 \pi \left(1 + \frac{y_0}{y_\text{max}} \right)^2 y_\text{max}^3 ,

a very slick result indeed. We can then obtain the final answer:

2*π*(1 + y0/ymax)^2 * ymax^3

1.8565328455275737e6

The Julia Solution

To solve the problem in a "code-y" way, we'll essentially perform the above calculation, but numerically. First, we'll get $B (x)$ by solving the optimization problem using Optim.jl. Then, we can integrate the result using QuadGK.jl.

using Optim
using QuadGK

boundary(x) = -1*optimize(θ -> -y(x, θ = θ), 0.0, 90.0).minimum
@time quadgk(x -> 2*π*x*boundary(x), 0, 2*ymax*sqrt(1 + y0/ymax))[1]

  0.314286 seconds (286.15 k allocations: 16.542 MiB, 5.87% gc time, 99.95% compilation time)



1.8565328455275744e6

It's just that easy folks.

Conclusion

Even though this problem doesn't end up being that hard, it's been living in my head more or less rent-free for over ten years so it's nice to finally write it down. My Project Euler career has been over for a long time, but there's still some good stuff in there. Let me know if solving them gets you rich and famous.

Creating Nice, Simple Geo-Plots From Scratch using Makie.jl

Gage Bonner — Fri, 28 Jul 2023 17:11:25 +0000

Introduction

It's quite a bit of work to make a plot that looks good. Phrases like "publication-quality" come up a lot in online discussions and one wonders what at what point a plot actually counts as publication-quality. Different journals, websites and people are all going to have different standards but, in my mind, if I have a plot that

has equations, labels and text that are all consistent with default LaTeX rendering, formats and fonts and ...
can be rendered using vector graphics if needed and ...
just in general isn't a complete eyesore ...

then I'm doing pretty well.

It's surprisingly hard to fulfill these requirements. Depending on the package, achieving items 1., 2. and 3. can be in essentially any order of difficulty. Makie.jl is a very promising plotting ecosystem, and in this short tutorial, I will demonstrate how the stock blocks can be used to make an equirectangular plot of Earth's landmasses at high resolution. More to the point, this plot will fulfill items 1., 2. and 3. with below-average inconvenience!

Packages

We will use CairoMakie as the Makie backend to give us vector graphics capabilities right away. Note that this post has been exported from a notebook, so the figures are PNGs, but it is easy to save a figure fig as, for example, a PDF via save("my_fig.pdf", fig). Item 2. already done!

We will also use GeoMakie and in particular GeoMakie.GeoJSON to make parsing the landmass data very easy.

using Makie, CairoMakie
using GeoMakie, GeoMakie.GeoJSON

Building The Axis

The general Makie workflow hierarchy is that plots (lines, heatmaps, contours ...) are stuck an Axis which is stuck on a Figure. We want a nice big figure with large text and plenty of padding on the sides to make room for labels, so let's start with that. In fact, we'll be making this figure several times, so let's put it in a function for easy regeneration.

function default_fig()
    Figure(
        resolution = (1920, 1080), 
        fontsize = 50,
        figure_padding = (100, 100, 100, 100)
    )
end

Arguably not very impressive yet, so let's add an Axis. All that is required of an Axis is that it is placed somewhere specific on a Figure. We want our axis in the first row and the first column of fig, so we simply pass fig[1, 1] to Axis.

fig = default_fig()
ax = Axis(fig[1, 1])
fig

At this point, we should ask a few things of our axis:

It should have a title
It should have the correct limits
The tick marks should be placed where we want
All numbers and text should have the correct formatting as per issue 1.

Fortunately, there is an easy way to demand LaTeX formatting, via L"my-latex-string". Therefore, to get a good-looking title, we can pass to axis something like title = L"\text{My Title}" with math and so on if needed.

What about limits? The world is truly our oyster here, but let's make a neat little box around the equator and Prime meridian. The region of interest will be from 20°W to 20°E by longitude and 20°S to 20°N by latitude. Hence, we can pass limits = (-20, 20, -20, 20).

And tick marks? Every 10° seems like a good separation. We'll use the default range constructor with xticks and yticks.

fig = default_fig()
ax = Axis(
    fig[1, 1],
    title = L"\text{The World, By Land}",
    limits = (-20, 20, -20, 20),
    xticks = -20:10:20,
    yticks = -20:10:20
)
fig

Now we're getting somewhere. The tick marks still don't have the correct font, and they also don't seem to be correctly referring to longitudes and latitudes. That is, -20 on the x axis should actually read 20°W and so on. The way to deal with this is with the handy xtickformat and ytickformat. This allows us to provide a map between the values of the ticks and the string that we actually want displayed.

A basic pattern for this is

xtickformat = values -> [f(value) for value in values]

where f(value) represents we want to actually do to a particular tick mark whose value is value. What do we want in this case? Well, if xtick < 0, we want it to be mapped to abs(xtick) °W with the L"" formatting and similarly for xtick > 0 and the yticks. This is accomplished with a basic if/elseif chain, with the one potential hiccup being that we need to use the % character to interpolate variables inside of L"".

These functions will map x and y ticks to the appropriate string.

function xtick_to_lon(value::Real)
    if value > 0 
        L"%$(abs(value)) \, \degree \mathrm{E}" 
    elseif value == 0 
        L"0\degree"
    elseif value < 0
        L"%$(abs(value)) \, \degree \mathrm{W}" 
    end
end

function ytick_to_lat(value::Real)
    if value > 0 
        L"%$(abs(value)) \, \degree \mathrm{N}" 
    elseif value == 0 
        L"0\degree"
    elseif value < 0
        L"%$(abs(value)) \, \degree \mathrm{S}" 
    end
end

Now we can get the right formatting.

fig = default_fig()

ax = Axis(
    fig[1, 1],
    title = L"\text{The World, By Land}",
    limits = (-20, 20, -20, 20),
    xticks = -20:10:20,
    yticks = -20:10:20,
    xtickformat = values -> [xtick_to_lon(value) for value in values],
    ytickformat = values -> [ytick_to_lat(value) for value in values]
)

fig

Hmm, the bottom left corner looks a bit cramped. It's always something, isn't it? Fortunately, we have the xticklabelpad and yticklabelpad to push the labels back a bit.

fig = default_fig()

ax = Axis(
    fig[1, 1],
    title = L"\text{The World, By Land}",
    limits = (-20, 20, -20, 20),
    xticks = -20:10:20,
    yticks = -20:10:20,
    xtickformat = values -> [xtick_to_lon(value) for value in values],
    ytickformat = values -> [ytick_to_lat(value) for value in values],
    xticklabelpad = 10.0,
    yticklabelpad = 10.0
)

fig

Add padding to taste, but it's looking pretty good!

Land spotted!

We will use Natural Earth Data since it is in the public domain and it has high resolution. The relevant repository is here. We're looking for the .geojson files; a 50 meter resolution should be enough. Simply download the raw data and pipe the result into GeoJSON.read to have a Makie-ready FeatureCollection.

url = "https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/"
landmass = download(url * "ne_50m_land.geojson") |> GeoJSON.read
# FeatureCollection with 1420 Features

A FeatureCollection is easily plottable using the poly! command. We can choose the color using the many built-in named colorant options. Let's go for a classy gray.

poly!(ax, landmass, color = colorant"gray50")
fig

All of a sudden, Africa appears! Looks like the Gulf of Guinea to me, and we even have enough resolution to see the Cameroon Line.

The Axis we have cooked up forms a great baseline for more complicated plots. You can now easily add lines, polygons, and other plots on top of this for basic visualization.

Conclusion

We have created a basic template for equirectangular geo-plots with landmasses already filled in, and were only somewhat inconvenienced. Issues 1., 2. and 3. have been handled to my satisfaction. Job's done, good luck out there!