Vector-Algebra Algorithms to Draw the Curve of Alignment, the Great Ellipse, the Normal Section, and the Loxodrome

,


Introduction
In a geographic information system (GIS), linear and polygonal features are usually defined by only their vertices, so the programmer is assigned the task of deciding which curve to use to connect the points; that is, the task of defining the implied edge [1].If the edge is to be a straight line in "the real world", then it could be found conceptually by stretching a chalk line between the endpoints and snapping it down, leaving the trace of the curve on the ground-or, for geomatics, on the reference ellipsoid that is used by the geodetic datum for the mapping.Such a curve is called a geodesic, being "the curve that does not veer to the left nor to the right".Then, the geodesic would be projected onto the map using the forward map-projection formulas.
The formulas for calculating a geodesic curve on an ellipsoid (hereafter just called a geodesic) are very complicated [2][3][4][5][6][7][8], and there are many special cases.For example, there are four geodesics connecting any two antipodal points on the Equator: two geodesics following the Equator itself, one to the east and the other to the west, and two others following the northwards and southwards meridians connecting the points.There are other peculiar cases as well, such as the "lift-off" condition [9] for which two points on the Equator are connected by either one or three geodesics depending on how far apart they are.Since there is not always exactly one geodesic connecting two points, and given the programming challenge of implementing them, a programmer might want to consider other curves that are very nearly geodesics but simpler to implement: normal sections, great ellipses, loxodromes, and curves of alignment.
This paper presents a single, unified approach for drawing these curves that is based on a parametric representation using vector algebra.(It would be ideal to provide an algorithm for the geodesic as well; however, there does not appear to be a way to implement general geodesics with this approach.)Parametric solutions have already been found for the great ellipse [10,11], but these approaches use geodetic coordinates for their parameters, which is very helpful for navigation but complicates the formulas if the goal is only to draw the curve.The direct and indirect problems for loxodromes were provided by Meyer and Rollins [12], but their formulations were also parameterized by geodetic coordinates, resulting in similarly complicated formulas [3] and Deakin [13] report formulas for points along a curve of alignment.Testing revealed that the method below produces the same points as Deakin [13], and the method below is simpler, in the opinion of the author.
Notice that drawing these curves between two points is a separate problem from the typical "direct" and "indirect" problems of geodesy.The former takes a point, an azimuth, and a distance, and the implied forward point is to be determined.The latter takes two points, and the forward and backwards azimuths, and their separation is to be determined.One can use the indirect problem to find the forward azimuth and the separation between the points, and then use the direct problem's solution parametrically as a function of distance to enumerate the points along the curve to be drawn.
The terminology and notation for the remainder of the paper are now given.The lengths of an ellipsoid's semi-major and semi-minor axes are a and b, and they are related by the first eccentricity as b and geodetic height is h.See Appendix A for other formulas.The geocentric Cartesian coordinate axes are X, Y, and Z, with Z being the Z-axis of the reference ellipsoid.The Xand Y-axes are perpendicular to Z and to each other, and they also define the Equatorial Plane.The Prime Meridian is the section formed by the X-Z plane.(The reader will recall that a section of an ellipsoid is the curve formed by the intersection of a plane and the ellipsoid.) The general idea of the algorithms below is to produce a set of XYZ points along the straight-line segment (i.e., the chord) connecting the endpoints of the curve to be drawn.There are simple relationships between these points, which are in the interior of the ellipsoid, and their images on the ellipsoid.The relationship depends upon which curve is desired.An XYZ point, or vector, will be denoted by P = X P , Y P , Z p .A point in geodetic longitude, latitude, and height (LBH) will be denoted by P = λ P , ϕ P , h P .

