Hostname: page-component-745bb68f8f-s22k5 Total loading time: 0 Render date: 2025-02-06T14:47:22.579Z Has data issue: false hasContentIssue false
  • English
  • Français

Circle Approximations on Spheroids

Published online by Cambridge University Press:  12 September 2011

Young Joon Ahn ,
Jian Cui  and
Christoph Hoffmann
Young Joon Ahn*
Affiliation:
(Chosun University, South Korea)
Jian Cui
Affiliation:
(Purdue University, USA)
Christoph Hoffmann
Affiliation:
(Purdue University, USA)
*
(Email: ahn@Chosun.ac.kr)
Rights & Permissions [Opens in a new window]

Abstract

We present an approximation method for geodesic circles on a spheroid. Our ap­proximation curve is the intersection of two spheroids whose axes are parallel, and it interpolates four points of the geodesic circle. Our approximation method has two merits. One is that the approximation curve can be obtained algebraically, and the other is that the approximation error is very small. For example, our approximation of a circle of radius 1000 km on the Earth has error 1·13 cm or less. We analyze the error of our approximation using the Hausdorff distance and confirm it by a geodesic distance computation.

Keywords

geodesic circlespheroidgeodesic distancegeodesic path
Type
Research Article
Information
The Journal of Navigation , Volume 64 , Issue 4 , October 2011 , pp. 739 - 749
Copyright
Copyright © The Royal Institute of Navigation 2011

1. INTRODUCTION

Global navigation, on the seas and in the air, approximates the earth's surface as an oblate spheroid. Spatial databases, from Microsoft and IBM, do likewise. Consequently, computing shortest distances and shortest paths on spheroids are of practical interest. Given two points on the spheroid, a shortest path between them is a geodesic, and geodesics constitute the basis for defining distances on the surface. Likewise, a circle on a spheroid would appropriately be defined as the locus of all points equidistant from a given point P, the centre of the circle.

The mathematics of geodesics on the spheroid is well understood; e.g., [Reference Kallay8,Reference Rollins13,Reference Schechter14,Reference Sjöberg16]. However, computing geodesics is sufficiently complicated that approximations of geodesics are of practical interest and are used in navigation algorithms.

The geodesy on spheroids has seen more recent development. Vector solutions for great ellipses, azimuth angles and the intersection of two circle of equal altitude on spheroids are presented in [Reference Earle4–Reference Gonzalez6]. The intersection of two geodesic paths on the sphere or the ellipsoid are calculated and approximated in [Reference Sjöberg15,Reference Sjöberg17]. The azimuth or arc length of a great ellipse, loxodrome, and geodesic path are analyzed and compared in [Reference Pallikaris and Latsas11,Reference Schechter14,Reference Tseng and Lee18]. Rollins presents an integral formula for the geodesic distance that is not an improper integral [Reference Rollins13].

Circles on spheroids arise as points at fixed distance from a given centre point. They can be used when enclosing a region on the spheroid by a buffer of fixed distance. A buffer could be used to define territorial water boundaries and is a generalization of Euclidean offsets in the plane or in Euclidean 3-space. Such offsets would be defined as a (one-sided) envelope, on the surface of the spheroid, of a family of fixed-radius circles centred on the boundary of the territory.

In this paper, we offer a particularly simple approximation of circles on the spheroid. We represent the circles as intersections with a second spheroid whose dimensions and coordinates are derived from the coordinates of the circle's centre and radius. This approximation is simple to compute and easy to evaluate using standard intersection algorithms; e.g., [Reference Hoffmann7,Reference Patrikalakis and Maekawa12]. We evaluate the quality of such an approximation in this note.

2. PRIOR WORK AND DEFINITIONS

We consider the spheroid with major axis a and minor axis b, centred at the origin with the minor axis on the z-axis:

(1)
$$S : \displaystyle{{x^2} \over {a^2}} + \displaystyle{{y^2} \over {a^2}} + \displaystyle{{z^2} \over {b^2}} = 1$$

The point P=(x, y, z) on the surface of the unit sphere is usually represented parametrically, using longitude λ and latitude ϕ, as:

