Example program of XLPack for Python and XLPack for Matplotlib: (3) Animation (Arenstorf orbits)
Consider two bodies with masses \(\mu\) and \(\mu’ (= 1 – \mu)\) undergoing planar orbital motion, and a third body whose mass is negligible compared to the other two, moving in the same plane. The equations of motion are given as follows (reference: E. Hairer et al., “Solving Ordinary Differential Equations I” (1987, 1993), Springer-Verlag).\[
\begin{align}
& y_1” = y_1 + 2y_2′ – \mu'(y_1 + \mu)/D_1 – \mu(y_1 – \mu’)/D_2 \\
& y_2” = y_2 – 2y_1′ – \mu’y_2/D_1 – \mu y_2/D_2 \\
& ただし, \\
& D_1 = ((y_1 + \mu)^2 + y_2^2)^{3/2} \\
& D_2 = ((y_1 – \mu’)^2 + y_2^2)^{3/2} \\
& \mu = 0.012277471, \mu’ = 1 – \mu \\
\end{align}
\] It is known that using the following initial conditions yields a periodic solution at \(t = 17.0652165601579625588917206249\).
\[
\begin{align}
& y_1(0) = 0.994, y_1′(0) = 0 \\
& y_2(0) = 0, y_2′(0) = -2.00158510637908252240537862224 \\
\end{align}
\] This is called as the Arenstorf orbits.
Here, we compute the numerical solution over one period and display it as an animation.
1. Overall structure
The main program is as follows.
Const Tend = 17.065216560158
Sub Main()
Dim Xt() As Double, Yt() As Double
Dim Dt As Double, Nt As Long, Info As Long
Dt = 0.1
Nt = Int(Tend / Dt) + 1
ReDim Xt(Nt), Yt(Nt)
Call ODE(Tend, Dt, Nt, Xt(), Yt(), Info)
Call Animate(Nt, Xt(), Yt(), Dt)
End SubThe subroutine ODE() solves the differential equation. The obtained solution is displayed by the subroutine Animate(). In Animate(), animation is performed using XLPack for Python or XLPack for Matplotlib.
The execution result will be as follows.
2. Solutions to differential equations
Since this is a second-order ordinary differential equation, we transform it into a system of first-order ordinary differential equations.
\[
\begin{align}
& Y(0) = y_1 \\
& Y(1) = y_2 \\
& Y(2) = y_1′ \\
& Y(3) = y_2′ \\
& Y(0)’ = y_1′ = Y(2) \\
& Y(1)’ = y_2′ = Y(3) \\
& Y(2)’ = y_1” = Y(0) + 2Y(3) – μ'(Y(0) + μ)/D1 – μ(Y(0) – μ’)/D2 \\
& Y(3)’ = y_2” = Y(1) – 2Y(2) – μ’Y(1)/D1 – μY(1)/D2 \\
\end{align}
\]
We solve this using a method for solving first-order systems of ordinary differential equations (Derkfa from XLPack).
Const Xinit = 0.994
Sub F(N As Long, T As Double, Y() As Double, Yp() As Double)
Const Mu = 0.012277471, Mup = 1 - Mu
Dim D1 As Double, D2 As Double
Yp(0) = Y(2)
Yp(1) = Y(3)
D1 = ((Y(0) + Mu) ^ 2 + Y(1) ^ 2) ^ (3# / 2#)
D2 = ((Y(0) - Mup) ^ 2 + Y(1) ^ 2) ^ (3# / 2#)
Yp(2) = Y(0) + 2 * Y(3) - Mup * (Y(0) + Mu) / D1 - Mu * (Y(0) - Mup) / D2
Yp(3) = Y(1) - 2 * Y(2) - Mup * Y(1) / D1 - Mu * Y(1) / D2
End Sub
Sub ODE(Tend As Double, Dt As Double, Nt As Long, Xt() As Double, Yt() As Double, Info As Long)
Const N = 4, Mode = 3, Tol = 0.0000000001 '1.0e-10
Dim Y(N - 1) As Double, T As Double, Tout As Double
Dim RTol(0) As Double, ATol(0) As Double, I As Integer
RTol(0) = Tol
ATol(0) = Tol
T = 0
Y(0) = Xinit: Y(1) = 0: Y(2) = 0: Y(3) = -2.00158510637908
Info = 0
Xt(0) = Y(0)
Yt(0) = Y(1)
For I = 0 To Nt - 1
Tout = (I + 1) * Dt
If Tout > Tend Then Tout = Tend
Call Derkfa(N, AddressOf F, T, Y(), Tout, Tend, RTol(), ATol(), Mode, Info)
If Info < 0 Or Info > 10 Then Exit Sub
Xt(I + 1) = Y(0)
Yt(I + 1) = Y(1)
Next
End SubThe obtained solution is saved in Xt() and Yt().
3. Animation of the Solution
The obtained solution can be directly visualized using Python code written with using XLPack for Python.
Since array data cannot be passed directly, the data is first written to a worksheet and then read on the Python side using xlc_get().
Sub Animate(Nt As Long, Xt() As Double, Yt() As Double, Dt As Double)
Const Col = 2, Row = 6
Dim I As Long
For I = 0 To Nt
Cells(Row + I, Col) = Xt(I)
Cells(Row + I, Col + 1) = Yt(I)
Next
'-- Python codes
PY "import matplotlib.pyplot as plt"
PY "import matplotlib.animation as animation"
PY "from XLPackPy import *"
PY "fig = plt.figure()"
PY "ax = fig.gca()"
PY "ax.grid()"
PY "ax.plot(0, 0, '-o', color = 'orange')"
PY "ax.plot(1, 0, '-o', color = 'yellow')"
PY "ax.set_xlim(-1.5, 1.5)"
PY "ax.set_ylim(-1.5, 1.5)"
PY "xt = xlc_get(5, 1," + Str(5 + Nt) + ")"
PY "yt = xlc_get(5, 2," + Str(5 + Nt) + ")"
PY "xtrace = []; ytrace = []; artist = []"
PY "for i in range(len(xt)):" + vbNewLine + _
" a_x = [xt[i]]; a_y = [yt[i]]" + vbNewLine + _
" a, = ax.plot(a_x, a_y, 'o-', color = 'b')" + vbNewLine + _
" xtrace.append(xt[i]); ytrace.append(yt[i])" + vbNewLine + _
" trace, = ax.plot(xtrace, ytrace, ':c', lw=1)" + vbNewLine + _
" artist.append((a, trace))"
PY "anim = animation.ArtistAnimation(fig, artist, interval = 100, blit = True, repeat = False)"
PY "plt.show()"
End SubUsing XLPack for Matplotlib, the same processing can be written entirely in VBA. Internally, the processing is carried out by calling Python.
Dim Plt As New PyPlot
Dim Animation As New Animation
Sub Animate(Nt As Long, Xt() As Double, Yt() As Double, Dt As Double)
Dim Ex(0) As Double, Ey(0) As Double, Mx(0) As Double, My(0) As Double
Dim X(0) As Double, Y(0) As Double, I As Long
Dim Fig As Figure, Ax As Axs, Artist() As PyObject
Dim Emarker As PyObject, Mmarker As PyObject
Dim Anim As ArtistAnimation
Ex(0) = 0: Ey(0) = 0
Mx(0) = Xinit: My(0) = 0
ReDim Artist(Nt, 1)
Set Fig = Plt.Figure()
Set Ax = Fig.Gca()
Call Ax.Grid
Call Ax.Set_xlim(-1.5, 1.5)
Call Ax.Set_ylim(-1.5, 1.5)
Set Emarker = Ax.Plot(1, Ex(), Ey(), "", "color=orange,marker=o")
Set Mmarker = Ax.Plot(1, Mx(), My(), "", "color=yellow,marker=o")
For I = 0 To Nt
X(0) = Xt(I): Y(0) = Yt(I)
Set Artist(I, 0) = Ax.Plot(1, X(), Y(), "", "color=b,marker=o")
Set Artist(I, 1) = Ax.Plot(I + 1, Xt(), Yt(), ":c", "lw=1")
Next
Set Anim = Animation.ArtistAnimation(Fig, Nt + 1, 2, Artist(), "interval =" + Str(Dt * 1000) + ", repeat = False")
Call Plt.Show
End Sub