Curve of Alignment
The curve of alignment can be visualized in this way (c.f.[13] (Deakin, 2009)).Consider two points of interest, A and B, on the ellipsoid, and imagine connecting them with a straight line that penetrates through the ellipsoid.Now, imagine another line segment that is perpendicular to the connecting line and that rises outwards, emerging perpendicularly (normal) to the ellipsoid.This rising line is free to slide along the connecting line and it remains normal to the ellipsoid as it moves.The curve of alignment is the locus of points at the intersection of the rising line and the near side of the ellipsoid.

Properties
The curve of alignment is undefined if A and B are antipodal.Since A and B are on an ellipsoid, their curve of alignment spans less than 180 • of longitude.
The curve of alignment from A to B is the same as the curve of alignment from B to A. Any portion less than 180 • of the Equator is a curve of alignment and is also a geodesic, as is any such portion of a meridian.A portion of any parallel other than the Equator is not a curve of alignment.(The curve of alignment connecting two points at the same latitude will pass polewards of their common parallel; see Figure 1).
The curve of alignment is not in general a planar curve.In fact, it is the intersection of a hyperbolic paraboloid and an ellipsoid [3], p. 67.

Algorithm
The algorithm used to draw curve of alignment is based on geodetic coordinate-system conversions.The algorithm depends on the idea that changing a point's geodetic height only-not geodetic longitude or latitude-moves the point in a direction normal to the ellipsoid.The connecting line segment is given by the difference of A and B in geocentric Cartesian coordinates, and we formulate this segment parametrically.The rising segment is found by converting a point on the connecting segment from XYZ into geodetic longitude, latitude, and height, and then setting the geodetic height to zero.This determines a point on the ellipsoid as the other endpoint of the rising segment such that the rising segment is perpendicular to the ellipsoid.The following six steps make this algorithm explicit.
1. Some ellipsoid is specified, such as a geodetic reference ellipsoid.This provides the length of the semi-major axis a and the value of the (first) eccentricity squared  .(See Appendix A for the formulas and parameter values.)2. Two points of interest, A and B, not antipodal and on the ellipsoid are chosen.(Note: A and B are on the ellipsoid so their geodetic heights are zero.)Their coordinates are required in geocentric Cartesian coordinates:  = (  ,   ,   ) and  = (  ,   ,   ).
If given in geodetic longitude and latitude, they can be readily converted to geocentric (see Appendix A). 3. Compute the vector v = B-A.4. The points of the straight-line segment from A to B are given in parametric form by () =  +   for 0 ≤  ≤ 1. 5. Denote the conversion from geocentric Cartesian coordinates to geodetic longitude, latitude, and height by .Then, () = (, , ℎ).This conversion depends on a and  .Define a variant of  called ̅ such that ̅ () = (, , 0), where, for some h, () = (, , ℎ), i.e., ̅ () returns the same as () but with h set to zero.
6.The curve of alignment is the curve ̅ () for 0 ≤  ≤ 1.In other words, the curve of alignment is the set of geodetic points given by the conversion of all the points p(t) from geocentric Cartesian to geodetic, setting their geodetic heights to zero.
An implementation in Mathematica™ produced the image in Figure 1, which shows the curve of alignment (red) and the geodesic (green) from (0° E, 45° N) to (165° E, 40° S) on the WGS 84 reference ellipsoid.The arc length of the geodesic is 18,669,335.84m, and the arc length of the curve of alignment is 18,671,843.56m, which is 2507.72 m (134 ppm) longer.The arc length of the geodesic was computed using [5], and the arc length of the curve of alignment was computed by recursively bisecting the curve into linear segments and summing the lengths of the segments until the length converged to submillimeter levels.In this example, the curve of alignment starts out north of the geodesic, crosses exactly once near the middle, and is south thereafter.

Algorithm
The algorithm used to draw curve of alignment is based on geodetic coordinate-system conversions.The algorithm depends on the idea that changing a point's geodetic height only-not geodetic longitude or latitude-moves the point in a direction normal to the ellipsoid.The connecting line segment is given by the difference of A and B in geocentric Cartesian coordinates, and we formulate this segment parametrically.The rising segment is found by converting a point on the connecting segment from XYZ into geodetic longitude, latitude, and height, and then setting the geodetic height to zero.This determines a point on the ellipsoid as the other endpoint of the rising segment such that the rising segment is perpendicular to the ellipsoid.The following six steps make this algorithm explicit.

