Update blog post

_posts/guides/interactive-blog/2025-02-13-interactive-julia-plotting.md

@@ -25,10 +25,11 @@ which provides a step-by-step guide on how to create interactive plots on a web
 backend for the [Makie.jl](https://docs.makie.org/stable/) plotting library. 
 
 Aaron's blog post is now 4 years old and the method described there unfortunately no 
-longer works. I was determined, but since I barely know what HTML is I needed a little help. 
-I found out yet again how great, supportive and helpful the Julia community is. After 
-posting on the [Julia Discourse](https://discourse.julialang.org/)
-I got [a response from Simon Danisch](https://discourse.julialang.org/t/exporting-figures-to-static-html/125896/16?u=langestefan) who is the creator of Makie.jl. 
+longer works. I was determined, but I needed a little help, and I found out again how 
+supportive the Julia community is. 
+[Simon Danisch](https://discourse.julialang.org/t/exporting-figures-to-static-html/125896/16?u=langestefan) 
+who is the creator of Makie.jl, provided guidance on how to set things up and this blog 
+post is a result of that.
 
 {% alert note %}
 Besides WGLMakie.jl we will also need <a href="https://github.com/SimonDanisch/Bonito.jl">Bonito.jl</a>
@@ -37,10 +38,6 @@ to create the HTML descriptions, which will enable us to embed the plot in a blo
 
 ## A first example
 
-This blog post, which can be read as a tutorial or recipe, is a direct translation from 
-Simon's response into a working example. I hope it helps you to create interactive plots
-for your blog posts.
-
 First, we need a location to store the script that generates our plots. I prefer to 
 group all files for a specific post in a single folder. For this post, I created a 
 folder called `_posts/guides/interactive-blog/` and saved the script as `plots.jl`. As
@@ -132,7 +129,8 @@ $\vec{x}(t)$.
 There are a bunch of interesting questions we can ask in this general problem context.
 A few that come to mind are:
 - How does the trajectory of the ball change when we change the initial velocity $\vec{v}_0$?
-- Given some initial velocity $\vec{v}_0$, what is the maximum horizontal distance the ball can travel?
+- Given some initial velocity $\vec{v}_0$, what is the angle of the kick that maximizes 
+  the horizontal distance the ball will travel?
 - Given (noisy) observations of a ball's trajectory, can we estimate what the initial 
   position $\vec{p}_0$ and velocity $\vec{v}_0$ were? Or even more interesting, can we
   predict where the ball will land while it is still in the air?
@@ -172,7 +170,7 @@ $$
 \end{equation}
 $$
 
-Where $C_L$ is the lift coefficient, $\theta$ is the angle of attack and $f(\theta)$ is a
+Where $C_L$ is the lift coefficient and $f(\theta)$ is a
 function that depends on the spin angle of the ball. Let's assume that the ball is spinning
 around the $z$-axis, then $f(\theta) = [-1, 1, 0]^T$ where we assume that the dependence
 on the angular velocity $\omega$ is already included in $C_L$.
@@ -193,15 +191,6 @@ constant $H = \frac{\rho A}{2m}$, which will simplify the equations.
 Using the fact that $\frac{dx}{dt} = v_x, \frac{dy}{dt} = v_y, \frac{dz}{dt} = v_z$ the 
 system of differential equations can then be written as:
 
-<!-- $$
-\begin{align}
-    H &= \frac{\rho A}{2m} \\
-    \left| \vec{v} \right| &= \sqrt{v_x^2 + v_y^2 + v_z^2}
-    % C_L &= \left| \vec{v} \right|^{-1} \omega r \\ % C_L &= \frac{\omega r}{\left| \vec{v} \right|}\\
-    % \omega &= \omega_0 \cdot \exp\left(-\frac{t}{\tau}\right)    
-\end{align}
-$$ -->
-
 $$
 \begin{align}
     \vec{a} &= -\left | \vec{v}\right | H      
@@ -214,20 +203,34 @@ $$
 $$
 
 The value of $C_L$ is usually derived from experiments and depends on the speed and 
-angular velocity of the ball. We will assume that $C_L$ takes the following form:
+angular velocity of the ball. We will assume that $C_L$ takes the following simplified 
+form:
 
 $$
 \begin{align}
     C_L &= \frac{\omega r}{\left| \vec{v} \right|}\\
-    \omega &= \omega_0 \cdot \exp\left(-\frac{t}{\tau}\right)
+    \omega &= \omega_0 \cdot \exp\left(-\frac{t}{7}\right)
 \end{align}
 $$
 
-Where $\omega_0$ is the initial angular velocity and $\tau$ is the time constant. We 
-will solve this system numerically to get the trajectory of the ball. 
+Where $\omega_0$ is the initial angular velocity just after the kick. We will solve this 
+system numerically to obtain the trajectory of the ball after kicking it.
 </div>
 
-Now we can finally start writing code! Which should be a breeze after because we have 
-access to [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/).
+For those playing along at home, a table with all problem constants is given below:
+
+| Symbol | Description | Value |
+|--------|-------------|-------|
+| $m$ | mass of the ball | 0.43 kg |
+| $g$ | acceleration due to gravity | 9.81 m/s$^2$ |
+| $\rho$ | air density | 1.225 kg/m$^3$ |
+| $A$ | cross-sectional area of the ball | 0.013 m$^2$ |
+| $C_D$ | drag coefficient | 0.2 |
+| $r$ | radius of the ball | 0.11 m |
+| $\omega_0$ | initial angular velocity | 88 rad/s |
+
+
+Now we can finally start writing code! Which should be a breeze because we have 
+access to [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/),
+a package that provides a high-level interface for solving differential equations. 
 
-{% include_code file="_posts/guides/interactive-blog/plots.jl" lang="julia" start="30" end="50" %}

_posts/guides/interactive-blog/plots.jl

@@ -35,184 +35,104 @@ end
 
 # plot 3 - differential equation
 
-# using DifferentialEquations
-
-# # struct to hold problem constants
-# struct SoccerConst
-#     m::Float64 # mass of the ball
-#     g::Float64 # acceleration due to gravity
-#     r::Float64 # radius of the ball
-#     ρ::Float64 # density of the air
-#     C_L::Float64 # lift coefficient
-#     C_D::Float64 # drag coefficient
-# end
-
-# # inital conditions
-# struct SoccerIC
-#     x::Float64 # initial x position
-#     y::Float64 # initial y position
-#     z::Float64 # initial z position
-#     v_x::Float64 # initial x velocity
-#     v_y::Float64 # initial y velocity
-#     v_z::Float64 # initial z velocity
-# end
-
-# open("_posts/examples/diffeqviz/contour.html", "w") do io
-#     app = App() do 
-#         markersize = Bonito.Slider(range(10, stop=100, length=100))
-
-#         # Create a scatter plot
-#         fig, ax = meshscatter(rand(3, 100), markersize=markersize)
-
-#         # Return the plot and the slider
-#         return Bonito.record_states(session, DOM.div(fig, markersize))
-#     end;
-
-#     as_html(io, session, app)
-# end
-
-# # plot 2 - volume plot
-# open("_posts/examples/diffeqviz/sinc_surface.html", "w") do io
-#     # create sub session
-#     sub = Session(session)
-
-#     app = App(volume(rand(10, 10, 10)))
-#     as_html(io, sub, app)
-# end
-
-# open("_posts/examples/diffeqviz/diffeq.html", "w") do io
-#     # Create a Bonito session
-#     sub = Session(session)
-
-#     # Define the Lorenz system
-#     function lorenz!(du, u, p, t)
-#         x, y, z = u
-#         σ, ρ, β = p
-#         du[1] = σ * (y - x)
-#         du[2] = x * (ρ - z) - y
-#         du[3] = x * y - β * z
-#     end
-
-#     # Time span and initial condition
-#     tspan = (0.0, 40.0)
-#     u0    = [1.0, 0.0, 0.0]
-
-#     # Default parameters (classic Lorenz with adjustable β)
-#     σ₀ = 10.0
-#     ρ₀ = 28.0
-#     β₀ = 2.666
-#     p₀ = (σ₀, ρ₀, β₀)
-
-#     # Create the initial ODE solution
-#     prob = ODEProblem(lorenz!, u0, tspan, p₀)
-#     sol  = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, saveat=0.05)
-
-#     # Create a Makie figure with two rows:
-#     # Row 1: 3D scene (LScene) with the Lorenz attractor.
-#     # Row 2: A slider grid controlling the parameters.
-#     fig = Figure(size = (900, 600))
-
-#     # --- Row 1: 3D Plot ---
-#     ax = LScene(fig[1, 1], show_axis=true)
-#     lineplot = lines!(ax,
-#         sol[1, :], sol[2, :], sol[3, :],
-#         linewidth = 2
-#     )
-
-#     # --- Row 2: Slider Grid ---
-#     sgrid = SliderGrid(fig[2, 1],
-#         (label = "σ", range = LinRange(0, 20, 100)),
-#         (label = "ρ", range = LinRange(0, 50, 100)),
-#         (label = "β", range = LinRange(0, 10, 100))
-#     )
-#     σ_slider, ρ_slider, β_slider = sgrid.sliders
-
-#     # Set the slider default values.
-#     σ_slider.value[] = σ₀
-#     ρ_slider.value[] = ρ₀
-#     β_slider.value[] = β₀
-
-#     # Function to re-solve the ODE and update the plot.
-#     function update_plot!()
-#         # Get the current parameter values from the sliders
-#         σ = σ_slider.value[]
-#         ρ = ρ_slider.value[]
-#         β = β_slider.value[]
-#         new_p = (σ, ρ, β)
-#         new_prob = remake(prob, p = new_p)
-#         new_sol  = solve(new_prob, Tsit5(), reltol=1e-8, abstol=1e-8, saveat=0.1)
-#         # Update the line plot with the new solution
-#         lineplot[1][] = Point3f0.(new_sol[1, :], new_sol[2, :], new_sol[3, :])
-#     end
-
-#     # Connect each slider’s observable to update_plot!
-#     for slider in (σ_slider, ρ_slider, β_slider)
-#         on(slider.value) do _
-#             update_plot!()
-#         end
-#     end
-
-#     # --- Bonito App ---
-#     # Wrap the entire Makie figure (which includes the slider grid)
-#     # in a recorded DOM container so that slider changes propagate.
-#     app = App() do
-#         # By calling `Bonito.record_states` on the container,
-#         # all reactive states (including slider values) are captured.
-#         Bonito.record_states(sub, DOM.div(fig))
-#     end
-
-#     # Write the Bonito app to an HTML file.
-#     as_html(io, sub, app)
-# end
-
-
-
-# # plot 3 - differential equation
-# open("_posts/examples/diffeqviz/diffeq.html", "w") do io
-
-#     # function lorenz(du, u, p, t)
-#     #     x, y, z = u
-#     #     sigma, rho, beta = p
-#     #     du[1] = sigma * (y - x)
-#     #     du[2] = x * (rho - z) - y
-#     #     du[3] = x * y - beta * z
-#     # end
-
-#     # t_begin = 0.0
-#     # t_end = 10.0
-#     # tspan = (t_begin, t_end)
-#     # u_begin = [1.0, 0.0, 0.0]
-
-#     app = App() do session
-
-#         # create a slider
-#         a = Bonito.Slider(range(1, stop=3, length=3))
-
-#         # setup the ODE problem 
-#         # sol = lift(a) do a_val
-#         #     p = [10.0, 28.0, a_val]
-#         #     prob = ODEProblem(lorenz, u_begin, tspan, p)
-#         #     solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8)
-#         # end
-
-#         # Create a 3D makie plot with the solution
-#         # fig = Figure(size=(800, 400))
-#         # ax = Axis3(fig[1, 1])
-#         # x = Vector{Float64}(sol[][:][1, :])
-#         # y = Vector{Float64}(sol[][:][2, :])
-#         # z = Vector{Float64}(sol[][:][3, :])
-
-#         # fig, ax = meshscatter(
-#         #     x,
-#         #     y,
-#         #     z,
-#         #     markersize=1
-#         # )
-
-#         # Create a scatter plot
-#         fig, ax = meshscatter(rand(3, 100), markersize=markersize)
-
-#         # Return the plot and the slider
-#         return Bonito.record_states(session, DOM.div(fig, a))
-#     end
-# end
+struct C
+    m::Float64 # mass of the ball
+    g::Float64 # acceleration due to gravity
+    ρ::Float64 # density of the air
+    A::Float64 # cross-sectional area of the ball
+    C_D::Float64 # drag coefficient
+    r::Float64 # radius of the ball
+    ω₀::Float64 # initial angular velocity
+end
+
+struct IC
+    x::Float64 # initial x position
+    y::Float64 # initial y position
+    z::Float64 # initial z position
+    v_x::Float64 # initial x velocity
+    v_y::Float64 # initial y velocity
+    v_z::Float64 # initial z velocity
+end
+
+
+using DifferentialEquations
+
+C₁ = C(0.43, 9.81, 1.225, 0.013, 0.2, 0.11, 88)
+IC₁ = IC(0, 0, 0, 20, 0, 0)
+
+# Define the ODE function. The state vector u = [x, y, z, vx, vy, vz]
+function projectile!(ddu, du, u, p, t)
+    # Unpack state variables
+    x, y, z = u
+    vx, vy, vz = du
+
+    # Unpack parameters
+    m, g, ρ, A, C_D, r, ω₀ = p
+
+    # Define constant H = (ρ * A) / (2 * m)
+    H = (ρ * A) / (2 * m)
+
+    # Compute speed (magnitude of velocity)
+    v = sqrt(vx^2 + vy^2 + vz^2)
+
+    # Compute the decaying angular velocity ω(t)
+    ω = ω₀ * exp(-t/7)
+
+    # Compute lift coefficient, avoiding division by zero.
+    C_L = (v == 0.0 ? 0.0 : (ω * r / v))
+
+    # Compute acceleration components
+    ddu[1] = -v * H * (C_D * vx + C_L * vy)
+    ddu[2] = -v * H * (C_D * vy - C_L * vx)
+    ddu[3] = -v * H * (C_D * vz - g)
+end
+
+# Problem constants:
+# m: mass (kg), g: gravity (m/s²), rho: air density (kg/m³),
+# A: cross-sectional area (m²), C_D: drag coefficient,
+# r: ball radius (m), omega0: initial angular velocity (rad/s)
+p = [C₁.m, C₁.g, C₁.ρ, C₁.A, C₁.C_D, C₁.r, C₁.ω₀]
+
+# Initial conditions.
+# Here we assume the ball is kicked from the origin with an initial velocity.
+# For example: initial position (0, 0, 0) and initial velocity (30, 0, 30) m/s.
+dx₀ = [IC₁.v_x, IC₁.v_y, IC₁.v_z]
+x₀ = [IC₁.x, IC₁.y, IC₁.z]
+
+# Time span for the simulation
+tspan = (0.0, 50.0)
+
+# Define the ODE problem
+prob = SecondOrderODEProblem(projectile!, dx₀, x₀, tspan, p)
+sol = solve(prob, Tsit5())
+
+# load image 
+pitch = Makie.FileIO.load(expanduser("_posts/guides/interactive-blog/soccer_pitch.png"))
+
+# soccer pitch size 
+x_size = (0, 68)
+y_size = (0, 105)
+z_size = (0, 20)
+
+# make a plot for the initial conditions and a field 
+open(output_folder * "ode_ic.html", "w") do io
+    sub = Session(session)
+    app = App() do
+        fig = Figure(size=(500, 500))
+        ax = Axis3(fig[1, 1])
+        p0 = Point3f(x₀)
+        v0 = Vec3f(dx₀)
+
+        xlims!(ax, x_size)
+        ylims!(ax, y_size)
+        zlims!(ax, z_size)
+        arrows!(ax, [p0], [v0], color=:red)
+
+        # set picture as xy plane
+        surface!(ax, [0, 0, 68, 68], [0, 105, 105, 0], [0, 0, 0, 0], 
+                    color=:lightgreen, transparency=true)
+
+        fig
+    end
+    as_html(io, sub, app)
+end