8  Using packages

Prerequisites

Before reading this chapter, you are recommended to have read Chapter 2.

8.1 What is a package?

8.1.1 Modules

The programs we’ve written so far (and indeed, will continue to write for the purposes of this guide) are relatively short. However, in real applications, code may span thousands of lines, and it can be useful to group up sections of code into large chunks, to be used when necessary. This is the idea behind the module.

Modules give a name to encompass all of the functions, types, etc. contained within. This means the whole lot can be included in one go, and no part will be included without its dependents. Code is put into a module simply by surrounding it by a module-end block as the example below shows, with the name of the module going after module:

module MyChessEngine

const PIECES = (:pawn, :bishop, :knight, :rook, :queen, :king)
const COLOURS = (:white, :black)

struct ChessPiece
    piece::Symbol
    colour::Symbol
    row::Int64
    file::Int64
end

global turn = :white

[...]

end

Convention

Modules are named in upper camel-case, just like types. The code inside them should not be indented like other blocks, as there’s no need to indent the whole file.

Good examples of modules are Base and Core. Core is the module of crucial functionality that has to be written in another programming language (specifically, in C), while Base comprises the rest of the inbuilt functionality that Julia has, written entirely in Julia. These are special, in that the code within them is automatically accessible whenever you start Julia. Other modules exist too in the standard library (that is, the set of modules that comes preinstalled with Julia), but you have to tell Julia to specifically make them usable, with the keyword import. For instance, we can use the LinearAlgebra module (which adds more options when calculating with Arrays) as follows:

import LinearAlgebra

Then, we can use the functions from LinearAlgebra, by referring to them by their module, followed by a ., and then their name:

A = [[1, 2] [3, 4]]

2×2 Matrix{Int64}:
 1  3
 2  4

LinearAlgebra.det(A)

-2.0

More commonly, instead of using import, we use using, which allows for exported functions, types, variables, etc. to be referred to without the module name. The exported names are listed within the module’s code, after the keyword export, and are usually placed at the top of the module’s code for reference. det is exported by LinearAlgebra, hence:

using LinearAlgebra

det(A)

-2.0

The downside to using is that exported names from different modules could clash with each other. This can be remedied by simply including the module name as before, just as you would have done with import. Moreover, import can be used to select individual functions or types (which we used in Chapter 7), whereas using has to select entire modules.

8.1.2 Packages

While Julia’s Core and Base modules are extensive, and various modules exist in the standard library, they doesn’t necessarily do everything that you might want them to do. This is deliberate, after all every time you open a new REPL instance, both need to be loaded in, and it would be redundant to load in millions of lines of code which would never be any use to you. Instead, Core is as minimal as possible, and Base is written to be generally applicable, so for more specific applications, we need to turn to packages.

Packages are, for the most part, modules that others have developed for a specific purpose, which can be installed into Julia from the internet. They are to Julia as expansion packs are to games (video or board, take your pick), or as subscription services are to television, although one handy difference is that Julia’s packages are free! Some packages may just add functions and types, like Distributions.jl which allows for working with statistical distributions. More complicated packages can do much more, such as Pluto.jl which provides an entirely new environment to program Julia in.

Convention

Packages are often named online (such as in the GitHub repositories that contain them) with the file extension .jl, since that is the file name used for some or all of its code, and specifies that these files are for use with Julia. Within Julia, we drop the .jl, since the underlying module doesn’t have that as part of its name. When referring to the package, both forms are used pretty interchangeably.

We won’t consider here how modules and packages can be written, uploaded, and distributed, as this is unnecessary for small, beginner projects. However, being able to use packages that already exist ranges from helpful to essential in such projects, some examples of which are CSV (detailed in Chapter 10) and Plots (the focus of Chapter 11). We’ll also use a couple of other packages in examples, including NLsolve later in this chapter, and Colors, ImageShow, and ImageIO from Chapter 9 to display Julia and Mandelbrot sets.

8.2 Installing packages with Pkg

Pkg, short for package, is one of the modules available in Julia’s standard library, and is the package manager of Julia, meaning it is the interface through which you’ll most likely install, modify, and remove packages. Like any other module, we can use it through syntax like using Pkg, but it is important enough to have its own special place in the REPL.

If we open the REPL and press ], then this enters package mode (much as ? enters help mode), with the prompt changing to blue, and being the Julia version number followed by pkg. For example, in Julia 1.9, it looks a little like:

(@v1.9) pkg>

There are three commands in particular that we’ll focus on, add, update, and remove. You can probably guess what each of these do, but lets see how to use them anyway.

8.2.1 add

To install a new package, we need to use add. Simply following add by the name of the package will do just that. As an example, we’ll use the Colors package, which adds various support for dealing with colours and is widely used by other packages like Plots:

(@v1.9) pkg> add Colors

Warning

Beware of different spellings! If you’re used to ‘colour’ with a ‘u’, instead of ‘color’ without, then you might be greeted with error messages when trying to use the Colors package. Across programming languages, this is a near universal standard with the word color/colour in particular, but at the time of writing, one package does exist for Julia with the ‘u’ spelling, called PerceptualColourMaps. Similar issues will arise in package, module, function, and type names for other variably spelt words.

Note

Pkg can just be used like any other module, of course. For instance, using Pkg followed by Pkg.add("Colors") in the normal REPL mode has the same effect as add Colors does in package mode.

Colors is a relatively small package, so this should be very quick. Once it’s done, you’ll see a list of all the dependencies, that is the other packages and modules that the Colors package relies on to work, all of which are automatically downloaded if they aren’t already.

Now, it’s ready to use, so we can exit Pkg mode with ← Backspace, deleting the invisible ] at the start of the line to return to the normal REPL (notice how, different to help mode, the REPL stays in Pkg mode after running a line, since you may want to do multiple actions in Pkg mode in a row). Then, as before with modules, either using Colors or import Colors can be entered to load the additional functionality that the Colors package provides.

This begs the question: how do we actually find the packages that we need to use? There are several ways:

• The most obvious, and often easiest way, is searching online for ‘Julia package to […]’. This works most of the time

• Examples of Julia code frequently use packages, as we have done in this book. Anywhere where such examples can be found, packages can be found

• The website JuliaHub lists all the packages that have been made publically available and are accessible through Pkg. In particular, the search at https://juliahub.com/ui/Packages is useful in finding the right package

• Tab-completion works in Pkg mode. For example, typing add C and pressing Tab ⇆ lists out all the packages that Pkg knows about starting with C (note this is case-sensitive)

8.2.2 update

Just like Julia, most packages are being continually updated and improved over time. These updates are not installed automatically, and so weeks or months after first using them, you may find your packages to be outdated. Luckily, updating all the packages is as simple as one Pkg command:

(@v1.9) pkg> update

You can update specific packages if you wish by following this with the name of a package, and a version number if you want to be even more specific.

8.2.3 remove

If there’s a package you’re no longer using, you may want to delete it. This is done simply by the remove command, followed by the package name:

(@v1.9) pkg> remove Colors

Note

When Pkg installs packages, they will only be available for that environment alone. An environment is a folder in your file system (specifically, it’s a subdirectory of ~/.julia/environments/), containing configuration files called Project.toml and Manifest.toml that define which packages are available for use. If you’re working entirely within the REPL, or in an IDE which runs code in the REPL, the default environment has the name of the Julia version that you’re running, which is why that is shown in the Pkg mode prompt.

However, be aware that using Pluto notebooks creates a separate environment, which therefore has its own set of installed packages. Packages also work differently in Pluto, as it uses its own package manager, which doesn’t require direct use of Pkg, and Project.toml and Manifest.toml exist within the raw .jl file that Pluto produces. This allows for using and import to be used straight away, with packages downloaded as and when they are needed.

Also, updating Julia will create a new environment for the new version, so packages often need to be reinstalled afresh(depending on how you update Julia). Pkg provides further tools for working with different environments, but we’ll leave these alone.

8.3 Example: Lagrange points with NLsolve

The real use of packages, in particular from a scientific standpoint, comes in providing solutions to very specific use cases that others have come across before. In the following example, we’ll derive some equations that we can’t solve, but fortunately for us, a package will exist to find the solution. While this example itself may be quite specific, the scenario of having such a particular problem and requiring a particular package to solve it is not at all unusual.

8.3.1 What are Lagrange points?

Note

The method for finding the position of Lagrange points in this section mirrors a document written by Neil J. Cornish of NASA, to coincide with the 1998 launch of the WMAP probe to L₂.