1.
Some ellipsoid is specified, such as a geodetic reference ellipsoid.This provides the length of the semi-major axis a and the value of the (first) eccentricity squared ϵ 2 .(See Appendix A for the formulas and parameter values.)2.
Two points of interest, A and B, not antipodal and on the ellipsoid are chosen.(Note: A and B are on the ellipsoid so their geodetic heights are zero.)Their coordinates are required in geocentric Cartesian coordinates: If given in geodetic longitude and latitude, they can be readily converted to geocentric (see Appendix A).

3.
Compute the vector v = B-A.

4.
The points of the straight-line segment from A to B are given in parametric form by Denote the conversion from geocentric Cartesian coordinates to geodetic longitude, latitude, and height by γ.Then, γ(p) = (λ, ϕ, h).This conversion depends on a and ϵ 2 .Define a variant of γ called γ such that γ(p) = (λ, ϕ, 0), where, for some h, γ(p) = (λ, ϕ, h), i.e., γ(p) returns the same as γ(p) but with h set to zero.6.
The curve of alignment is the curve γ(p(t)) for 0 ≤ t ≤ 1.In other words, the curve of alignment is the set of geodetic points given by the conversion of all the points p(t) from geocentric Cartesian to geodetic, setting their geodetic heights to zero.
An implementation in Mathematica™ produced the image in Figure 1, which shows the curve of alignment (red) and the geodesic (green) from (0 • E, 45 • N) to (165 • E, 40 • S) on the WGS 84 reference ellipsoid.The arc length of the geodesic is 18,669,335.84m, and the arc length of the curve of alignment is 18,671,843.56m, which is 2507.72 m (134 ppm) longer.The arc length of the geodesic was computed using [5], and the arc length of the curve of alignment was computed by recursively bisecting the curve into linear segments and summing the lengths of the segments until the length converged to submillimeter levels.In this example, the curve of alignment starts out north of the geodesic, crosses exactly once near the middle, and is south thereafter.

Notes:
• The vector v is equally well defined as v = A − B, but then p(t) = B + t v.The same curve of alignment results either way.

•
Pseudocode for the curve of alignment is given below.The algorithm assumes there is a function named XYZ_to_LBH that converts XYZ coordinates to geodetic longitude, latitude, and height LBH.This code is assumed to be written in a language with built-in vector algebra operators such as Mathematica or Matlab, or with a vector-algebra library such as NumPy.

Normal Sections
Normal sections can be defined with a single point on an ellipsoid and an azimuth [14].There is a vector that is normal to the ellipsoid at this point.The sectioning plane contains the point and the normal vector, and the sectioning plane is oriented by the azimuth.The normal section is the intersection of that plane and the ellipsoid.This definition is not well suited for the purposes of this paper, so an alternative based on two points is provided.A normal section from A to B is a section formed from the intersection of a plane containing A, B, and another point V, where VA is normal to the ellipsoid at A (see Figure 2).There are many choices for V, but it is convenient and helpful to choose V to be the point on the Z-axis intercepted by the inwards-directed normal vector at A (details are given below). Notes:

•
The vector v is equally well defined as v = A − B, but then () =  +  .The same curve of alignment results either way.

•
Pseudocode for the curve of alignment is given below.The algorithm assumes there is a function named XYZ_to_LBH that converts XYZ coordinates to geodetic longitude, latitude, and height LBH.This code is assumed to be written in a language with built-in vector algebra operators such as Mathematica or Matlab, or with a vector-algebra library such as NumPy.