(2)
$$x(\lambda, {\varphi} ) = {\cos} \ \lambda \ {\cos} \ \varphi, y(\lambda, {\varphi} ) = {\sin} \ \lambda \ {\cos} \ \varphi, z(\lambda, \varphi ) = {\sin} \ \varphi $$

where we assume that the equator is the xy-plane and the point (1, 0, 0) is the intersection of the null-meridian with the equator. The coordinates of the point P are also the coordinates of the unit normal of the unit sphere at P. A point P s=(xs,ys,zs) on the spheroid has the unit normal:

(3)
$$N\left( {x_s, y_s, z_s} \right) = \displaystyle{{\left( {{\raise0.7ex\hbox{${x_s} $} \!\mathord{\left/ {\vphantom {{x_s} {a^2,}}} \right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${a^2,} $}}{\raise0.7ex\hbox{${y_s} $} \!\mathord{\left/ {\vphantom {{y_s} {a^2,}}} \right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${a^2,} $}}{\raise0.7ex\hbox{${z_s} $} \!\mathord{\left/ {\vphantom {{z_s} {b^2,}}} \right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${b^2} $}}} \right)} \over {\sqrt {{\raise0.7ex\hbox{${x_s^2} $} \!\mathord{\left/ {\vphantom {{x_s^2} {a^4}}} \right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${a^4} $}} + {\raise0.7ex\hbox{${y_s^2} $} \!\mathord{\left/ {\vphantom {{y_s^2} {a^4}}} \right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${a^4} $}} + {\raise0.7ex\hbox{${z_s^2} $} \!\mathord{\left/ {\vphantom {{z_s^2} {b^4}}} \right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${b^4} $}}}}} $$

thus, we can parameterize the points of the spheroid; e.g, [Reference Kallay9].

In [Reference Rollins13] a simplified integral is derived for computing the geodesic between two points on a spheroid numerically. The formulation depends on Clairaut's constant that is to be determined, also numerically, before the geodesic can be integrated. Reference [Reference Sjöberg16] proposes a method to compute the constant numerically.

The fact that a numerical computation is needed, coupled with the fact that the eccentricity of the planet is slight, motivates approximating geodesics, and Kallay gives an elegant solution in [Reference Kallay9]. Kallay's construction is as follows: Given two points P 1 and P 2 on the spheroid with normals N 1 and N 2, consider the geodesic between N 1 and N 2 on the Gauss sphere, a greatest circle. Using the inverse Gauss map, this greatest circle is mapped to the spheroid and constitutes the approximate geodesic. Kallay's approximate geodesic can be shown to be a greatest ellipse on the spheroid, that is, the intersection of the spheroid with a plane through the origin; [Reference Kallay9]. This also implies that the approximate geodesic can be extended and that a subdivision of the approximate geodesic is also an approximate geodesic.

In this paper, we propose an approximation that samples the geodesic circle at four distinguished points and constructs an axis-parallel spheroid that passes through those four points. We note that for approximating circles on an oblate spheroid the intersecting spheroid is prolate, and, conversely, if the circle of a prolate spheroid is to be approximated, then the intersecting spheroid is oblate. The intersection of the oblate spheroid and prolate spheroid then serves as an approximation of the geodesic circle. The approximation is remarkably ac­curate, even for highly eccentric spheroids.

