Constellations in Mathematica

Posted on 4 December 2023 by John

Mathematica has data on stars and constellations. Here is Mathematica code to create a list of constellations, sorted by the declination (essentially latitude on the celestial sphere) of the brightest star in the constellation.

constellations = EntityList["Constellation"]
sorted = SortBy[constellations, -#["BrightStars"][[1]]["Declination"] &]

We can print the name of each constellation with

Map[#["Name"] &, sorted]

This yields

{"Ursa Minor", "Cepheus", "Cassiopeia", "Camelopardalis", 
…, "Hydrus", "Octans", "Apus"}

We can print the name of the constellation along with its brightest star as follows.

Scan[Print[#["Name"], ", " #["BrightStars"][[1]]["Name"]] &, sorted]

This prints

Ursa Minor, Polaris
Cepheus, Alderamin
Cassiopeia, Tsih
Camelopardalis, β Camelopardalis
…
Hydrus, β Hydri
Octans, ν Octantis
Apus, α Apodis

Mathematica can draw star charts for constellations, but when I tried

Entity["Constellation", "Orion"]["ConstellationGraphic"]

it produced extraneous text on top of the graphic.

A disk around Paris

Posted on 1 November 2023 by John

The other day I saw an image of a large disk centered on Paris subjected to the Mercator projection. I was playing around in Mathematica and made similar images for different projections. Each image below is a disk of radius 4200 km centered on Paris (latitude 49°, longitude 2°).

All images were produced with the following Mathematica code, changing the GeoProjection argument each time.

    GeoGraphics[GeoDisk[GeoPosition[{49, 2}],
       Quantity[4200, "Kilometers"] ],
       GeoProjection -> "...", 
       GeoRange -> "World"]

Robinson projection

    … GeoProjection -> "Robinson", …

Robinson projection

Winkel-Snyder projection

    … GeoProjection -> "WinkelSnyder", …

Winkel-Snyder projection

Orthographic projection

    … GeoProjection -> "Orthographic", …

Orthographic projection

Lambert Azimuthal projection

    … GeoProjection -> "LambertAzimuthal", …

Lambert Azimuthal projection

Peirce Quincuncial projection

    … GeoProjection -> "PeirceQuincuncial", …

Peirce Quincuncial projection

This last projection has some interesting mathematics and history behind it. See this post for the backstory.

Curvature: principal, Gauss, and mean

Posted on 24 October 2023 by John

This post will compute the center of curvature for an object described in the previous post. In order to do that we first need to describe principle curvature and Gauss curvature, and we’ll throw in mean curvature while we’re at it.

Let S be a surface sitting in three dimensional space. No need for more generality to get where we want to go. At any point p on S we can draw curves on the S through p and compute the curvature at p. The curvature is the reciprocal of the radius of the kissing circle.

If we draw curves through p in every direction, some may have larger or smaller curvature than others. Let k₁ and k₂ be the minimum and maximum curvatures at p. These re the principal curvatures. The product k₁ k₂ of the principle curvatures is the Gaussian curvature and their average ( k₁ + k₂)/2 is the mean curvature. Incidentally, when principle curvatures are not equal. the directions in which the curvature is minimized and maximized are orthogonal.

In the previous post I said that the center of curvature for the surface with equation

is finite because the curvature is always positive. In particular, we wanted to know the center of curvature at the bottom, where x = y = 0 and z = −h.

The calculation to get there is messy, but the end result is simple: the principle curvatures are equal by symmetry, and both equal 1 − s². Therefore the center of curvature is at z = 1/(1 − s²).

Calculation details

The following Mathematica code calculates the (signed) curvature of a curve of the form F(x, y) = 0.

k[f_, x_, y_] := (D[f[x, y], y]^2 D[f[x, y], {x, 2}] - 
    2 D[f[x, y], x] D[f[x, y], y] D[f[x, y], {x, 1}, {y, 1}] + 
    D[f[x, y], x]^2 D[f[x, y], {y, 2}]) / (D[f[x, y], x]^2 + 
     D[f[x, y], y]^2)^(3/2)

Define

g[x_, y_] := x^2 + y^2 - s^2 x^2 y^1 - 1

and replace y with z/h. When we evaluate the curvature at x = 0 and simplify

Simplify[k[g, x, y] /. { y -> z/h, x -> 0}, Assumptions -> {h > 0}]

we get

(h – s² z) / |z|.

When z = −h we find that the unsigned curvature is 1 − s². In the previous post we assumed h > 1, but the calculation above shows that if h < 1 it’s possible for the curvature to be 0. (Recall s must be between 0 and 1.)

Checksum polynomials

Posted on 6 September 2023 by John

A large class of checksum algorithms have the following pattern:

Think of the bits in a file as the coefficients in a polynomial P(x).
Divide P(x) by a fixed polynomial Q(x) mod 2 and keep the remainder.
Report the remainder as a sequence of bits.

In practice there’s a little more to the algorithm than this, such as appending the length of the file, but the above pattern is at the heart of the algorithm.

There’s a common misconception that the polynomial Q(x) is irreducible, i.e. cannot be factored. This may or may not be the case.

CRC-32

Perhaps the most common choice of Q is

Q(x) = x³² + x²⁶ + x²³ + x²² + x¹⁶ + x¹² + x¹¹ + x¹⁰ + x⁸ + x⁷ + x⁵ + x⁴ + x³ + x² + x + 1

This polynomial is used in the cksum utility and is part of numerous standards. It’s know as CRC-32 polynomial, though there are other polynomials occasionally used in 32-bit implementations of the CRC algorithm. And it is far from irreducible as the following Mathematica code shows. The command

    Factor[x^32 + x^26 + x^23 + x^22 + x^16 + x^12 + 
           x^11 + x^10 + x^8 +  x^7 + x^5 + x^4 + 
           x^3 + x^2 + x + 1, Modulus -> 2]

shows that Q can be factored as

(1 + x)⁵ (1 + x + x³ + x⁴ + x⁶) (1 + x + x² + x⁵ + x⁶)
(1 + x + x⁴ + x⁶ + x⁷) (1 + x + x⁴ + x⁵ + x⁶ + x⁷ + x⁸)

(Mathematica displays polynomials in increasing order of terms.)

Note that the factorization is valid when done over the field with 2 elements, GF(2). Whether a polynomial can be factored, and what the factors are, depends on what field you do your arithmetic in. The polynomial Q(x) above is irreducible as a polynomial with real coefficients. It can be factored working mod 3, for example, but it factors differently mod 3 than it factors mod 2. Here’s the factorization mod 3:

(1 + 2 x² + 2 x³ + x⁴ + x⁵) (2 + x + 2 x² + x³ + 2 x⁴ + x⁶ + x⁷)
(2 + x + x³ + 2 x⁷ + x⁸ + x⁹ + x¹⁰ + 2 x¹² + x¹³ + x¹⁵ + 2 x¹⁶ + x¹⁷ + x¹⁸ + x¹⁹ + x²⁰)

CRC-64

The polynomial

Q(x) = x⁶⁴ + x⁴ + x³ + x + 1

is known as CRC-64, and is part of several standards, including ISO 3309. This polynomial is irreducible mod 2 as the following Mathematica code confirms.

    IrreduciblePolynomialQ[x^64 + x^4 + x^3 + x + 1, Modulus -> 2]

The CRC algorithm uses this polynomial mod 2, but out of curiosity I checked whether it is irreducible in other contexts. The following code tests whether the polynomial is irreducible modulo the first 100 primes.

    Table[IrreduciblePolynomialQ[x^64 + x^4 + x^3 + x + 1, 
        Modulus -> Prime[n]], {n, 1, 100}]

It is irreducible mod p for p = 2, 233, or 383, but not for any other primes up to 541. It’s also irreducible over the real numbers.

Since Q is irreducible mod 2, the check sum essentially views its input P(x) as a member of the finite field GF(2⁶⁴).

Curvature at Cairo

Posted on 27 August 2023 by John

I was flipping through Gravitation [1] this weekend and was curious about an illustration on page 309. This post reproduces that graph.

The graph is centered at Cairo, Egypt and includes triangles whose side lengths are the distances between cities. The triangles are calculated using only distances, not by measuring angles per se.

The geometry of each triangle is Euclidean: giving the three edge lengths fixes all the features of the figure, including the indicated angle. … The triangles that belong to a given vertex [i.e. Cairo], laid out on a flat surface, fail to meet.

I will reproduce the plot in Python because I’m more familiar with making plots there. But I’ll get the geographic data out of Mathematica, because I know how to do that there.

Geographic information from Mathematica

I found the distances between the various cities using the GeoDistance function in Mathematica. The arguments to GeoDistance are “entities” which are a bit opaque. When using Mathematica interactively, you can use ctrl + = to enter the name of an entity. There’s some guesswork, e.g. whether I meant New York City or the state of New York when I entered “New York”, but Mathematica guessed correctly. The following code lists the city entities explicitly.

    cities = {
        Entity["City", {"Cairo", "Cairo", "Egypt"}], 
        Entity["City", {"Delhi", "Delhi", "India"}], 
        Entity["City", {"Moscow", "Moscow", "Russia"}], 
        Entity["City", {"Brussels", "Brussels", "Belgium"}], 
        Entity["City", {"Reykjavik", "Hofudhborgarsvaedhi", "Iceland"}], 
        Entity["City", {"NewYork", "NewYork", "UnitedStates"}], 
        Entity["City", {"CapeTown", "WesternCape", "SouthAfrica"}], 
        Entity["City", {"PortLouis", "PortLouis", "Mauritius"}] }

Most of these are predictable, but I would not have guessed the code for Reykjavik or Cape Town. I found these by using the command InputForm and entering the entities as above.

I found the distance from Cairo to each of the other cities with

    Table[GeoDistance[cities[[1]], cities[[i]]], {i, 2, 8}]

and the distances from the cities to their neighbors with

    Table[GeoDistance[cities[[i]], cities[[i + 1]]], {i, 2, 7}]
    GeoDistance[cities[[8]], cities[[2]]]

Drawing the plot

Now that we’ve got the data, how do we draw the plot?

Let’s put Cairo at the origin. First we draw a line from Cairo to Delhi. We might as well put Delhi on the x-axis to make things simple.

Next we need to plot Moscow. We know the distance R₁ from Cairo to Moscow, and the distance R₂ from Delhi to Moscow. So imagine drawing a circle of radius R₁ centered at Cairo and a circle of radius R₁ centered at Delhi. Moscow is located where the two circles intersect. The previous post shows how to find the intersection of circles.

The two circles intersect in at two points, so which do we choose? We choose the intersection point that preserves the orientation of the original graph (and the globe). As we go through the cities in counterclockwise order, the cross product of the previous line to the next line should have positive z component.

This shows that the original graph was not to scale, though the gap between triangles was approximately to scale. In hindsight this should have been obvious: Brussels and Reykjavik are much closer to each other than Capetown and New York are.

The gap

Why the gap? Because the earth is curved at Cairo (and everywhere else). If the earth were flat, the triangles would fit together without any gaps.

There’s no gap when you take spherical triangles on the globe. But even though the triangles preserve length when projected to the plane, they cannot preserve angles too. The sum of the angles in a spherical triangle adds up to more than 180°, and the amount by which the sum exceeds 180° is proportional to the size of the spherical triangle. Since the angles of triangles in the plane do add up to 180°, each flat triangle fails to capture a bit of the corresponding spherical triangles, and the failures add up to the gap we see in the image.

[1] Gravitation by Misner, Thorne, and Wheeler. 1973.

Creating a Traveling Salesman Tour of Texas with Mathematica

Posted on 10 August 2023 by John

A Traveling Salesman tour visits a list of destinations using the shortest path. There’s an obvious way to find the shortest path connecting N points: try all N! paths and see which one is shortest. Unfortunately, that might take a while.

Texas has 254 counties, and so calculating a tour of Texas counties by brute force would examine 254! paths, over 10⁵⁰⁰ paths. In theory, large Traveling Salesman problems are unsolvable. In practice they can often be solved quickly. As is often the case, the key is to give yourself just a little slack and look for solutions that are close to optimal.

I’ve used the example of a Traveling Salesman tour of Texas before because it makes a nice visual. People asked me for the code that made the image, but I didn’t save the code and didn’t remember offhand how to re-create it. So here’s the code for future reference.

Incidentally, computing the tour itself took only a second or two. Creating the visualization took several seconds.

    texas = Entity["AdministrativeDivision", "Texas"]; 
    counties = texas["Subdivisions"];
    tour = FindShortestTour[texas["Subdivisions"]];
    GeoGraphics[{Thick, Red, GeoPath[counties[[tour[[2]]]]]}]

Here counties is a list of objects representing Texas counties, sorted by alphabetical order, from Anderson County to Zavala County.

The tour object is a pair of a distance and a list of integers. The distance, 6780.74 nautical miles, is the length of the tour. The integers are the indexes of the counties in the tour.

{6780.74 nmi, {1, 107, 234, …, 201, 37, 1}}

The tour starts with the first county, Anderson County. It’s got to start somewhere, and I expect it always starts with the first item in the list. Next it goes to the 107th county, Henderson County, and so on. Because FindShortestTour returns a closed tour, the tour ends where it started, in Anderson County.

Related posts: Traveling Salesman tours of Africa, Americas, Eurasia and Oceania.

Cayley graphs in Mathematica

Posted on 30 June 2023 by John

The previous post mentioned the group S₄, the group of all permutations of a set with four elements. This post will show a way to visualize this group.

The Mathematica command

    CayleyGraph[
        SymmetricGroup[4], 
        VertexLabels -> Placed["Name", Center],
        VertexSize -> 0.4]

generates the graph below.

Cayley graph of alternating group S4

This is an interesting image, but what does it mean?

The elements of S₄ are represented by the circled numbers. The numbers correspond to the permutations of four elements, listed in lexicographical order. If you label the four elements a, b, c, and d then the permutations are listed in alphabetical order. Permutation 1 is [1, 2, 3, 4] to itself and Permutation 24 is its reverse [4, 3, 2, 1].

In the Mathematica application, mousing over a number shows which permutation it represents, though the static image above doesn’t have this feature.

The blue arrows represent the permutation that swaps the first two elements. So the blue arrow between node 1 and node 7 says that swapping the first two elements of Permutation 1 gives you Permutation 7, which is [2, 1, 3, 4]. The blue arrow going back from 7 to 1 says that the same swapping operation applied to Permutation 7 returns you to Permutation 1.

All the blue arrows come in pairs because swapping is its own inverse.

The green arrows represent a rotation. For example, the green arrow from 1 to 10 says that rotation turns [1, 2, 3, 4] into [2, 3, 4, 1]. The rotation operation is not its own inverse, so the arrows only go in one direction. But every green arrow is part of a diamond because applying the rotation operation four times sends you back where you started.

You can get from any permutation to any other permutation by repeatedly either swapping the first two elements or applying a rotation. In group theoretical terminology, these two permutations generate the group S₄.

Drag equation exponent variation

Posted on 26 June 2023 by John

The motion of a falling body of mass m is given by

$m \frac{dv}{dt} = mg - kv^r$

where the term −kv^r accounts for drag due to air resistance. One can derive r = 2 under simple physical assumptions, but if I remember correctly other values of r may be more realistic in certain circumstances. I don’t know much about the physics here; if you know about the use of other values of r, please let me know by leaving a comment.

Terminal velocity

When r = 1 or r = 2 the differential equation above can be solved in terms of elementary functions, but otherwise it cannot. Nevertheless one can show that for all values of r the object reaches a terminal velocity, and calculate that velocity without explicitly solving the differential equation. William Waterhouse demonstrated this in a one-page article [1]. He rewrites the equation to look at time as a function of velocity rather than velocity as a function of time

$\frac{dt}{dv} = \frac{1}{g} \frac{1}{1 - (k/mg)v^r}$

and concludes

$t = \frac{1}{g} \int_0^v \frac{dv}{1 - (k/mg)v^r} + t_0$

He notes that the integral diverges as v approaches

$\left(\frac{mg}{k}\right)^{1/r}$

and so that is the terminal velocity, i.e. it takes an infinite amount of time to achieve this velocity. Waterhouse recommends this derivation as “a good example of deriving information about a problem without knowing an explicit solution.”

I would add that such an approach is the norm, not an exception. A closed-form solution to a differential equation in nice when you can get it, but usually not possible. And even when you can find a closed-form solution, you may be able to achieve your goal more directly by not using it.

Hypergeometric solution

I suspected the differential equation could be solved for general values of r using special functions, and that is the case. Mathematica was able to evaluate the integral for t as a function of v in terms of a hypergeometric function.

$v \, _2F_1\left(1,\frac{1}{r};1+\frac{1}{r};\frac{g k v^r}{m}\right)$

When I asked Mathematica to solve the differential equation directly, it said that the solution is the inverse function of the function above. Apparently Waterhouse and Mathematica agree that it is easier to think of t as a function of v rather than the original formulation.

The notation ₂F₁ indicates there are two upper parameters and one lower parameter. In our application, the upper parameters are 1 and 1/r, the lower parameter is 1 + 1/r, and the function is evaluated at gkv^r/m. You can find a brief introduction to hypergeometric functions here. A hypergeometric function ₂F₁ has a singularity at 1, and so we could derive the terminal velocity from the explicit solution.

Mathematica implementation

Let c = gk/m. Then we can express velocity as a function of time in Mathematica by

    f[r_, c_, v_] := InverseFunction
    [
        #1 Hypergeometric2F1[1, 1/r, 1 + 1/r, c #1^r] &
    ][t]

and use this to make plots of the velocity for various values of c and r.

The following sets c = 2 and varies r over 1, 1.1, 1.2, … 2.

    Plot[Table[f[2, d/10, t], {d, 10, 20}], {t, 0, 4}, PlotRange -> All]

Here’s the output.

The terminal velocity decreases as r increases. The opposite is true for c < 1.

[1] William C. Waterhouse. A Fact about Falling Bodies. Mathematics Magazine, Vol. 44, No. 1 (Jan., 1971), pp. 33–34. The article straddles two pages, but takes up less than half of each page.

Double duals of polyhedra

Posted on 6 April 2023 by John

The previous post mentioned the dual of a tetrahedron is another tetrahedron. The dual of a cube is an octahedron and the dual of an octahedron is a cube. And the dual of a dodecahedron is an icosahedron, and the dual of an icosahedron is a dodecahedron.

So if you take the dual of a regular solid twice, you get back another regular solid of the same type. But you do not get the same regular solid back. To see this, let’s look at how you form the dual of a polyhedron.

To find the dual of a polyhedron, you make a new polyhedron whose vertices are the centroids of the faces of the original polyhedron. The vertices of a regular polyhedron all lie on the surface of a sphere. A face of the polygon lies inside the sphere, so the centroids of the solid are all inside the sphere. Since the centroids form the vertices of a regular solid, they also all lie on a sphere, but a smaller sphere inside the first sphere.

Size of duals

If you take the double dual of a regular solid, the dual of the dual, then you get a smaller solid of the same type. How much smaller? We can get Mathematica to calculate this for us.

The function DualPolyhedron does what the name implies. Let’s apply this to the default tetrahedron. When we enter

    DualPolyhedron[ Tetrahedron[] ]

we get back

    Tetrahedron[{2 Pi/3}, Pi}, 1/3]

The default tetrahedron, returned by calling Tetrahedron[], has vertices

(1, 0, 0)
(0, 0, 1)
(1, 0, 1)
(1, 1, 1)

Each edge has length 1. Its dual has been rotated by 2π/3 about the z axis and by π about the y axis, and has edges of length 1/3. Taking the dual of a tetrahedron rotates the tetrahedron and shrinks it by a factor of 3 (in each dimension).

Size of double duals

If we take the double dual of a tetrahedron, we get a tetrahedron nine times smaller. Let’s see what happens to the rest of the regular solids.

    DualPolyhedron[DualPolyhedron[#]] & /@ 
        {Tetrahedron[], Octahedron[], Cube[],
         Dodecahedron[], Icosahedron[]}

This tells us again what we found above for the double dual of a tetrahedron. It also tells us that taking the double dual of an octahedron or cube results in an octahedron or cube 3 times smaller. And it tells us that taking the double dual of a dodecahedron or icosahedron reduces its size by a factor of

(5 + 2√ 5)/15 = 0.6314…

Taking the double dual shrinks a regular solid, and the amount of shrinkage goes down as the number of faces goes up.

Orientation of double duals

The process of taking the double dual may also rotate the solid. The second dual of the standard tetrahedron, for example, is not just 9 times smaller. It’s also been rotated. We can see this by calling PolyhedronCoordinates. The code

    PolyhedronCoordinates[ 
        DualPolyhedron[ DualPolyhedron[ Tetrahedron[] ] ] ]

shows that the new vertices are

(7/9, 2/9, 6/9)
(7/9, 2/9, 7/9)
(7/9, 3/9, 7/9)
(6/9, 2/9, 7/9)

So the second dual of the tetrahedron is not just 9 times smaller.

However, taking the double dual of the rest of the regular solids does not change the orientation. You could verify this for the cube and dodecahedron using Mathematica. It then follows that taking the double dual of octahedra and icosahedra does not change the orientation.

Sixth dual

If you take the dual of a tetrahedron 6 times you do return to the original orientation, but not if you take the dual fewer times. You can verify this by running the following:

    NestList[DualPolyhedron, Tetrahedron[], 6]