Normal Sections
Normal sections can be defined with a single point on an ellipsoid and an azimuth [14].There is a vector that is normal to the ellipsoid at this point.The sectioning plane contains the point and the normal vector, and the sectioning plane is oriented by the azimuth.The normal section is the intersection of that plane and the ellipsoid.This definition is not well suited for the purposes of this paper, so an alternative based on two points is provided.A normal section from A to B is a section formed from the intersection of a plane containing A, B, and another point V, where VA is normal to the ellipsoid at A (see Figure 2).There are many choices for V, but it is convenient and helpful to choose V to be the point on the Z-axis intercepted by the inwards-directed normal vector at A (details are given below).

Properties
The normal section is a planar curve.On an ellipsoid, the normal section is the entire ellipse (a circle for the Equator), but it is common practice to consider the normal section to be the shorter limb from A to B. The algorithm below produces the shorter limb, which is consistent with the goal of this paper.
A, B, and V will not define a normal section if they are collinear, but selecting B to be antipodal to A is not a problem (generally), as will now be shown.Point C in Figure 2 is the point where VA is extended to intersect the ellipsoid on the hemisphere opposite A. C is collinear with A and V, but C is not antipodal to A unless A is on the Equator or is a Pole (antipodal points are on a line through the origin).VA is normal to the ellipsoid at A, but VA does not generally contain the origin due to the ellipsoid having non-zero eccentricity.Therefore, the normal section is uniquely defined if A and B are antipodal but not on the Equator and not Poles.In that case, there is a singularity at the exact midpoint, which is at the origin and thus does not possess a unique normal vector.
The normal section from A to B is not the same as the normal section from B to A unless A and B are at the same latitude or at the same longitude.
The Equator and the meridians are normal sections.Any parallel other than the Equator is not a normal section (the normal section connecting two distinct points at the same latitude other than zero will pass polewards of their common parallel; see Figure 3).

Properties
The normal section is a planar curve.On an ellipsoid, the normal section is the entire ellipse (a circle for the Equator), but it is common practice to consider the normal section to be the shorter limb from A to B. The algorithm below produces the shorter limb, which is consistent with the goal of this paper.
A, B, and V will not define a normal section if they are collinear, but selecting B to be antipodal to A is not a problem (generally), as will now be shown.Point C in Figure 2 is the point where VA is extended to intersect the ellipsoid on the hemisphere opposite A. C is collinear with A and V, but C is not antipodal to A unless A is on the Equator or is a Pole (antipodal points are on a line through the origin).VA is normal to the ellipsoid at A, but VA does not generally contain the origin due to the ellipsoid having non-zero eccentricity.Therefore, the normal section is uniquely defined if A and B are antipodal but not on the Equator and not Poles.In that case, there is a singularity at the exact midpoint, which is at the origin and thus does not possess a unique normal vector.
The normal section from A to B is not the same as the normal section from B to A unless A and B are at the same latitude or at the same longitude.
The Equator and the meridians are normal sections.Any parallel other than the Equator is not a normal section (the normal section connecting two distinct points at the same latitude other than zero will pass polewards of their common parallel; see Figure 3).The normal section from A to B does not have the same arc length as the normal section from B to A unless  =  or  = − .

Algorithm
The algorithm used to draw a normal section is similar to that for curve of alignment, but it is a little more complex, so it warrants a figure to help explain it (see Figure 2).Some of the steps given here are expounded upon below in the notes.
1. Some ellipsoid is specified, such as a geodetic reference ellipsoid.This provides the length of the semi-major axis a and the value of the (first) eccentricity squared  .(See Appendix A for the formulas and parameter values.)2. There is a point V on the Z-axis at a distance NA from A, where  =  1 −  sin  ⁄ is the radius of curvature in the prime vertical at A. V's geocentric coordinates are (0, 0,   −   sin   ) = (0, 0, −   sin   ) and VA is normal to the ellipsoid at A. The plane AVB contains VA and so AVB is the normal sectioning plane at A.  The normal section from A to B does not have the same arc length as the normal section from B to A unless ϕ A = ϕ B or ϕ A = −ϕ B .

