An almost perfect Babylonian grid
If there were a perfect Babylonian grid, it would
be a tesselation of the interior surface of a celestial sphere for which every patch
was a spherical triangle and every vertex was
surrounded by six identical patches. Not possible, in part because of
Euler's formula,
F - E + V = 2 , (1)
which holds for any selection of polygonal patches. For present purposes a spherical triangle will be taken to be
a tile that is bounded by three geodesic arcs (great circle arcs) of finite length
( > 0 and < pR ) that intersect
at three vertices with interior angles greater than 0 and less than 180°. When all the patches are such triangles
(call that a triangulated graph) each vertex has an associated number of (three or more) connected neighbors,
the number of vertices is V=N3+N4+N5+. . ., the number
of faces is F=(3N3+4N4+. . .)/3, and the number of edges is
E=(3N3+4N4+. . .)/2. The substitution of those in eq(1) gives
3N3 + 2N4 + N5 = 12 + N7 + 2N8
+ 3N9 + . . .(2)
So N6 is not limited by eq(2), but some of N3,
N4 and N5 have to be present. There are nineteen ways to satisfy eq(2) when
Nk=0 for k>6, but there may be some that can't be realized as inscribed
polyhedrons. The most familiar of the polyhedrons are N3=4 (tetrahedrons),
N4=6 (octahedrons) and N5=12 (icosahedrons). The
regular octahedron is
easiest to begin with because it is relatively easy to visualize an octant of the sphere and, at the same
time, an inscribed equilateral triangle that lies beneath it and another equilateral triangle that is
related to it in an entirely different way.
Before anything more is said about spherical triangles, this may be a good place to
mention
Babylonian trigonometry, where
sin(q°) = sin(pq/180) and cos(q°) = cos(pq/180) ,(3)
and the presence or absence of a '°' serves to indicate degrees or radians, where necessary.
In such terms the latitude at the north or south pole is 90° (north or south).
Either of the poles can be taken to be a point
where a triangle touches the octant above it, and the longitudes of the other two points
where the triangle touches the octant on the equator can be taken to be f=0 and f=90°.
The edges of the octant are on great circles whose centers are the center of the sphere.
They intersect one another at right angles, so the octant can be called a 90-90-90° spherical
triangle. The triangle beneath it is a 60-60-60° plane triangle, called equilateral, but equiangular
might be better in this context. For the regular tetrahedron and icosahedron the angles are
120-120-120° and 72-72-72°, respectively, and a 75-60-75° spherical triangle will appear presently.
For an arc that extends for q radians on a circle
of radius R, the length of the arc is Rq, and the length of the chord beneath it is
L(q,R) = 2Rsin(q/2) ,
(4)
in modern trigonometric notation. It's pretty certain the Babylonians didn't have a formula
like that, but they were quite familiar with arcs and chords and the Pythagorian theorem, so
it is quite likely they knew about the result, RÖ_2, for the
arcs at the edges of an octant.
Only a few of the dozens of mappings of a spherical earth onto a plane map will be
mentioned here. The first is the projection where any point in an inscribed
polygon (often a triangle) and the associated point on the sphere are at intersections of the ray
through the center of the sphere. Straight lines on the polygon are mapped to arcs of great
circles, and such maps of plane triangles are spherical triangles. Among faces of
a triangulated, inscribed polyhedron, congruent triangles define spherical patches
that have equal areas. In general, congruent triangles within a larger inscribed polygon
are not projected to congruent spherical triangles; the ones whose centroids are nearer to
the centroid of the polygon are projected to larger areas on the sphere.
A different mapping of an octant to a triangle is simply to use latitude and longitude
to define points on the sphere, and then to use triangular regions of latitude vs. longitude
to refer to parts of the sphere that are not necessarily within spherical triangles. So let
l be latitude and f be longitude. For triangular maps
it is convenient now to introduce a complimentary measure of longitude where y = 90°- f.
Then the local picture of any octant that has two vertices on the equator can be represented
as one or the other of
n=9 w=9 * * * s=0 * * * e=9
* * * * * * * * * * *
* * * * * * * * * * *
* * * * e * * * * * 0
0 * * * w o = * * * * =
= * * * * = r 0 * * * w
e * * * * * 0 * * * *
* * * * * * * * * * *
* * * * * * * * * * *
w=9 * * * n=0 * * * e=9 s=9
A few notes now, in no particular order:
- The scale of latitudes is uniform, and dl = D° from row to row, where D is a divisor of 90°.
D is 10° here. so these can be called pieces of a 10-degree computational grid. The number of D°
intervals is M = 90°/D (=9 here).
- The scale of longitudes is also uniform, either way, but the coordinate axes are oblique,
with the l-axes aligned at 60° to the (rightward) f-axis and at 60° to the
(leftward) y-axis.
- On the equator e + w = M, at the poles e = w = 0, and on
the entire triangle e + w + n = M or e + w + s = M. Another way to put that is to
say that the fixed span of longitudes (90°) is divided in M subintervals at the equator, and the
number of subintervals decreases uniformly to 0 at n = M or s = M.
This kind of graph is sometimes
called a ternary plot,
and the coordinates are called
triangular coordinates.
- The number of interior patches is M2 (sum the first M
odd integers). A square or rhombic grid (e.g. l vs f or l vs y)
whose interior patches were each bisected by a single diagonal would have 2M2 triangles.
- The pole at (X,Y,Z)=(0,0,R) could equally well be located at (0,R,0) or (R,0,0), and those
arrangements of 'latitudes' vs 'longitudes' are also shown by the diagrams. The upright ones, with
a mapmaker's definition of spherical coordinates, will be used to compute positions and
define tangent planes on the sphere. Thus,
Z = Rsin(l(n.)) , X = Rcos(l)cos(f(l,e.))
and Y = Rcos(l)cos(y(l,w.)), (5)
where f + y = 90° and -M -< n -< M.
(Note, s = -n in this view of things.)
- Points in the triangles can be identified by the coordinate pairs (e,w), accompanied by
whichever of n = M - e - w or s = M - e - w is appropriate. With the notation
'r.r' for integer- and fractional-parts of real numbers, the gridpoints are at (e.,w.,n.)
or (e.,w.,s.). The integer parts (e.,w.) specify a pair of interior triangles,
and the fractional parts can be used to determine whether the point is in the north-facing or
south-facing one and to define triangular coordinates, with .e + .w + .n = 1
or .e + .w + .s = 1, within it.
The big picture is a spherical earth with the usual North and South poles and with four
'equatorial poles', the Atlantic, East, Pacific and West poles. A typical piece of a
computational grid (AE, say) looks more or less like this:
N
* * '15e°' is longitude (east)
0 * w '15w°' is longitude (west)
= * * = '15n°' is latitude (north)
e * * * 0
* * * * * *
A n = s = 0 E p/2D (radians) = 90°/D (degrees) and D = 15
* * * * * *
e * * * 0
= * * = '15s°' is latitude (south)
0 * w
* *
S
One of three similar views of the entire computational grid (for D = 15) looks like this:
N N N N
* * * * * * * *
* * * * * * * * * * * *
* * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * *
o * * * * o=o * * * * o=o * * * * o=o * * * * o
A * * + * * E * * + * * P * * + * * W * * + * * A
o * * * * o=o * * * * o=o * * * * o=o * * * * o
* * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * *
* * * * * * * * * * * *
* * * * * * * *
S S S S
In this the vertices marked 'N', 'S', 'A', 'E', 'P' & 'W' are poles, and there is only one value of any
scalar function there. Vertices marked 'o' are on edges that have been cut to make it possible
to lay the grid out on the page, and they appear in pairs for which
values of scalar functions are the same.
Similar pictures are there for the other two diagrams, and the results are almost perfect grids,
in which every computational vertex is the center of a regular hexagon. The challenge is to find an
unbiased way to treat the phantom patches that are indicated by 'o=o' and the neighbors
of all of the vertices on edges that have been cut.
First let it be noted that it is not the poles that are a problem; values of whatever needs to be
computed there can be assigned without conflicts. It is the interior
patches along the edges of the octants that require an unbiased treatment of data from
neighboring patches outside of the octant they inhabit.
A diagram of it all for a quadrant (NEP & PES (say) for D = 15) looks like this
NEP
o o o o o o o o o o
o n o o w * * * * * e o
o * * o o * * PES * * o
NAE o * * * o NPW o * * * * * o
o * * * * o o * * * * o
o * * * * * o EAS o * * * o WPS
o * * NEP * * o o * * o
o w * * * * * e o o s o
o o o o o o o o o o
PES
Note that within any of the octants 'e' and 'w' merely mean east and west of one another, regardless
of which of the equatorial poles they are associated with. Likewise, 'n' and 's' are merely north
and south of one another, but they are always associated with the polar poles.
The sets of eight vertices marked oooooooo have boundary values, which are
to be assigned as or computed from data in the exterior faces.
There are at least two ways to proceed, and the first is by interpolation of external
data, as indicated by this diagram for the we-edge,
* * NEP * * * * NEP * * o o o NEP o o o
w * * * * * e => w * * * * * e and w * * * * * e
* * PES * * o o o PES o o o * * PES * *
At the right, below and above the we-edge, values of scalar data from shorter edges at the left are
interpolated to produce data for an edge that is two segments longer. The process, call it
stitching the patches, consists of twenty-four interpolations of data on edges
that are neighbors of the twelve edges of the octahedron.
Another way to look at neighbors near the poles and on edges is to dismantle the
octahedron in an entirely different way that begins with a cut around the 'equator', relative to
any of the poles. The result is two pyramids (hemispheres) that can be
cut from equator to pole in several different ways. Given a face that is
N's context, for example, two cuts can be made as shown below for the face, NEP (in
a previous view, the the upper part of the second diamond from the left):
*
* *
A * * * W A W
* * * *
* * * * * *
+ - + -
A * * + N - * * W <=> * N *
* * * * * *
* * * * * * * *
* * * * * *
E * * * P E P
* *
*
The diagram at the left shows the patch from ANW that points to the one beneath N in
NEP, and the diagram to the right is a polar projection that eliminates the phantom
patches that don't point to anything in NEP. The important fact is that there are no
interpolations here, so this
is definitely plan-A, with a plan-B to be held in reserve, just in case.
After the icosahedron, the next useful polyhedron is an irregular one, with
hexagons at the poles instead of pentagons, and six pairs of pentagons wrapped around the equator.
For this grid D is a divisor of 60°, M is 60°/D, and the twenty-four faces can be laid out as
N N N N N N 90°
* * * * * * * * * * * *
* * * * * * * * * * * *
* o=o o=o o=o o=o o=o *
1 + 2 + 3 + 4 + 5 + 6 + 1 30°
* * * * * * * * * * * * *
* * * * * * * * * * * * *
* * * * * * * * * * * * *
3' + 4' + 5' + 6' + 1' + 2' + 3' -30°
* o=o o=o o=o o=o o=o *
* * * * * * * * * * * *
* * * * * * * * * * * *
S S S S S S -90°
As intended here, horizontal lines are lines of constant latitude, and rays that connect N
or S to a '+' or a numbered vertex are lines of constant longitude. This grid is
topologically the same as a
hexagonal antiprism.
A Digression:
Also
see ref 1 and
search for 'hexagonal antiprism'. The polyhedron there is an inflated hexagonal
antiprism. If the twenty-four faces in the diagram were meant to represent congruent, inscribed,
isosceles triangles they would be narrower than indicated above, and the levels marked ±30°
would be the latitudes of midpoints (marked '+') of pairs of perpendicular great circle arcs.
On the sphere the tiles would be congruent, 75-60-75° spherical triangles.
|
The phantom patches work quite differently here, and there is only one (marked o=o)
that is directed toward central patches at the numbered vertices. On the sphere,
the vertical arcs from N or S bisect neighboring arcs, and the use of equal values of scalars
for o=o is an unbiased way to get around Euler's formula. The polar views
are also different, as shown below.
5 * * * 4 = 4 * * * 6' 5 * * * 4 * * * 5'
* * * * * * * * * * * * *
* * * * * * * * * * * * *
* * * o=o * * * * * * - *
6 * * * N * * * 3 - * * 5' and 6 * * * N * * * 3 + * * 4'
* * * * + * * * * o=o * *
* * * * * * * * * * * * *
* * * * * * * * * * * * *
1 * * * 2 * * * 4' 1 * * * 2 = 2 * * * 3'
This time, there are no phantom patches at the poles, and
there are two ways to locate just one phantom patch at each of the twelve singularities
that aren't poles in this grid. The situation at the vertex marked '3' in the above
diagram can also be displayed as
4 * * * 5'
3 * *
+ - * * * *
* * * * * - *
use * * twice in N * * 3 ± #
4'* * * 5' * * * + *
* * * *
* *
2 * * * 4'
where '±' and '#' stand for the other vertex on the same external patch.
That and the phantom patches in the first diagram of this grid indicate unbiased
ways to surround all of the real patches without any interpolations. That is Plan A
for this grid, and there is also a Plan B, based upon interpolation on edges,
just in case.
The first thing to notice about this grid is that the number of vertices in
every central row where
|l| <- 30°
is 360/D, and their spacing is uniform (all three ways). Thus the positions of vertices are
more conveniently described in terms of ordinary spherical coordinates
in the tropical latitudes. To keep the triangular symmetry,
just use the rhombic coordinates for either of l vs f or l vs y.
The triangles are still there, but with a uniform spacing of vertices within. If M is even the
south-polar cap can be rotated by 30°, and diagrams that show vertices and some of the
special cases for M = 4 are:
P P P 90°
* * * * * *
* * * * * * * V
* o o * * * * * *
V * * V * * V * * * * 30° * * *
* * + - * * * * * * * -
* * * * * * * * and P * * V ±
* * * * - + * * * * * +
* * * * V * * V * * V 30° * * *
* * * * * o o * *
* * * * * * * V
* * * * * *
P P P -90°
Notes about the second grid, again in no particular order:
- Because vertices are treated differently in polar and tropical latitudes here,
this is an example of what are sometimes called hybrid grids. Some of them have names,
but no name other than Babylonian has been found for this one, so far. The six
polar triangles fold down flat, so if Babylonian seems like a frivolous name, this
could also be called a hexagonal hat-box.
- The earth map that is most like this one is the (area preserving)
sinusoidal projection in
which latitudes are horizontal lines and longitudes are sine-waves _
e.g. y(l) = l and x(l,f) = fcos(l)
on the map. Shapes are relatively true for the sector, |f|<30°, and the grid shown
above is just that on three of six segments, slightly deformed.
- In the tangent planes, dy = Rdl, dx = Rcos(l)df, and the element
of area is
dA = dxdy = R2cos(l)dldf =
R2d(sin(l))df.
(6)
- The direction that is 'north' or 'south' at vertices within the polar faces is that of a ray that
contains the vertex and the pole. Except for vertices on e = 0 or w = 0, the ray does not intersect
the edge above or below at another vertex, and different interpolations are required to
define y-derivatives in different tangent planes. At vertices within the central region, it is
always the same interpolation (at the midpoint between vertices). At the latitudes, ±30°,
both ways happen.
- In the event interpolations in the patches are employed, the functions
are f(e.e,w.w,t),
with specified values of f(e.,w.,t). By the construction of the grid, the interpolation formula
is exactly the same everywhere on it.
- The grid looks somewhat like a bucky dome, but in fact the only latitude-line that
is a geodesic on the sphere is the equator. The grid is not derived from a projection of the hexagonal
antiprism, and it is more like the segmented domes at Biosphere 2.
In map-making,
area-preserving maps are often preferred (but globes are better). For latitude-longitude
grids, area can't be addressed without specifying precisely what kinds of lines are to be
drawn to
connect the vertices on the sphere. If the uniformity of the oblique coordinate systems is
maintained, latitude varies linearly with longitude on the connecting lines, and the areas of
of the l-f triangles on the sphere simply are not all equal to one another.
So what remains to be seen is what almost perfect means here.
There are two different effects in play: they are (1) the separations of lines
of constant latitude
and (2) the areas of two different kinds of triangles within zonal bands. Consider first
the idea of shifting the latitude lines at '60°/D' to lines at latitudes lk
that make the zonal averages of the areas of patches (on the sphere) equal.
From eq(6) it follows that the latitudes of vertices are
shifted, proceeding from |l| = 90° at either pole to |l| = 30°
in twelve polar faces, according to
sin|lk| = sin|lk-1| -
(2k-1)/2M2 for
1 <- k <- M
with sin|l0| = 1 .(7a)
The corresponding shifts in the twenty-four triangular faces in the tropics are
sin(lk) = sin(lk-1) +-
1/M for
1 <- k <- M
with sin(l0) = -+1/2 .(7b)
The partial sums of the recurrence relations are:
sin|lk| = 1 - (k/M)2/2 ,
1 <- k <- M
(8a)
for |l| in [90°,30°], and
sin|lk| = |1/2 - k/M| ,
1 <- k <- M
(8b)
for |l| in [30°,0°].
With these choices, the averages of areas of patches (of two kinds) in zonal bands
(defined as within the intervals, [k,k±1]), are all equal to
(p/6)(R/M)2.
This grid defines computational vertices whose images on the sphere are more
nearly uniformly distributed than they were before. To compare it with the first grid where
|lk| = |30 - kD°| for |l| < 30° and
|lk| = |90 - kD°| for 30° < |l| < 90°,
we observe that eqs(8) define a hybrid approximation of sin(l) that is
linear for |l|<30°, quadratic for |l|>30°, and correct at
0, 30° and 90°.
The approximation, expressed in radians, is
sin|l| »
3|l|/p for |l| -< p/6
(1/2)(1 +
(3/p)2(|l| - p/6)(5p/6 - |l|))
for |l| -> p/6
(9)
The slope of the approximation is approximately 0.955 where the exact result decreases from 1 to
Ö_3/2,
and the magnitude of the second derivative is approximately 0.912 where the exact result increases
from Ö_3/2, to 1. The approximation in eq(9) is significantly better than
the (2/p)2|l|(p-|l|) that comes from a similar
treatment of the faces of the octahedron. The zonal averages can be made uniform by simply displacing
vertices according to eqs(8). Because eq(9) is not qualitatively wrong, it is an extra computation
that may not be necessary. However, if the exact l's are
incorporated, the Babylonian connection is all but severed, and the grid might better be described
as almost area-preserving.
As for the different areas within zonal bands, so what? The vertices are spaced uniformly
on latitude lines, and areas of individual patches will be treated as a detail that can easily be dealt
with if ever it becomes important. A comparison of the two grids is summarized below:
latitude
| | Babylonian
(degrees) |
area-preserving
(radians) | range |
|
| |l| | = | 30°(3-2(e+w)/M) |
sin-1(1-(e+w)2/2M2) | |l| > 30° |
| 30°|1-2(e+w)/M| | sin-1|1/2-(e+w)/M| |
|l| < 30° |
| |
More comments:
- In semi-Lagrangian algorithms for transport equations, linear interpolation,
on edges or in patches, isn't good enough ! The error is systematic, and the resolution
of steep fronts is always seriously compromised by numerical dissipation
(a.k.a. numerical viscosity or diffusion).
- The systematic error can be eliminated, and the dissipation can even be reversed in stable
algorithms that use higher order interpolations (e.g. Hermite interpolation). For patches, the
triangular coordinates help to simplify the polynomial bases
for interpolation in a triangle that is surrounded by six identical neighboring triangles.
- After the evaluation of the values for the basis, operation counts (per function) run from
3'x's , 2'+'s (linear), through 9'x's , 8'+'s (plus six cubic
elements), to 12'x's , 11'+'s (plus three quartics whose sum is the last cubic, ew(1-e-w)).
Having identical elements of whichever basis is adopted eliminates most of the unnecessary
arithmetic (for thousands of patches).
- In a grid of triangular patches that are implemented as structures, including pointers to neighbors,
the drawing of contour lines is about as easy as it ever gets.
- When selective grid bisections are introduced, their rather simple family tree is
n n
* * * *
* * * *
* * <==> w * * e
* * n * * n
* * * * * * * *
w * * * * * e * * s * *
w * * e w * * e
And finally: The Babylonians had no love for multiplication or division, except with numbers
that have a few common factors, and they were remarkably good at avoiding those two operations.
Why? Just think about the fourth and fifth grades and numbers that have sixty degrees, rather
than a mere ten digits.
A few references
Online:
hexagonal prism grid,
Euler's formula,
StAndrews-Babylonians.
examples of earth maps,
Also, from wickipedia:
chord length
triangular coordinates,
ternary plot,
barycentric coords,
hexagonal antiprism.
Numerical Weather and Climate Prediction, Thomas Tomkins Warner, Cambridge U. Press, 2011.