Let P be a point on the spheroid. For its latitude we use the reduced latitude β, as shown in Figure 1 (left). It is related to the geodesic latitude ϕ by: $\beta = \arctan \left( {\sqrt {1 - \varepsilon} \tan \varphi} \right)$, where $\varepsilon = \sqrt {1 - \left( {{\raise0.7ex\hbox{$b$} \!\mathord{\left/ {\vphantom {b a}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$a$}}} \right)^2} $ is the eccentricity of the spheroid.

Figure 1. Left: Normal angle ϕ and reduced latitude β. Right: Azimuth α at point P of the geodesic path between two point P and Q.

When following a path from P, the direction in which we head is quantified by α. The angle α is the azimuth of geodesic path at the starting point P. It is the deviation from the meridian through P, with zero the direction to the (north) pole. See also Figure 1 (right). Rollins [Reference Rollins13] presents formulae of proper integrals which solve the direct and indirect problems for geodesic paths connecting two points or azimuth headings.

The Clairaut constant is the cosine of the maximum reduced latitude of a geodesic path,

(4)
$$c = {\cos} \ \beta _{{\rm max}} = {\cos} \ \beta\ {\sin} \ \alpha, $$

where α and β are the reduced latitude and azimuth along the geodesic path. Sjoberg [Reference Sjöberg16] gives a numerical method to determine the constant for a geodesic path connecting two given points. He points out that, having determined the constant, accurate relations between the latitude, azimuth and longitude can be computed at any point along the geodesic.

3. APPROXIMATION OF GEODESIC CIRCLES

In this section, we present our approximation of geodesic circles on a spheroid by a Spheroid-Spheroid Intersection (SSI) curve. The geodesic circle on a spheroid S of radius r and centred at the point x0=(x 0,y0,z0)∊S is the set of points x=(x, y, z)∊S such that:

(5)
$$d_g ({\bf x},{\rm} {\bf x}_0 ) = r$$

where d g( x, x0) is the geodesic distance between two points x and x0. For each point x on the geodesic circle, there is a geodesic path connecting two points x and x0, which has an azimuth α at the starting point x0, as shown in Figure 1. Thus each point x depends on the heading azimuth α,

(6)
$${\bf x} = {\bf x}\left( \alpha \right),{\rm} 0 \les\alpha \les{\rm 2}\pi. $$

Figure 2 shows the geodesic circles on a spheroid obtained numerically by vary­ing α from 0 to 2π.

Figure 2. Spheroid $S = \displaystyle{{x^2} \over {20^2}} + \displaystyle{{y^2} \over {20^2}} + \displaystyle{{z^2} \over {10^2}} = 1$ and the five geodesic circles centred at β 0=0°, 30°, 60°, 80° and 90° with geodesic radius=2.

The geodesic circle is not a planar curve unless the centre x0 lies on the North or the South Pole. The two lines x(0)x(π) and ${\bf x}\left( {\displaystyle{\pi \over 2}} \right)\cdot {\bf x}\left( {\displaystyle{{3\pi} \over 2}} \right)$ are orthogonal and skew. Thus there are two parallel planes containing the lines. We call the distance between the two parallel planes the thickness of the geodesic circle x(α). Figure 3 shows the parallel planes with a bounding box of the required thickness.

Figure 3. The geodesic circles (green) x(α) and our approximation curves (red) r(z) with the bounding box (blue) for a=20, b=10, and geodesic radius r=2.

Table 1 shows the thickness of a sample of circles on a spheroid with major axis 20 and minor axis 10. When decreasing the radius of the circle by one half, the thickness decreases to almost 1/4, as shown in the table. The smaller is the radius of geodesic circle, the closer the decreasing ratio is to one quarter.

Table 1. Thickness of geodesic circles x(α) of radius r=2 centred at β 0 for a=20 and b=10. The last row represents the decreasing ratio of thickness when the radius decreases by half.

In our approximation method, the approximation curve is the intersection of the spheroid S with another spheroid S′ chosen to interpolate the geodesic circle at the points ${\bf x}\left( {\displaystyle{k\pi \over 2}} \right)$, where k=0, …, 3; see Figure 2. The new spheroid S′ is centred at x1=(x1,y1,z1) and has the major and minor axes R/a and R/b:

(7)
$$S' : \displaystyle{{\left( {x - x_1} \right)^2} \over {\left( {{\raise0.7ex\hbox{$R$} \!\mathord{\left/ {\vphantom {R a}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$a$}}} \right)^2}} + \displaystyle{{\left( {y - y_1} \right)^2} \over {\left( {{\raise0.7ex\hbox{$R$} \!\mathord{\left/ {\vphantom {R a}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$a$}}} \right)^2}} + \displaystyle{{\left( {z - z_1} \right)^2} \over {\left( {{\raise0.7ex\hbox{$R$} \!\mathord{\left/ {\vphantom {R b}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$b$}}} \right)^2}} = 1$$

Note that the ratio of the major and minor axes of S′ is inverse of that of S, as shown in Figure 4. Thus, if S is oblate, then S′ is prolate and vice versa. If S is a sphere, then so is S′. Both spheroids S and S′ are symmetric with respect to the plane through the three points, x0, the centre of S′, and the origin. Thus the longitude λ0 of the point x0 is equal to that of the point x1. Cutting the spheroids S and S′ by the ρz-plane of symmetry, two ellipses E and E′ on the ρz-plane are obtained:

(8)
$$E:\displaystyle{{\rho ^2} \over {a^2}} + \displaystyle{{z^2} \over {b^2}} = 1\;{\rm and}\;E':\displaystyle{{\left( {\rho - \rho _1} \right)^2} \over {\left( {{\raise0.7ex\hbox{$R$} \!\mathord{\left/ {\vphantom {R a}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$a$}}} \right)^2}} + \displaystyle{{\left( {z - z_1} \right)^2} \over {\left( {{\raise0.7ex\hbox{$R$} \!\mathord{\left/ {\vphantom {R b}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$b$}}} \right)^2}} = 1$$

where $\rho = \sqrt {x^2 + y^2} $ and $\rho _i = \sqrt {x_i^2 + y_i^2} $, i=0, 1, as shown in Figure 5.

Figure 4. The major and minor axes of the prolate (elongated) spheroid S′ are R/a=5·57 and R/b=11·14 when a=20, b=10, β 0=0◦, and r=2.

Figure 5. Two ellipses are the cross-sections of two spheroids S (blue) and S′ (red) with a plane passing through the origin, x0 and the north pole. The geodesic circle x(α) is centred at (ρ 0,z0) and the spheroid S′ is centred at (ρ 1,z1). The two ellipses have intersection points at x(0) and x(π)(green points).

We present the formulae of coefficients of the spheroid S′ as follows:

Proposition 3.1.

The spheroid S′ in Equation (7) passing through four points, ${\bf x}\left( {\displaystyle{{k\pi} \over 2}} \right)$, k=0, …, 3, is unique unless the geodesic circle is centred at the North or South Pole, and its coefficients are:

(9)
$$z_1 = \left\{ {\matrix{ {\left( {{\raise0.7ex\hbox{${A_3} $} \!\mathord{\left/ {\vphantom {{A_3} 2}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$2$}} - A_0 A_2} \right)/ \left\{ {\left( {{\raise0.7ex\hbox{$b$} \!\mathord{\left/ {\vphantom {b a}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$a$}}} \right)^2 \left. {\left( {z\displaystyle{\pi \over 2}} \right) - z\left( 0 \right)} \right) + \left. {A_1 A_2} \right\Bigg\}\quad if\,\beta _0 \ne 0} \right.} \cr {0 \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ \, if\,\beta _0 = 0}\hfill \cr}} \right.$$
(10)
$$\rho _1 = \left\{ {\matrix{ {A_1 z_1 + A_0 \qquad\qquad\qquad\qquad\qquad\qquad\qquad if\,\beta _0 \ne 0}\hfill \cr {\left\{ {x\left( {\displaystyle{\pi \over 2}} \right)^2 + y\left( {\displaystyle{\pi \over 2}} \right)^2 - \rho \left( 0 \right)^2 - \left( {{\raise0.7ex\hbox{$b$} \!\mathord{\left/ {\vphantom {b a}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$a$}}} \right)^2 \left. {z\left( 0 \right)^2} \right\Bigg\}/2A_2 \,\,\,\,if\,\beta _0 = 0} \right.} \cr}} \right.$$
(11)
$$R = \sqrt {a^2 \left( {\rho \left( 0 \right) - \rho _1} \right)^2 + b^2 \left( {z\left( 0 \right) - z_1} \right)^2} $$

where:

(12)
$$\eqalign{& A_0 = \displaystyle{{b^2 \,\left( {z\left( \pi \right) - z\left( 0 \right)} \right)\left( {z\left( \pi \right) + z\left( 0 \right)} \right)} \over {a^2 \left( {\rho \left( \pi \right) - \rho \left( 0 \right)} \right)2}} + \displaystyle{{\rho \left( \pi \right) + \rho \left( 0 \right)} \over 2} \cr & A_1 = - \displaystyle{{b^2 \left( {z\left( \pi \right) - z\left( 0 \right)} \right)} \over {a^2 \left( {\rho \left( \pi \right) - \rho \left( 0 \right)} \right)}} \cr & A_2 = x\left( {\displaystyle{\pi \over 2}} \right)\cos \lambda _0 + y\left( {\displaystyle{\pi \over 2}} \right)\sin \lambda _0 - \rho \left( 0 \right) \cr & A_3 = x\left( {\displaystyle{\pi \over 2}} \right)^2 + y\left( {\displaystyle{\pi \over 2}} \right)^2 + \displaystyle{{b^2} \over {a^2}} z\left( {\displaystyle{\pi \over 2}} \right)^2 - \left( {x\left( 0 \right)^2 + y\left( 0 \right)^2 + \displaystyle{{b^2} \over {a^2}} z\left( 0 \right)^2} \right)} $$

Proof.

Let the centre of the geodesic circle lie on the equator, i.e., β 0=0. By symmetry of the geodesic circle with respect to the equator, z 1 and $z\left( {\displaystyle{\pi \over 2}} \right)$ are zero. Since the ellipse E′ in Equation (8) interpolates the points (ρ(0), z(0)), we obtain Equation (11). Substituting the point x(π/2) into Equation (7) and using:

(13)
$$x_1 {\rm} = {\rm} \rho _1 {\rm cos}\left( {\lambda _0} \right),y_1 {\rm} = {\rm} \rho _1 {\rm sin}\left( {\lambda _0} \right){\rm} $$

and Equation (11), we have Equation (10).

Let β 0≠0. Substituting the point (x(0), z(0)) into Equation (8), we obtain Equation (11), and R depends on (ρ 1,z1). Since E′ interpolates two points (ρ(0),z(0)) and (ρ(π),z(π)), Equation (8) yields in this case:

(14)
$$a^2 \left( {\rho \left( 0 \right) - \rho _1} \right)^2 + b^2\left( {z\left( 0 \right) - z_1} \right)^2 = a^2 \left( {\rho \left( \pi \right) - \rho _1} \right)^2 + b^2\left( {z\left( \pi \right) - z_1} \right)^2 $$

which implies Equation (10), where ρ 1 depends on z 1. Inserting the point x(π/2) into Equation (7) and using Equations (10), (11) and (13), we have:

(15)
$$\displaystyle{{b^2} \over {a^2}} \left( {z\left( {\displaystyle{\pi \over 2}} \right) - z\left( 0 \right)} \right)z_1 + A_2 \rho _1 - \displaystyle{{A_3} \over 2} = 0$$

which is a linear equation of z 1, since ρ 1=A1z1 + A0 is too. Thus it can be solved easily to yield Equation (9). □

The coefficients x 1 and y 1 are obtained by Equation (13) from Proposition 3.1. Table 2 lists the coefficients (R/a, ρ 1,z1) of spheroid S′ for a=20, b=10. Now, the parametric equation of intersection of two spheroids S and S′ can be obtained algebraically, which is our approximation curve for the geodesic circle.

Table 2. The coefficients (R/a, ρ1,z1) of the spheroid S′ when the geodesic circle centred at β0=0°, 30°, 60°, 80°, with radius r for a=20, b=10.

Proposition 3.2.

The approximation curve of the geodesic circle is:

(16)
$$\left( {x,{\rm} y,{\rm} z} \right) = \left( {X(z} \right)\ {\rm cos}\left( {\lambda _0} \right){\rm} - {\rm} Y(z) \ {\sin}\left( {\lambda _0} \right),X(z) \ {\sin}\left( {\lambda _0} \right) + Y(z) \ {\cos}\left( {\lambda _0} \right),z)$$

for z(π)⩽z⩽z(0), where:

(17)
$$X\left( z \right) = \displaystyle{1 \over {2\rho _1}} \left( {\displaystyle{{b^2} \over {a^2}} \left( {z - z_1} \right)^2 - \displaystyle{{a^2} \over {b^2}} z^2 + a^2 - \displaystyle{{R^2} \over {a^2}} + \rho _1^2} \right)$$
(18)
$$Y\left( z \right) = \pm \sqrt {a^2 - \displaystyle{{a^2} \over {b^2}} z^2 - X\left( z \right)^2} $$

Proof.

Without loss of generality we assume the circle is centred on the 0-meridian. Equivalently, we rotate the axes of x, y by angle λ0 with respect to the z-axis. Let (X,Y,z) be the new coordinates.

(19)
$$\left( {\matrix{ x \cr y \cr}} \right) = \left( {\matrix{ {\cos \left( {\lambda _0} \right)\ \ - \sin \left( {\lambda _0} \right)} \cr {\sin \left( {\lambda _0} \right) \ \ \cos \left( {\lambda _0} \right)} \cr}} \right)\left( {\matrix{ X \cr Y \cr}} \right)$$

Then the spheroids S and S′ in Equations (1) and (7) yield:

(20)
$$X^2 + Y^2 + \displaystyle{{a^2} \over {b^2}} z^2 = a^2 $$
(21)
$$\left( {X - \rho _1} \right)^2 + Y^2 + \displaystyle{{b^2} \over {a^2}} \left( {z - z_1} \right)^2 = \displaystyle{{R^2} \over {a^2}} $$

Subtracting Equation (20) from Equation (21) we get Equations (17) and (18). Therefore the assertion follows from Equations (17), (18), and (19), which is the equation for the intersection of two spheroids S and S′. □

Let r(z)=(x(z),y(z),z), for z(π) ⩽ z ⩽ z(0), be our approximation curve. As an example, we approximate several geodesic circles on the spheroid with a=20, b=10. For the geodesic circles x(α) centred at β 0=0° , 30° , 60° and 80°, with the radius r=2, which are shown green in Figure 2, their approximation curves r(z) plotted by red colour together, as shown in Figure 3.

We use the Hausdorff distance d H (x, r) to quantify the maximum difference between two curves x(α) and r(z). The Hausdorff distance is widely used to measure approximation error in Computer Aided Geometric Design, as a fair measure for the distance between two curves or surfaces; e.g., [Reference Ahn1–Reference Barton, Hanniel, Elver and Kim3,Reference Kim, Oh, Yoon, Kim and Elber10]. We sample points along the curves. For the sample points x(α), α=j°, j=0, …, M=180, we calculate the Hausdorff distance numerically by:

(22)
$$d_H \left( \bf{x,r} \right) = \mathop {\max} \limits_{\,j = 0, \cdots, M} \mathop {\min} \limits_{z\left( \pi \right) \les z \les z\left( 0 \right)} \left| { {\bf x}\left( {\,j\,\ring} \right) - { \bf r}\left( z \right)} \right|$$

The resulting distance values d H (x, r) for the approximations of the geodesic circles, centred at β 0=0° , 30° , 60°, 80° and with radius r=2, 1, 0·5, 0·25, are shown in Table 3. We observe that the approximation error decreases with ratios between $\displaystyle{1 \over {7{\cdot}09}}$ and $\displaystyle{1 \over {13{\cdot}55}}$.

Table 3. The Hausdorff distance dH (x, r) between the geodesic circle x(α) and our approximation curve r(z) for a=20, b=10.

4. QUALITY OF APPROXIMATION FOR GEOGRAPHIC APPLICATIONS

In this section, we use our method to find geodesic circles on the earth's surface. Customarily, the earth is considered an oblate spheroid with axes a=6378·137 km and b=6356·7523 km. We computed the thickness of the geodesic circles of radius r=1000 km, 500 km, 250 km, 125 km, a measure of the non-planarity, shown in Table 4. The decreasing ratios are between ${\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {3{\cdot}96}}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${3{\cdot}96}$}}$ and ${\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {4{\cdot}00}}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${4{\cdot}00}$}}$, and the smaller is the radius of geodesic circle, the closer the decreasing ratio is to a quarter, as was the case for the highly eccentric spheroid in Table 1.

Table 4. Thickness of the geometric circle of radius r=1000 km, 500 km, 250 km, 125 km, and centered at β 0=0°, 45°, on the earth, a=6378·137 km and b=6356·7523 km.

For the geodesic circles centred on the equator (β 0=0°) and with radii r=1000∼ 125 km, the Hausdorff distance measuring the approximation error for our curve r(z) is 1·13 cm or less, as shown in the middle column of Table 5. Also, for the geodesic circle centred at β 0=45°, the Hausdorff distance is 0·54 cm or less, as shown in the third column of Table 5. The approximation accuracy is remarkable.

Table 5. Hausdorff distance dH (x, r) on the earth (a=6378·137 km and b=6356·7523 km).

The decreasing error ratios are between ${\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {7{\cdot}83}}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${7{\cdot}83}$}}$ and ${\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {10{\cdot}38}}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${10{\cdot}38}$}}$, and the smaller is the radius of the geodesic circle, the closer the decreasing ratio is to ${\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 8}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$8$}}$, similar to the results in Table 3.