Algorithm
The algorithm used to draw a normal section is similar to that for curve of alignment, but it is a little more complex, so it warrants a figure to help explain it (see Figure 2).Some of the steps given here are expounded upon below in the notes.

1.
Some ellipsoid is specified, such as a geodetic reference ellipsoid.This provides the length of the semi-major axis a and the value of the (first) eccentricity squared ϵ 2 .(See Appendix A for the formulas and parameter values.)

2.
There is a point V on the Z-axis at a distance N A from A, where N A = a/ 1 − ϵ 2 sin 2 ϕ A is the radius of curvature in the prime vertical at A. V's geocentric coordinates are (0, 0, Z A − N A sin ϕ A ) = 0, 0, −N A ϵ 2 sin ϕ A and VA is normal to the ellipsoid at A. The plane AVB contains VA and so AVB is the normal sectioning plane at A.

3.
A point of the straight-line segment (chord) from A to B is given in parametric form by p(t) = A + t (B − A) for 0 ≤ t ≤ 1.

4.
Any scalar multiple of a vector from V to p(t) is in the normal sectioning plane defined by AVB.Let q denote a vector in the normal sectioning plane defined by q(t, u) = V + u(p(t) − V), where u is a positive real number.

5.
There is one value for u, say u*, such that the geodetic height of γ(q(t, u * )) is zero for u ≥ 0. The point q(t, u * ) is on the normal-section curve from A to B because q(t, u * ) is in the sectioning plane and the geodetic height of γ(q(t, u * )) is zero.Notice that this is different from the curve of alignment: there, the geodetic height is set to zero, but here, a value for u must be found that produces a geodetic height of zero.
Setting the height to zero does not work for normal sections because q is not generally normal to the ellipsoid.The points on an ellipsoid satisfy x 2 a 2 + y 2 a 2 + z 2 (1−ϵ 2 )a 2 = 1, so q(t, u * ) must also.Define n = p(t) − V. Substituting in q(t, u * )'s components and recalling that V x = V y = 0 gives: which is a quadratic in u * .After some algebra, u * can be computed from where The normal section from A to B is the curve γ(q(t, u * )) for 0 ≤ t ≤ 1.In other words, the normal section from A to B is the set of points given by the set of all the points q(t, u * ) converted from geocentric Cartesian to geodetic coordinates.Notes: • The normal-section curve is parameterized by t, the same as the curve of alignment.

•
In general, each t will have its own u*.• The point on the normal section is given by V + u * (p(t) − V).
• The vector q(0, u) is normal to the ellipsoid but no other vector q(t, u) for 0 < t ≤ 1 is also normal unless A and B have the same latitude, and then q(0, u) and q(1, u) are both normal (at A and B).

•
Pseudocode for the normal section is given below, with the same assumptions as for the curve of alignment.
A := (XA, YA, ZA) # starting point in XYZ B := (XB, YB, ZB) # ending point in XYZ t: floating point scalar in [0, 1] The two normal sections (black) and the curve of alignment (red) between (0 • E, 45 • N) and (165 • E, 40 • S) on the WGS 84 reference ellipsoid are shown in Figure 3.The normal sections have arc lengths of (N-S) 18,669,545.69m and (S-N) 18,670,163.62 m, making the N-S section 617.94 m longer than the other.The lengths of the normal sections were computed with recursive subdivision, as above.

Great Ellipses
For non-antipodal points A and B, a great ellipse from A to B is a section formed from the intersection of a plane containing A, B, and the origin; therefore, the great ellipse is available from the normal-section algorithm simply by setting V = (0,0,0).The equation for q(t, u*) simplifies to q(t, u * ) = u * p(t).The value for: The algorithm is otherwise unchanged.

1.
The great ellipse is a planar curve.

