net.sf.geographiclib

## Class Geodesic

• ```public class Geodesic
extends Object```
Geodesic calculations.

The shortest path between two points on a ellipsoid at (lat1, lon1) and (lat2, lon2) is called the geodesic. Its length is s12 and the geodesic from point 1 to point 2 has azimuths azi1 and azi2 at the two end points. (The azimuth is the heading measured clockwise from north. azi2 is the "forward" azimuth, i.e., the heading that takes you beyond point 2 not back to point 1.)

Given lat1, lon1, azi1, and s12, we can determine lat2, lon2, and azi2. This is the direct geodesic problem and its solution is given by the function `Direct`. (If s12 is sufficiently large that the geodesic wraps more than halfway around the earth, there will be another geodesic between the points with a smaller s12.)

Given lat1, lon1, lat2, and lon2, we can determine azi1, azi2, and s12. This is the inverse geodesic problem, whose solution is given by `Inverse`. Usually, the solution to the inverse problem is unique. In cases where there are multiple solutions (all with the same s12, of course), all the solutions can be easily generated once a particular solution is provided.

The standard way of specifying the direct problem is the specify the distance s12 to the second point. However it is sometimes useful instead to specify the arc length a12 (in degrees) on the auxiliary sphere. This is a mathematical construct used in solving the geodesic problems. The solution of the direct problem in this form is provided by `ArcDirect`. An arc length in excess of 180° indicates that the geodesic is not a shortest path. In addition, the arc length between an equatorial crossing and the next extremum of latitude for a geodesic is 90°.

This class can also calculate several other quantities related to geodesics. These are:

• reduced length. If we fix the first point and increase azi1 by dazi1 (radians), the second point is displaced m12 dazi1 in the direction azi2 + 90°. The quantity m12 is called the "reduced length" and is symmetric under interchange of the two points. On a curved surface the reduced length obeys a symmetry relation, m12 + m21 = 0. On a flat surface, we have m12 = s12. The ratio s12/m12 gives the azimuthal scale for an azimuthal equidistant projection.
• geodesic scale. Consider a reference geodesic and a second geodesic parallel to this one at point 1 and separated by a small distance dt. The separation of the two geodesics at point 2 is M12 dt where M12 is called the "geodesic scale". M21 is defined similarly (with the geodesics being parallel at point 2). On a flat surface, we have M12 = M21 = 1. The quantity 1/M12 gives the scale of the Cassini-Soldner projection.
• area. The area between the geodesic from point 1 to point 2 and the equation is represented by S12; it is the area, measured counter-clockwise, of the geodesic quadrilateral with corners (lat1,lon1), (0,lon1), (0,lon2), and (lat2,lon2). It can be used to compute the area of any simple geodesic polygon.

The quantities m12, M12, M21 which all specify the behavior of nearby geodesics obey addition rules. If points 1, 2, and 3 all lie on a single geodesic, then the following rules hold:

• s13 = s12 + s23
• a13 = a12 + a23
• S13 = S12 + S23
• m13 = m12 M23 + m23 M21
• M13 = M12 M23 − (1 − M12 M21) m23 / m12
• M31 = M32 M21 − (1 − M23 M32) m12 / m23

The results of the geodesic calculations are bundled up into a `GeodesicData` object which includes the input parameters and all the computed results, i.e., lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, M12, M21, S12.

The functions `Direct`, `ArcDirect`, and `Inverse` include an optional final argument outmask which allows you specify which results should be computed and returned. If you omit outmask, then the "standard" geodesic results are computed (latitudes, longitudes, azimuths, and distance). outmask is bitor'ed combination of `GeodesicMask` values. For example, if you wish just to compute the distance between two points you would call, e.g.,