To confirm our numerical results, we compute the geodesic distance d g(x0, r(z)) from the approximation curve r(z), z(π)⩽z⩽z(0) to the centre x0. For the sample points:

(23)
$$z = z\left( \varsigma \right) = \displaystyle{{z\left( 0 \right) + z\left( \pi \right)} \over 2} + \displaystyle{{\left( {z\left( 0 \right) - z\left( \pi \right)} \right)\cos \left( \varsigma \right)} \over 2}$$

ζ=(j −0·5)π/M, j=1, …, M, we plot the error function d g(x0, r(ζ)) −r as shown in Figure 6. The fact that the maximum differences $\mathop {\max} \limits_{0 \les\varsigma \les\pi} \left| {d_g \left( {{\bf x}_0, {\bf r}\left( z \right)} \right) - r} \right|$ in Table 6 are the almost same as the Hausdorff distances d H (x, r) in Table 5 corroborates the values of Table 6.

Figure 6. On the earth, (a=6378·137 km and b=6356·7523 km), for β 0=0°(left column), 45°(right column), and r=1000, 500, 250, 125 km (from top to bottom), the error function as Equation (23).

Table 6. The maximum difference $\mathop {\max} \limits_{0 \les\varsigma \les\pi} \left| {d_g \left( {{\bf x}_0, {\bf r}\left( \rm z \right)} \right) - r} \right|$ between the radius r and the geodesic distance dg(x0, r(z)) from approximation curve r(z) to the centre x0.