2.
The great ellipse is the entire ellipse (a circle for the Equator).Every great ellipse contains A's antipodal point by construction; however, the antipodal point cannot be used to define the great ellipse because, then, A, B, and V would be collinear.
The algorithm below produces the shorter limb between A and B, which is consistent with the goal of this paper.

3.
The great ellipse from A to B is the same as the great ellipse from B to A.

4.
The Equator and the meridians are great ellipses.Any parallel other than the Equator is not a great ellipse (the great ellipse connecting two distinct points at the same latitude other than zero will pass poleward of their common parallel; see Figure 4).
sections have arc lengths of (N-S) 18,669,545.69m and (S-N) 18,670,163.62 m, making the N-S section 617.94 m longer than the other.The lengths of the normal sections were computed with recursive subdivision, as above.

Great Ellipses
For non-antipodal points A and B, a great ellipse from A to B is a section formed from the intersection of a plane containing A, B, and the origin; therefore, the great ellipse is available from the normal-section algorithm simply by setting V = (0,0,0).The equation for q(t, u*) simplifies to (,  * ) =  * ().The value for: The algorithm is otherwise unchanged.Properties 1.The great ellipse is a planar curve.2. The great ellipse is the entire ellipse (a circle for the Equator).Every great ellipse contains A's antipodal point by construction; however, the antipodal point cannot be used to define the great ellipse because, then, A, B, and V would be collinear.The algorithm below produces the shorter limb between A and B, which is consistent with the goal of this paper.3. The great ellipse from A to B is the same as the great ellipse from B to A. 4. The Equator and the meridians are great ellipses.Any parallel other than the Equator is not a great ellipse (the great ellipse connecting two distinct points at the same latitude other than zero will pass poleward of their common parallel; see Figure 4).
Figure 4 shows the great ellipse, the curve of alignment, and the normal sections between (0° E, 45° N) and (165° E, 40° S) on the WGS 84 reference ellipsoid.The arc length of the great ellipse is 18,669,406.160m, as computed by recursive subdivision.For this example, the curve of alignment starts out north of the great ellipse, then crosses once near the middle, and remains south thereafter, which is also true for the geodesic.

Loxodromes
The loxodrome is a curve of constant azimuth [12], so all the parallels of constant latitude (i.e., "parallels") are loxodromes.A loxodrome is drawn as a straight line on a Mercator-projection map, which permits a method of drawing loxodromes in geodetic coordinates using this paper's approach.Figure 4 shows the great ellipse, the curve of alignment, and the normal sections between (0 • E, 45 • N) and (165 • E, 40 • S) on the WGS 84 reference ellipsoid.The arc length of the great ellipse is 18,669,406.160m, as computed by recursive subdivision.For this example, the curve of alignment starts out north of the great ellipse, then crosses once near the middle, and remains south thereafter, which is also true for the geodesic.

Loxodromes
The loxodrome is a curve of constant azimuth [12], so all the parallels of constant latitude (i.e., "parallels") are loxodromes.A loxodrome is drawn as a straight line on a Mercator-projection map, which permits a method of drawing loxodromes in geodetic coordinates using this paper's approach.

1.
The parallels and meridians are loxodromes.

2.
A loxodrome is a spiral if its starting azimuth is not a cardinal direction.In that case, loxodromes spiral infinitely many times around both Poles but do not reach them-yet, they have finite length.For example, for the WGS 84 ellipsoid and azimuth 45 • , the total length works out to be 28,289,831.17m [12].

3.
The loxodrome from A to B is the same as the loxodrome from B to A.

4.
The loxodrome is not a planar curve in general.