One of the oldest problems in dynamics is the three-body problem, which asks how the orbit of three bodies in space evolves through time, when each is pulling gravitationally on the others. In its most general form, it is analytically unsolveable, and can only be approximated numerically, although there are some special cases where exact results can be deduced.

In particular, we’ll focus on a simplified version called the restricted three-body problem, where one of the three bodies is far smaller than the other two, so much so that we can ignore the impact of its gravitational pull on them. In this scenario, five remarkable points of stability exist, called the Lagrange points, where this third body can lie indefinitely. These are laid out approximately as follows:

Figure 8.1: The Lagrange points of the Sun-Earth system

• L₁ lies between the two larger bodies
• L₂ lies beyond the smaller of the two larger bodies
• L₃ lies beyond the larger of the two larger bodies
• L₄ lies in the path of the smaller of the two larger bodies, approximately 60 degrees ahead of it
• L₅ lies in the path of the smaller of the two larger bodies, approximately 60 degrees behind it

In this illustration, we’ve used the Sun as the largest body, and the Earth as the smaller large body. Indeed, this is quite an apt choice, since the Lagrange points of the Earth-Sun system are the most useful ones to scientists at the moment, most notably perhaps in the case of the James Webb Space Telescope which has orbited in and around L₂ since 2022. Moreover, while orbits around L₁, L₂, and L₃ are unstable and require periodic correction to keep satellites in place, orbits around L₄ and L₅ do not, and so asteroids called trojans have been found orbiting at L₄ (in fact, many more trojans exist at the L₄ and L₅ points of the Sun-Jupiter system, which is where they were first discovered).

In that case, why doesn’t the Earth run into the trojans at L₄? What isn’t really shown on this diagram is the fact that the whole system is rotating around a common centre of mass, somewhere between the Earth and the Sun. Obviously, the fact that the Earth orbits the Sun shouldn’t be too surprising (sorry Pope Urban VIII), but this is better described as both the Earth and the Sun orbiting each other, with the point that they orbit around being much closer to the Sun since the Sun is much more massive. Additionally, the Lagrange points move along with this orbit, and this will be crucial to us being able to find them.

8.3.2 Deriving an equation for the Lagrange points

Warning

For brevity, we won’t explain everything about what we’re doing here, and you’re encouraged to skip straight to the programming if you’d prefer!

Before we get to finding them, we need to do a bit of mathematics to give us some equations to solve.

We’ll call the masses of the two larger bodies \(M_1\) and \(M_2\), and the mass of the smaller body \(m\) (although as we’ll see, this won’t actually matter). With the centre of mass (and hence the centre of rotation) of the system at \(\mathbf{0}\), we’ll let the larger bodies lie at \(\mathbf{r}_1\) and \(\mathbf{r}_2\) respectively. By choice of appropriate basis vectors and unit of length, we can take \(\mathbf{r}_1 = r_1 \hat{\mathbf{x}}\) and \(\mathbf{r}_2 = r_2 \hat{\mathbf{x}}\), with \(|r_1 - r_2| = 1\). Then, as we’re ignoring any effect of the smaller mass on these two, their positions can be solved for as the solution of the two-body problem, namely: \[ \mathbf{r}_1 = -\alpha \hat{\mathbf{x}}, \quad \mathbf{r}_2 = (1 - \alpha) \hat{\mathbf{x}} \qquad \text{for } \alpha = \frac{M_2}{M_1 + M_2} \]

We also need to consider how fast the orbit of these two bodies is, and in particular find its angular velocity, which we’ll call \(\Omega\). Since the centripetal force on each body consists solely of the gravitational force from the other, we find that: \[ - M_1 \Omega^2 \mathbf{r}_1 = - \frac{G M_1 M_2}{|\mathbf{r}_1 - \mathbf{r}_2|^3} (\mathbf{r}_1 - \mathbf{r}_2) \]

We can substitute in the values we already know to give: \[ M_1 \Omega^2 \frac{M_2}{M_1 + M_2} \hat{\mathbf{x}} = G M_1 M_2 \hat{\mathbf{x}} \]

hence: \[ \Omega^2 = G(M_1 + M_2) \tag{8.1}\]