ACKNOWLEDGEMENT

This study was supported by research funds from Chosun University, 2007.

References

REFERENCES

[1]Ahn, Y. J.. (2005) Helix approximation with conic and quadratic Bezier curves. Comp. Aided Geom. Desi., 22, 551–565.CrossRefGoogle Scholar
[2]Ahn, Y. J.. (2010) Approximation of conic sections by curvature continuous quartic Bezier curves. Comp. Math. Appl., 60, 1986–1993.CrossRefGoogle Scholar
[3]Barton, M., Hanniel, I., Elver, G., and Kim, M.-S.. (2010) Precise Hausdorff distance computation between polygonal meshes. Comp. Aided Geom. Desi., 27, 580–591.CrossRefGoogle Scholar
[4]Earle, M. A.. (2000). A vector solution for navigation on a great ellipse. The Journal of Navigation 53, 473–481.CrossRefGoogle Scholar
[5]Earle, M. A.. (2008). Vector solution for azimuth. The Journal of Navigation 61, 537–545.CrossRefGoogle Scholar
[6]Gonzalez, A. R.. (2008). Vector solution for the intersection of two circles of equal altitude. The Journal of Navigation, 61, 355–365.CrossRefGoogle Scholar
[7]Hoffmann, C. M.. (1989). Geometric and Solid Modeling. Morgan Kaufmann, San Mateo, CA, 1989. http://www.cs.purdue.edu/homes/cmh/distribution/books/geo.html.Google Scholar
[8]Kallay, M.. (2007). Defining edges on a round earth. In ACM GIS '07, 2007. Article 63.CrossRefGoogle Scholar
[9]Kallay, M.. (2008). Geometric algorithms on an ellipsoid earth model. In ACM GIS '08, 2008. Article 42.CrossRefGoogle Scholar
[10]Kim, Y.-J., Oh, Y.-T., Yoon, S.-H., Kim, M.-S., and Elber, G.. (2010). Precise Haussdorff distance computation for planar freeform curves using biarcs and depth buffer. The Visual Computer, 26, 1007–1016.CrossRefGoogle Scholar
[11]Pallikaris, A. and Latsas, G.. (2009). New algorithm for great elliptic sailing (GES). The Journal of Navigation 62, 497–507.CrossRefGoogle Scholar
[12]Patrikalakis, N. M. and Maekawa, T.. (2002). Shape Interrogation for Computer Aided Design and Manufacturing. Springer Verlag, New York.CrossRefGoogle Scholar
[13]Rollins, C.. An integral for geodesic length. (2010). Survey Review, 42, 20–26.CrossRefGoogle Scholar
[14]Schechter, M.. (2007). Which way is Jerusalem? Navigating on a spheroid. The College Mathematics Journal, 38, 96–105.CrossRefGoogle Scholar
[15]Sjöberg, L.. (2002). Intersection on the sphere and ellipsoid. J. Geod, 76, 115–120.Google Scholar
[16]Sjöberg, L.. (2007). Precise determination of the Clairaut constant in ellipsoidal geodesy. Survey Review, 39, 81–86.CrossRefGoogle Scholar
[17]Sjöberg, L.. (2008). Geodetic intersection on the ellipsoid. J. Geod, 82, 565–567.CrossRefGoogle Scholar
[18]Tseng, W.-K. and Lee, H.-S.. (2010). Navigation on a great ellipse. J. Mari. Scie. Tech., 18, 369–375.Google Scholar
Figure 0