5.
The loxodrome is defined if A and B are antipodal.
The algorithm used to draw loxodromes given below is based on the loxodromes being straight lines when drawn using a Mercator projection.See [15] or [3] for traditional implementations of the Mercator projection.However, Rollins and Meyer [16] provide a more succinct implementation, which is used here.Set the central meridian of the projection λ 0 to be the average of λ A and λ b , λ, which avoids issues that arise if the loxodrome crosses 180 • .However, the arithmetic average of two angles can fail, so let s = sin(λ A ) + sin(λ B ) and c = cos(λ A ) + cos(λ B ).Then, λ 0 = λ atan2(c, s), and where atan2 is a 2-argument arctangent function, atanh is the inverse hyperbolic tangent, x is the Mercator easting, y is the Mercator northing, and q is isometric latitude.The Mercator forward mapping formulas are applied to A and B, producing Mercator grid coordinates (x, y), and the vector v is given by their difference, as before.Points p(t) in the Mercator grid are computed as a multiple of v and then converted back to geodetic coordinates.The inverse formulas are from [16], as provided in Appendix A. The following six steps make this algorithm explicit.

1.
Some ellipsoid is specified, such as a geodetic reference ellipsoid.This provides the length of the semi-major axis a and the value of the (first) eccentricity squared ϵ 2 .(See Appendix A for the formulas and parameter values.)2.
Two points of interest on the ellipsoid are chosen.Their coordinates are required in geodetic coordinates; e.g., A = (λ A , ϕ A ) and B = (λ B , ϕ B ).If given in geocentric Cartesian coordinates, they can be readily converted to geodetic.

3.
Using the forward mapping equations (see above), convert A and B into Mercator coordinates a = (x a , y a ) and b = (x b , y b ). 4.
Compute the vector v = b − a.

5.
The points of the straight-line segment from a to b are given in parametric form by p(t) = a + t v for 0 ≤ t ≤ 1.

6.
Denote the inverse Mercator mapping equations (see Appendix A) by γ M .Then, The loxodrome is curve γ M (p(t)) for 0 ≤ t ≤ 1.In other words, the loxodrome is the set of points given by the conversion of all the points p(t) from Mercator grid coordinates to geodetic coordinates.

•
The vector v is equally well defined as v = a − b, but then () =  +  .The same loxodrome results either way.

•
The inverse formula for latitude requires finding a fixed point.This can be cast as a root-finding problem by finding the latitude such that the difference of the input latitude minus the output latitude equals zero.

•
Pseudocode for the loxodrome is given below.The reference-ellipsoid object RE is assumed to have a method ecc1 that returns the first eccentricity of the reference ellipsoid; otherwise, the same assumptions as above pertain.

Summary
Four new algorithms are presented that draw loxodromes, great ellipses, curves of alignment, and normal sections.The algorithms are parameterized so that a simple implementation based on vector algebra is available, given a supporting geometric geodesy library to calculate the coordinate system.There are many, more or less equivalent, choices for the conversion from geocentric Cartesian coordinates to geodetic longitude, latitude, and height [17,18].However, for long lines, the chord between the beginning and ending points can come near to the geocenter, which can cause some algorithms to fail.For such situations, a robust algorithm such as [19] is recommended.Also, the algorithms herein are language neutral and do not depend on external packages.However, it should be noted that efficient algorithms and packages exist for geodesics [20] when the actual geodesic is required.
The examples herein were computed using the WGS84 reference ellipsoid.Ref. [21] showed that the minute difference in the flattening of these ellipsoids cannot result in a difference in coordinates at the millimeter level.So, these examples are properly illustrative for GRS80, as well.Notes: • The vector v is equally well defined as v = a − b, but then p(t) = b + t v.The same loxodrome results either way.• The inverse formula for latitude requires finding a fixed point.This can be cast as a root-finding problem by finding the latitude such that the difference of the input latitude minus the output latitude equals zero.

•
Pseudocode for the loxodrome is given below.The reference-ellipsoid object RE is assumed to have a method ecc1 that returns the first eccentricity of the reference ellipsoid; otherwise, the same assumptions as above pertain.