Now, we can concentrate on the smaller body. We’ll work in a rotating frame of reference relative to the orbit of the larger bodies (i.e. at angular velocity \(\Omega\) around \(\mathbf{0}\)), which means that the Lagrange points will be stationary and so easier to find, however this does mean that we’ll have to add fictitious forces into the mix to compensate. The mass \(m\) will lie in the plane of the orbit at a position \(\mathbf{r} = x \hat{\mathbf{x}} + y \hat{\mathbf{y}}\), with the forces acting upon it as follows:

• The gravitational force from \(M_1\) on \(m\) is given by: \[ F_1 = - \frac{G M_1 m}{|\mathbf{r} - \mathbf{r}_1|^3} (\mathbf{r} - \mathbf{r}_1) = - \frac{M_1 m}{((x + \alpha)^2 + y^2)^\frac{3}{2}} ((x + \alpha) \hat{\mathbf{x}} + y \hat{\mathbf{y}}) \]

• The gravitational force from \(M_2\) on \(m\) is similarly: \[ F_2 = - \frac{G M_2 m}{|\mathbf{r} - \mathbf{r}_2|^3} (\mathbf{r} - \mathbf{r}_2) = - \frac{M_2 m}{((x + \alpha - 1)^2 + y^2)^\frac{3}{2}} ((x + \alpha - 1) \hat{\mathbf{x}} + y \hat{\mathbf{y}}) \]

• The Coriolis force is a fictitious force defined in terms of the velocity of \(m\) and the angular velocity vector \(\mathbf{\Omega} = \Omega \hat{\mathbf{z}}\): \[ F_{Cor} = - 2m \mathbf{\Omega} \times \dot{\mathbf{r}} = - 2m\Omega (\dot{x} \hat{\mathbf{y}} - \dot{y} \hat{\mathbf{x}}) \]

• The centrifugal force is a fictitious force that balances the inherent centripetal force needed to stay in place in a rotating frame of reference \[ F_{Cent} = - m \mathbf{\Omega} \times (\mathbf{\Omega} \times \mathbf{r}) = - m \Omega^2 (- x \hat{\mathbf{x}} - y \hat{\mathbf{y}}) \]

Then, force is mass times acceleration, with the total force being the sum of the above components: \[ m \ddot{\mathbf{r}} = F_1 + F_2 + F_{Cor} + F_{Cent} \]

which simplifies (using Equation 8.1) to: \[ \begin{aligned} \ddot{\mathbf{r}} &= \left[ \Omega^2 x - \frac{\Omega^2 (1 - \alpha) (x + \alpha)}{((x + \alpha)^2 + y^2)^\frac{3}{2}} - \frac{\Omega^2 \alpha (x + \alpha - 1)}{((x + \alpha - 1)^2 + y^2)^\frac{3}{2}} + 2 \Omega \dot{y} \right] \hat{\mathbf{x}} \\ &\quad + \left[ \Omega^2 y - \frac{\Omega^2 (1 - \alpha) y}{((x + \alpha)^2 + y^2)^\frac{3}{2}} - \frac{\Omega^2 \alpha y}{((x + \alpha - 1)^2 + y^2)^\frac{3}{2}} - 2 \Omega \dot{x} \right] \hat{\mathbf{y}} \end{aligned} \]

At a Lagrange point, both the velocity \(\dot{\mathbf{r}}\) and acceleration \(\ddot{\mathbf{r}}\) are zero. If we substitute this into the above, and cancel the factor of \(\Omega^2\), this means that the positions of the Lagrange points are the solutions \(\mathbf{r} = x \hat{\mathbf{x}} + y \hat{\mathbf{y}}\) to the vector equation: \[ \mathbf{F}(\mathbf{r}) = \begin{pmatrix} x - \frac{(1 - \alpha) (x + \alpha)}{((x + \alpha)^2 + y^2)^\frac{3}{2}} - \frac{\alpha (x + \alpha - 1)}{((x + \alpha - 1)^2 + y^2)^\frac{3}{2}} \\ y - \frac{(1 - \alpha) y}{((x + \alpha)^2 + y^2)^\frac{3}{2}} - \frac{\alpha y}{((x + \alpha - 1)^2 + y^2)^\frac{3}{2}} \end{pmatrix} = \mathbf{0} \tag{8.2}\]

8.3.3 Finding the Lagrange points with NLsolve

Prerequisites

For this section, you should know how to write a custom function (Chapter 6), and various Array operations (Chapter 9).

Let’s quickly review the equation that we’ve just found:

• The masses of the larger bodies are \(M_1\) and \(M_2\), with the crucial parameter \(\alpha\) defined by \(\alpha = \frac{M_2}{M_1 + M_2}\), i.e. the proportion of the total mass concentrated in the second large body

• The two larger bodies \(M_1\) and \(M_2\) lie at \((-\alpha, 0)\) and \((1-\alpha, 0)\) respectively

• The smaller body \(m\) lies at \((x,y)\). To find a Lagrange point, we need to solve Equation 8.2 for \(x\) and \(y\)

Unfortunately, Equation 8.2 looks horrific, and no amount of algebraic arrangement will spit us out a solution. Instead, we’ll have to solve it numerically, and this is where the package NLsolve comes in. NLsolve stands for non-linear solve, and its purpose is the solve equations just like ours. Let’s download it and load it in as we would with any other package:

# Alternative to Pkg mode in the REPL
# Simply use the Pkg module as you would any other
using Pkg
Pkg.add("NLsolve");

using NLsolve

Note

As with many packages, once NLsolve is loaded in, typing NLsolve into help mode in the REPL will give a run down of the full capabilities of the package, but we’ll just stick to what we’ll need for now.

First, we’ll need to write a mutating function to represent our vector function \(\mathbf{F}\). This must take two inputs, with the first being a Vector that we’ll mutate to the value of \(\mathbf{F}(\mathbf{r})\), and the second being a Vector representing the input \(\mathbf{r}\). Since we want to mutate the output Vector, we need to be sure to assign its values element by element, in order to not redefine the variable name:

function evalF!(F, r)
    x, y = r
    F[1] = x -
        (1 - α)*(x + α) / ((x + α)^2 + y^2)^(3/2) -
        α*(x + α - 1) / ((x + α - 1)^2 + y^2)^(3/2)
    F[2] = y -
        (1 - α)*y / ((x + α)^2 + y^2)^(3/2) -
        α*y / ((x + α - 1)^2 + y^2)^(3/2)
end

The problem with this function is that it depends on the parameter \(\alpha\), which wouldn’t be defined here. Therefore, we’ll instead write a function with an input of the parameter \(\alpha\), and return evalF! as its output:

function createevalF(α)
    return function evalF!(F, r)
        x, y = r
        F[1] = x -
            (1 - α)*(x + α) / ((x + α)^2 + y^2)^(3/2) -
            α*(x + α - 1) / ((x + α - 1)^2 + y^2)^(3/2)
        F[2] = y -
            (1 - α)*y / ((x + α)^2 + y^2)^(3/2) -
            α*y / ((x + α - 1)^2 + y^2)^(3/2)
    end
end

createevalF (generic function with 1 method)

To use this with the NLsolve package, we’ll need its nlsolve function. A variety of possibilities for inputs are allowed, but we’ll use the following:

• The first input needs to be the mutating function evalF! that we’ve just constructed

• The second input needs to be a starting point, ideally near to the solution. Because the Lagrange points are spread apart nicely, and the function \(\mathbf{F}\) gets very large in absolute value near \(\mathbf{r}_1\) and \(\mathbf{r}_2\), we’ll find conveniently that the starting points shown in Table 8.1 will work flawlessly. To actually use these in Julia, we’ll want to use a Vector of Float64s to represent them (e.g. [0.0, 0.0] for \((0, 0)\))

Table 8.1: The starting points we’ll using for each Lagrange point with nlsolve

Lagrange point	Starting position
L₁	\((0, 0)\)
L₂	\((1, 0)\)
L₃	\((-1, 0)\)
L₄	\((0, 1)\)
L₅	\((0, -1)\)

• The third input will be a keyword argument autodiff = :forward. This tells NLsolve to use automatic differentiation to find derivatives of our function \(\mathbf{F}\) where it needs them, which is particularly well suited to our function, as well as to the multiple dispatch system of Julia

Let’s try this out with a test value for \(\alpha\), say α = 0.1:

solution = nlsolve(createevalF(0.1), [0.0, 0.0], autodiff = :forward)

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [0.0, 0.0]
 * Zero: [0.609035110084086, 0.0]
 * Inf-norm of residuals: 0.000000
 * Iterations: 14
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 15
 * Jacobian Calls (df/dx): 15

This gives us quite a lot of information about the solving process, but we’re only really interested in the answer in this case, which can be found by the zero field:

solution.zero

2-element Vector{Float64}:
 0.609035110084086
 0.0

Now that we’ve worked out how to use NLsolve, we want to choose a list of values of \(\alpha\) to try. We’ll use:

# Can assume M₁ ≥ M₂, so no point using α > 0.5
αs = 0.001:0.001:0.5

0.001:0.001:0.5

Then, it’s just a matter of running NLsolve on all the possible values of \(\alpha\) for each of the five points in turn. We could do this with a for loop, or with array comprehension followed by a bit of re-Array-ngement:

const STARTING_POINTS = ([0.0, 0.0], [1.0, 0.0], [-1.0, 0.0], [0.0, 1.0], [0.0, -1.0])

lagrangepoints = [
    nlsolve(createevalF(α), r₀, autodiff = :forward).zero
    for α ∈ αs, r₀ ∈ STARTING_POINTS
] |> stack |> (A -> permutedims(A, (3, 2, 1)))

5×500×2 Array{Float64, 3}:
[:, :, 1] =
  0.931287   0.913219   0.900374  …   0.00282353   0.00141177   0.0
  1.06992    1.08786    1.10029       1.1991       1.19875      1.19841
 -1.00042   -1.00083   -1.00125      -1.19771     -1.19806     -1.19841
  0.499      0.498      0.497         0.002        0.001        0.0
  0.499      0.498      0.497         0.002        0.001        0.0

[:, :, 2] =
  0.0        0.0        0.0       …   0.0        0.0        0.0
  0.0        0.0        0.0           0.0        0.0        0.0
  0.0        0.0        0.0           0.0        0.0        0.0
  0.866025   0.866025   0.866025      0.866025   0.866025   0.866025
 -0.866025  -0.866025  -0.866025     -0.866025  -0.866025  -0.866025

The Array that we’re left with (lagrangepoints) has three dimensions. Specifying the index of the first chooses the number of the Lagrange point, while the index of the second dimension chooses the value of \(\alpha\) corresponding to that index in αs, and lastly the index of the second dimension is either 1 or 2 corresponding to \(x\) and \(y\) respectively. For instance, the \(x\)-coordinate of L₃ in a system with \(\alpha = 0.056\) is:

lagrangepoints[3, 56, 1]

-1.023323464216801

While we have got \(x\) and \(y\)-coordinates, the \(y\)-coordinate isn’t actually very interesting. From Figure 8.1, we can see that L₁, L₂, and L₃ will have \(y\)-coordinate exactly \(0\), which they do:

lagrangepoints[1:3, :, 2]

3×500 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0

Meanwhile, L₄ and L₅ are exactly the third vertex of an equilateral triangle with vertices at the two larger masses, so their \(y\)-coordinates are exactly \(\frac{\sqrt{3}}{2}\) and \(-\frac{\sqrt{3}}{2}\) respectively, again borne out by our calculations:

lagrangepoints[4:5, :, 2]

2×500 Matrix{Float64}:
  0.866025   0.866025   0.866025  …   0.866025   0.866025   0.866025
 -0.866025  -0.866025  -0.866025     -0.866025  -0.866025  -0.866025

 √3/2

0.8660254037844386

Instead, we’ll plot the \(x\)-coordinate as \(\alpha\) changes:

Figure 8.2: The \(x\)-coordinates of the Lagrange points as α varies

Note

In Figure 8.2, we’ve extended our data to \(\alpha > 0.5\). As this happens, the masses \(M_1\) and \(M_2\) are switch (since \(M_1\) is always the largest), and therefore so do the Lagrange points (which are defined relative to which body is larger). We’ve accounted for this, so the Lagrange point referred to as L₂ is actually L₃ when \(\alpha > 0.5\), and vice versa, and the same for L₄ and L₅.

Note

In the NASA document by Neil J. Cornish, the stability of the Lagrange points is already examined. We’ve already said that for the Sun-Earth system, orbits around L₁, L₂, and L₃ are unstable, while orbits around L₄ and L₅ are stable. Actually, the stability of L₄ and L₅ is depends on \(\alpha\) being small enough, specifically less than the bizzare quantity \(\frac{2}{27 + \sqrt{621}}\). Cornish proves this analytically, but we could do the same numerically in Julia by simulating orbits, for instance with the DifferentialEquations package.