``` ```
GeodesicData g = Geodesic.WGS84.Inverse(lat1, lon1, lat2, lon2,

Additional functionality is provided by the `GeodesicLine` class, which allows a sequence of points along a geodesic to be computed.

The shortest distance returned by the solution of the inverse problem is (obviously) uniquely defined. However, in a few special cases there are multiple azimuths which yield the same shortest distance. Here is a catalog of those cases:

• lat1 = −lat2 (with neither point at a pole). If azi1 = azi2, the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting [azi1, azi2] → [azi2, azi1], [M12, M21] → [M21, M12], S12 → −S12. (This occurs when the longitude difference is near ±180° for oblate ellipsoids.)
• lon2 = lon1 ± 180° (with neither point at a pole). If azi1 = 0° or ±180°, the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting [ azi1, azi2] → [−azi1, −azi2], S12 → − S12. (This occurs when lat2 is near −lat1 for prolate ellipsoids.)
• Points 1 and 2 at opposite poles. There are infinitely many geodesics which can be generated by setting [azi1, azi2] → [azi1, azi2] + [d, −d], for arbitrary d. (For spheres, this prescription applies when points 1 and 2 are antipodal.)
• s12 = 0 (coincident points). There are infinitely many geodesics which can be generated by setting [azi1, azi2] → [azi1, azi2] + [d, d], for arbitrary d.

The calculations are accurate to better than 15 nm (15 nanometers) for the WGS84 ellipsoid. See Sec. 9 of arXiv:1102.1215v1 for details. The algorithms used by this class are based on series expansions using the flattening f as a small parameter. These are only accurate for |f| < 0.02; however reasonably accurate results will be obtained for |f| < 0.2. Here is a table of the approximate maximum error (expressed as a distance) for an ellipsoid with the same equatorial radius as the WGS84 ellipsoid and different values of the flattening.

```     |f|      error
0.01     25 nm
0.02     30 nm
0.05     10 um
0.1     1.5 mm
0.2     300 mm ```

The algorithms are described in

Example of use:

``` ```
// Solve the direct geodesic problem.

// This program reads in lines with lat1, lon1, azi1, s12 and prints
// out lines with lat2, lon2, azi2 (for the WGS84 ellipsoid).

import java.util.*;
import net.sf.geographiclib.*;
public class Direct {
public static void main(String[] args) {
try {
Scanner in = new Scanner(System.in);
double lat1, lon1, azi1, s12;
while (true) {
lat1 = in.nextDouble(); lon1 = in.nextDouble();
azi1 = in.nextDouble(); s12 = in.nextDouble();
GeodesicData g = Geodesic.WGS84.Direct(lat1, lon1, azi1, s12);
System.out.println(g.lat2 + " " + g.lon2 + " " + g.azi2);
}
}
catch (Exception e) {}
}
}``````
• ### Field Summary

Fields
Modifier and Type Field and Description
`static Geodesic` `WGS84`
A global instantiation of Geodesic with the parameters for the WGS84 ellipsoid.
• ### Constructor Summary

Constructors
Constructor and Description
```Geodesic(double a, double f)```
Constructor for a ellipsoid with
• ### Method Summary

All Methods
Modifier and Type Method and Description
`GeodesicData` ```ArcDirect(double lat1, double lon1, double azi1, double a12)```
Solve the direct geodesic problem where the length of the geodesic is specified in terms of arc length.
`GeodesicData` ```ArcDirect(double lat1, double lon1, double azi1, double a12, int outmask)```
Solve the direct geodesic problem where the length of the geodesic is specified in terms of arc length and with a subset of the geodesic results returned.
`GeodesicLine` ```ArcDirectLine(double lat1, double lon1, double azi1, double a12)```
Define a `GeodesicLine` in terms of the direct geodesic problem specified in terms of arc length with all capabilities included.
`GeodesicLine` ```ArcDirectLine(double lat1, double lon1, double azi1, double a12, int caps)```
Define a `GeodesicLine` in terms of the direct geodesic problem specified in terms of arc length with a subset of the capabilities included.
`GeodesicData` ```Direct(double lat1, double lon1, double azi1, boolean arcmode, double s12_a12, int outmask)```
The general direct geodesic problem.
`GeodesicData` ```Direct(double lat1, double lon1, double azi1, double s12)```
Solve the direct geodesic problem where the length of the geodesic is specified in terms of distance.
`GeodesicData` ```Direct(double lat1, double lon1, double azi1, double s12, int outmask)```
Solve the direct geodesic problem where the length of the geodesic is specified in terms of distance and with a subset of the geodesic results returned.
`GeodesicLine` ```DirectLine(double lat1, double lon1, double azi1, double s12)```
Define a `GeodesicLine` in terms of the direct geodesic problem specified in terms of distance with all capabilities included.
`GeodesicLine` ```DirectLine(double lat1, double lon1, double azi1, double s12, int caps)```
Define a `GeodesicLine` in terms of the direct geodesic problem specified in terms of distance with a subset of the capabilities included.
`double` `EllipsoidArea()`
`double` `Flattening()`
`GeodesicLine` ```GenDirectLine(double lat1, double lon1, double azi1, boolean arcmode, double s12_a12, int caps)```
Define a `GeodesicLine` in terms of the direct geodesic problem specified in terms of either distance or arc length with a subset of the capabilities included.
`GeodesicData` ```Inverse(double lat1, double lon1, double lat2, double lon2)```
Solve the inverse geodesic problem.
`GeodesicData` ```Inverse(double lat1, double lon1, double lat2, double lon2, int outmask)```
Solve the inverse geodesic problem with a subset of the geodesic results returned.
`GeodesicLine` ```InverseLine(double lat1, double lon1, double lat2, double lon2)```
Define a `GeodesicLine` in terms of the inverse geodesic problem with all capabilities included.
`GeodesicLine` ```InverseLine(double lat1, double lon1, double lat2, double lon2, int caps)```
Define a `GeodesicLine` in terms of the inverse geodesic problem with a subset of the capabilities included.
`GeodesicLine` ```Line(double lat1, double lon1, double azi1)```
Set up to compute several points on a single geodesic with all capabilities included.
`GeodesicLine` ```Line(double lat1, double lon1, double azi1, int caps)```
Set up to compute several points on a single geodesic with a subset of the capabilities included.
`double` `MajorRadius()`
• ### Methods inherited from class java.lang.Object

`equals, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait`
• ### Field Detail

• #### WGS84

`public static final Geodesic WGS84`
A global instantiation of Geodesic with the parameters for the WGS84 ellipsoid.
• ### Constructor Detail

• #### Geodesic

```public Geodesic(double a,
double f)```
Constructor for a ellipsoid with

Parameters:
`a` - equatorial radius (meters).
`f` - flattening of ellipsoid. Setting f = 0 gives a sphere. Negative f gives a prolate ellipsoid.
Throws:
`GeographicErr` - if a or (1 − f ) a is not positive.
• ### Method Detail

• #### Direct

```public GeodesicData Direct(double lat1,
double lon1,
double azi1,
double s12)```
Solve the direct geodesic problem where the length of the geodesic is specified in terms of distance.

Parameters:
`lat1` - latitude of point 1 (degrees).
`lon1` - longitude of point 1 (degrees).
`azi1` - azimuth at point 1 (degrees).
`s12` - distance between point 1 and point 2 (meters); it can be negative.
Returns:
a `GeodesicData` object with the following fields: lat1, lon1, azi1, lat2, lon2, azi2, s12, a12.

lat1 should be in the range [−90°, 90°]. The values of lon2 and azi2 returned are in the range [−180°, 180°].

If either point is at a pole, the azimuth is defined by keeping the longitude fixed, writing lat = ±(90° − ε), and taking the limit ε → 0+. An arc length greater that 180° signifies a geodesic which is not a shortest path. (For a prolate ellipsoid, an additional condition is necessary for a shortest path: the longitudinal extent must not exceed of 180°.)

• #### Direct

```public GeodesicData Direct(double lat1,
double lon1,
double azi1,
double s12,
Solve the direct geodesic problem where the length of the geodesic is specified in terms of distance and with a subset of the geodesic results returned.

Parameters:
`lat1` - latitude of point 1 (degrees).
`lon1` - longitude of point 1 (degrees).
`azi1` - azimuth at point 1 (degrees).
`s12` - distance between point 1 and point 2 (meters); it can be negative.
`outmask` - a bitor'ed combination of `GeodesicMask` values specifying which results should be returned.
Returns:
a `GeodesicData` object with the fields specified by outmask computed.

lat1, lon1, azi1, s12, and a12 are always included in the returned result. The value of lon2 returned is in the range [−180°, 180°], unless the outmask includes the `GeodesicMask.LONG_UNROLL` flag.

• #### ArcDirect

```public GeodesicData ArcDirect(double lat1,
double lon1,
double azi1,
double a12)```
Solve the direct geodesic problem where the length of the geodesic is specified in terms of arc length.

Parameters:
`lat1` - latitude of point 1 (degrees).
`lon1` - longitude of point 1 (degrees).
`azi1` - azimuth at point 1 (degrees).
`a12` - arc length between point 1 and point 2 (degrees); it can be negative.
Returns:
a `GeodesicData` object with the following fields: lat1, lon1, azi1, lat2, lon2, azi2, s12, a12.

lat1 should be in the range [−90°, 90°]. The values of lon2 and azi2 returned are in the range [−180°, 180°].

If either point is at a pole, the azimuth is defined by keeping the longitude fixed, writing lat = ±(90° − ε), and taking the limit ε → 0+. An arc length greater that 180° signifies a geodesic which is not a shortest path. (For a prolate ellipsoid, an additional condition is necessary for a shortest path: the longitudinal extent must not exceed of 180°.)

• #### ArcDirect

```public GeodesicData ArcDirect(double lat1,
double lon1,
double azi1,
double a12,
Solve the direct geodesic problem where the length of the geodesic is specified in terms of arc length and with a subset of the geodesic results returned.

Parameters:
`lat1` - latitude of point 1 (degrees).
`lon1` - longitude of point 1 (degrees).
`azi1` - azimuth at point 1 (degrees).
`a12` - arc length between point 1 and point 2 (degrees); it can be negative.
`outmask` - a bitor'ed combination of `GeodesicMask` values specifying which results should be returned.
Returns:
a `GeodesicData` object with the fields specified by outmask computed.

lat1, lon1, azi1, and a12 are always included in the returned result. The value of lon2 returned is in the range [−180°, 180°], unless the outmask includes the `GeodesicMask.LONG_UNROLL` flag.

• #### DirectLine

```public GeodesicLine DirectLine(double lat1,
double lon1,
double azi1,
double s12)```
Define a `GeodesicLine` in terms of the direct geodesic problem specified in terms of distance with all capabilities included.

Parameters:
`lat1` - latitude of point 1 (degrees).
`lon1` - longitude of point 1 (degrees).
`azi1` - azimuth at point 1 (degrees).
`s12` - distance between point 1 and point 2 (meters); it can be negative.
Returns:
a `GeodesicLine` object.

This function sets point 3 of the GeodesicLine to correspond to point 2 of the direct geodesic problem.

lat1 should be in the range [−90°, 90°].

• #### DirectLine

```public GeodesicLine DirectLine(double lat1,
double lon1,
double azi1,
double s12,
int caps)```
Define a `GeodesicLine` in terms of the direct geodesic problem specified in terms of distance with a subset of the capabilities included.

Parameters:
`lat1` - latitude of point 1 (degrees).
`lon1` - longitude of point 1 (degrees).
`azi1` - azimuth at point 1 (degrees).
`s12` - distance between point 1 and point 2 (meters); it can be negative.
`caps` - bitor'ed combination of `GeodesicMask` values specifying the capabilities the GeodesicLine object should possess, i.e., which quantities can be returned in calls to `GeodesicLine.Position`.
Returns:
a `GeodesicLine` object.

This function sets point 3 of the GeodesicLine to correspond to point 2 of the direct geodesic problem.

lat1 should be in the range [−90°, 90°].

• #### ArcDirectLine

```public GeodesicLine ArcDirectLine(double lat1,
double lon1,
double azi1,
double a12)```
Define a `GeodesicLine` in terms of the direct geodesic problem specified in terms of arc length with all capabilities included.

Parameters:
`lat1` - latitude of point 1 (degrees).
`lon1` - longitude of point 1 (degrees).
`azi1` - azimuth at point 1 (degrees).
`a12` - arc length between point 1 and point 2 (degrees); it can be negative.
Returns:
a `GeodesicLine` object.

This function sets point 3 of the GeodesicLine to correspond to point 2 of the direct geodesic problem.

lat1 should be in the range [−90°, 90°].

• #### ArcDirectLine

```public GeodesicLine ArcDirectLine(double lat1,
double lon1,
double azi1,
double a12,
int caps)```
Define a `GeodesicLine` in terms of the direct geodesic problem specified in terms of arc length with a subset of the capabilities included.

Parameters:
`lat1` - latitude of point 1 (degrees).
`lon1` - longitude of point 1 (degrees).
`azi1` - azimuth at point 1 (degrees).
`a12` - arc length between point 1 and point 2 (degrees); it can be negative.
`caps` - bitor'ed combination of `GeodesicMask` values specifying the capabilities the GeodesicLine object should possess, i.e., which quantities can be returned in calls to `GeodesicLine.Position`.
Returns:
a `GeodesicLine` object.

This function sets point 3 of the GeodesicLine to correspond to point 2 of the direct geodesic problem.

lat1 should be in the range [−90°, 90°].

• #### GenDirectLine

```public GeodesicLine GenDirectLine(double lat1,
double lon1,
double azi1,
boolean arcmode,
double s12_a12,
int caps)```
Define a `GeodesicLine` in terms of the direct geodesic problem specified in terms of either distance or arc length with a subset of the capabilities included.

Parameters:
`lat1` - latitude of point 1 (degrees).
`lon1` - longitude of point 1 (degrees).
`azi1` - azimuth at point 1 (degrees).
`arcmode` - boolean flag determining the meaning of the s12_a12.
`s12_a12` - if arcmode is false, this is the distance between point 1 and point 2 (meters); otherwise it is the arc length between point 1 and point 2 (degrees); it can be negative.
`caps` - bitor'ed combination of `GeodesicMask` values specifying the capabilities the GeodesicLine object should possess, i.e., which quantities can be returned in calls to `GeodesicLine.Position`.
Returns:
a `GeodesicLine` object.

This function sets point 3 of the GeodesicLine to correspond to point 2 of the direct geodesic problem.

lat1 should be in the range [−90°, 90°].

• #### Inverse

```public GeodesicData Inverse(double lat1,
double lon1,
double lat2,
double lon2)```
Solve the inverse geodesic problem.

Parameters:
`lat1` - latitude of point 1 (degrees).
`lon1` - longitude of point 1 (degrees).
`lat2` - latitude of point 2 (degrees).
`lon2` - longitude of point 2 (degrees).
Returns:
a `GeodesicData` object with the following fields: lat1, lon1, azi1, lat2, lon2, azi2, s12, a12.

lat1 and lat2 should be in the range [−90°, 90°]. The values of azi1 and azi2 returned are in the range [−180°, 180°].

If either point is at a pole, the azimuth is defined by keeping the longitude fixed, writing lat = ±(90° − ε), taking the limit ε → 0+.

The solution to the inverse problem is found using Newton's method. If this fails to converge (this is very unlikely in geodetic applications but does occur for very eccentric ellipsoids), then the bisection method is used to refine the solution.

• #### InverseLine

```public GeodesicLine InverseLine(double lat1,
double lon1,
double lat2,
double lon2)```
Define a `GeodesicLine` in terms of the inverse geodesic problem with all capabilities included.

Parameters:
`lat1` - latitude of point 1 (degrees).
`lon1` - longitude of point 1 (degrees).
`lat2` - latitude of point 2 (degrees).
`lon2` - longitude of point 2 (degrees).
Returns:
a `GeodesicLine` object.

This function sets point 3 of the GeodesicLine to correspond to point 2 of the inverse geodesic problem.

lat1 and lat2 should be in the range [−90°, 90°].

• #### InverseLine

```public GeodesicLine InverseLine(double lat1,
double lon1,
double lat2,
double lon2,
int caps)```
Define a `GeodesicLine` in terms of the inverse geodesic problem with a subset of the capabilities included.

Parameters:
`lat1` - latitude of point 1 (degrees).
`lon1` - longitude of point 1 (degrees).
`lat2` - latitude of point 2 (degrees).
`lon2` - longitude of point 2 (degrees).
`caps` - bitor'ed combination of `GeodesicMask` values specifying the capabilities the GeodesicLine object should possess, i.e., which quantities can be returned in calls to `GeodesicLine.Position`.
Returns:
a `GeodesicLine` object.

This function sets point 3 of the GeodesicLine to correspond to point 2 of the inverse geodesic problem.

lat1 and lat2 should be in the range [−90°, 90°].

• #### Line

```public GeodesicLine Line(double lat1,
double lon1,
double azi1)```
Set up to compute several points on a single geodesic with all capabilities included.

Parameters:
`lat1` - latitude of point 1 (degrees).
`lon1` - longitude of point 1 (degrees).
`azi1` - azimuth at point 1 (degrees).
Returns:
a `GeodesicLine` object.

lat1 should be in the range [−90°, 90°]. The full set of capabilities is included.

If the point is at a pole, the azimuth is defined by keeping the lon1 fixed, writing lat1 = ±(90 − ε), taking the limit ε → 0+.

`public double MajorRadius()`
Returns:
a the equatorial radius of the ellipsoid (meters). This is the value used in the constructor.
• #### Flattening

`public double Flattening()`
Returns:
f the flattening of the ellipsoid. This is the value used in the constructor.
• #### EllipsoidArea

`public double EllipsoidArea()`
Returns:
total area of ellipsoid in meters2. The area of a polygon encircling a pole can be found by adding EllipsoidArea()/2 to the sum of S12 for each side of the polygon.