Summary
Four new algorithms are presented that draw loxodromes, great ellipses, curves of alignment, and normal sections.The algorithms are parameterized so that a simple implementation based on vector algebra is available, given a supporting geometric geodesy library to calculate the coordinate system.There are many, more or less equivalent, choices for the conversion from geocentric Cartesian coordinates to geodetic longitude, latitude, and height [17,18].However, for long lines, the chord between the beginning and ending points can come near to the geocenter, which can cause some algorithms to fail.For such situations, a robust algorithm such as [19] is recommended.Also, the algorithms herein are language neutral and do not depend on external packages.However, it should be noted that efficient algorithms and packages exist for geodesics [20] when the actual geodesic is required.
The examples herein were computed using the WGS84 reference ellipsoid.Ref. [21] showed that the minute difference in the flattening of these ellipsoids cannot result in a difference in coordinates at the millimeter level.So, these examples are properly illustrative for GRS80, as well.s = Limit of s (0) , s (1) , s (2) , . . .ϕ = asin(s) where: s (0) = tanh(q) s (i+1) = tanh(q + ϵatanh(ϵ s (i) ) and asin(. ..) is the inverse sine function.The iteration is avoidable by using truncated infinite series: "first calculate the conformal latitude chi based on the inverse of the spherical Mercator, then use the method found in https://arxiv.org/pdf/2212.05818.pdf to calculate geodetic latitude" (accessed 17 April 2024) (anonymous reviewer 2024).

Figure 1 .
Figure 1.The curve of alignment (orange) and the geodesic (blue) from (165 E, 40 S) and (0 E, 45 N) on the WGS 84 reference ellipsoid.The Equator appears in green and the Prime Meridian in orange.

Figure 1 .
Figure 1.The curve of alignment (orange) and the geodesic (blue) from (165 • E, 40 • S) and (0 • E, 45 • N) on the WGS 84 reference ellipsoid.The Equator appears in green and the Prime Meridian in orange.

Figure 2 .
Figure 2. The construction of the (shorter limb of the) normal section from A to B.Figure 2. The construction of the (shorter limb of the) normal section from A to B.

Figure 2 .
Figure 2. The construction of the (shorter limb of the) normal section from A to B.Figure 2. The construction of the (shorter limb of the) normal section from A to B.

Figure 3 .
Figure 3.The two normal sections (black) and the geodesic (blue) between (0° E, 45° N) and (165° E, 40° S) on the WGS 84 reference ellipsoid.The Equator appears in green and the Prime Meridian in orange.
3. A point of the straight-line segment (chord) from A to B is given in parametric form by () =  +  ( − ) for 0 ≤  ≤ 1.

Figure 3 .
Figure 3.The two normal sections (black) and the geodesic (blue) between (0 • E, 45 • N) and (165 • E, 40 • S) on the WGS 84 reference ellipsoid.The Equator appears in green and the Prime Meridian in orange.

Figure 5 .
Figure 5.The geodesic (blue) and loxodrome (black) from (0° E, 45° N) to (165° E, 40° S) on the WGS 84 reference ellipsoid.The Equator appears in green and the Prime Meridian in orange.

Figure 5 .
Figure 5.The geodesic (blue) and loxodrome (black) from (0 • E, 45 • N) to (165 • E, 40 • S) on the WGS 84 reference ellipsoid.The Equator appears in green and the Prime Meridian in orange.
t: floating point scalar in [0, 1] A := (λ A , ϕ A ) # starting point a := Mf(A) # forward Mercator mapping formulas B := (λ B , ϕ B ) # ending point b := Mf(B) # forward Mercator mapping formulas v := b − a # vector subtraction p t := a + tv # linear interpolation p L := Mr(p t ) # reverse Mercator mapping formula • N) to (165 • E, 40 • S) on the WGS 84 reference ellipsoid.The arc length of the geodesic is 18,669,335.840m, and the arc length of the loxodrome is 19,066,164.69m, as computed by recursive subdivision, which is 396,828.850m longer.