Figure 1. Left: Normal angle ϕ and reduced latitude β. Right: Azimuth α at point P of the geodesic path between two point P and Q.

Figure 1

Figure 2. Spheroid $S = \displaystyle{{x^2} \over {20^2}} + \displaystyle{{y^2} \over {20^2}} + \displaystyle{{z^2} \over {10^2}} = 1$ and the five geodesic circles centred at β0=0°, 30°, 60°, 80° and 90° with geodesic radius=2.

Figure 2

Figure 3. The geodesic circles (green) x(α) and our approximation curves (red) r(z) with the bounding box (blue) for a=20, b=10, and geodesic radius r=2.

Figure 3

Table 1. Thickness of geodesic circles x(α) of radius r=2 centred at β0 for a=20 and b=10. The last row represents the decreasing ratio of thickness when the radius decreases by half.

Figure 4

Figure 4. The major and minor axes of the prolate (elongated) spheroid S′ are R/a=5·57 and R/b=11·14 when a=20, b=10, β0=0◦, and r=2.

Figure 5

Figure 5. Two ellipses are the cross-sections of two spheroids S (blue) and S′ (red) with a plane passing through the origin, x0 and the north pole. The geodesic circle x(α) is centred at (ρ0,z0) and the spheroid S′ is centred at (ρ1,z1). The two ellipses have intersection points at x(0) and x(π)(green points).

Figure 6

Table 2. The coefficients (R/a, ρ1,z1) of the spheroid S′ when the geodesic circle centred at β0=0°, 30°, 60°, 80°, with radius r for a=20, b=10.

Figure 7

Table 3. The Hausdorff distance dH (x, r) between the geodesic circle x(α) and our approximation curve r(z) for a=20, b=10.

Figure 8

Table 4. Thickness of the geometric circle of radius r=1000 km, 500 km, 250 km, 125 km, and centered at β0=0°, 45°, on the earth, a=6378·137 km and b=6356·7523 km.

Figure 9

Table 5. Hausdorff distance dH (x, r) on the earth (a=6378·137 km and b=6356·7523 km).

Figure 10

Figure 6. On the earth, (a=6378·137 km and b=6356·7523 km), for β0=0°(left column), 45°(right column), and r=1000, 500, 250, 125 km (from top to bottom), the error function as Equation (23).

Figure 11

Table 6. The maximum difference $\mathop {\max} \limits_{0 \les\varsigma \les\pi} \left| {d_g \left( {{\bf x}_0, {\bf r}\left( \rm z \right)} \right) - r} \right|$ between the radius r and the geodesic distance dg(x0, r(z)) from approximation curve r(z) to the centre x0.