This is the command GeodSolve that can be run in the OnWorks free hosting provider using one of our multiple free online workstations such as Ubuntu Online, Fedora Online, Windows online emulator or MAC OS online emulator

**PROGRAM:**

**NAME**

GeodSolve -- perform geodesic calculations

**SYNOPSIS**

**GeodSolve**[

**-i**|

**-l**

__lat1__

__lon1__

__azi1__] [

**-a**] [

**-e**

__a__

__f__]

**-u**] [

**-d**|

**-:**] [

**-w**] [

**-b**] [

**-f**

] [

**-p**

__prec__] [

**-E**] [

**--comment-delimiter**

__commentdelim__] [

**--version**|

**-h**|

**--help**] [

**--input-file**

__infile__|

**--input-string**

__instring__] [

**--line-separator**

__linesep__] [

**--output-file**

__outfile__]

**DESCRIPTION**

The shortest path between two points on the 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

forward azimuths

__azi1__and

__azi2__at the two end points.

**GeodSolve**operates in one of three modes:

1. By default,

**GeodSolve**accepts lines on the standard input containing

__lat1__

__lon1__

__azi1__

__s12__and prints

__lat2__

__lon2__

__azi2__on standard output. This is the direct geodesic

calculation.

2. Command line arguments

**-l**

__lat1__

__lon1__

__azi1__specify a geodesic line.

**GeodSolve**then

accepts a sequence of

__s12__values (one per line) on standard input and prints

__lat2__

__lon2__

__azi2__for each. This generates a sequence of points on a single geodesic.

3. With the

**-i**command line argument,

**GeodSolve**performs the inverse geodesic

calculation. It reads lines containing

__lat1__

__lon1__

__lat2__

__lon2__and prints the

corresponding values of

__azi1__

__azi2__

__s12__.

**OPTIONS**

**-i**perform an inverse geodesic calculation (see 3 above).

**-l**line mode (see 2 above); generate a sequence of points along the geodesic specified by

__lat1__

__lon1__

__azi1__. The

**-w**flag can be used to swap the default order of the 2 geographic

coordinates, provided that it appears before

**-l**.

**-a**arc mode; on input

__and__output

__s12__is replaced by

__a12__the arc length (in degrees) on

the auxiliary sphere. See "AUXILIARY SPHERE".

**-e**specify the ellipsoid via

__a__

__f__; the equatorial radius is

__a__and the flattening is

__f__.

Setting

__f__= 0 results in a sphere. Specify

__f__< 0 for a prolate ellipsoid. A simple

fraction, e.g., 1/297, is allowed for

__f__. By default, the WGS84 ellipsoid is used,

__a__=

6378137 m,

__f__= 1/298.257223563.

**-u**unroll the longitude. Normally, on output longitudes are reduced to lie in

[-180deg,180deg). However with this option, the returned longitude

__lon2__is "unrolled"

so that

__lon2__-

__lon1__indicates how often and in what sense the geodesic has encircled

the earth. Use the

**-f**option, to get both longitudes printed.

**-d**output angles as degrees, minutes, seconds instead of decimal degrees.

**-:**like

**-d**, except use : as a separator instead of the d, ', and " delimiters.

**-w**on input and output, longitude precedes latitude (except that, on input, this can be

overridden by a hemisphere designator,

__N__,

__S__,

__E__,

__W__).

**-b**report the

__back__azimuth at point 2 instead of the forward azimuth.

**-f**full output; each line of output consists of 12 quantities:

__lat1__

__lon1__

__azi1__

__lat2__

__lon2__

__azi2__

__s12__

__a12__

__m12__

__M12__

__M21__

__S12__.

__a12__is described in "AUXILIARY SPHERE". The four

quantities

__m12__,

__M12__,

__M21__, and

__S12__are described in "ADDITIONAL QUANTITIES".

**-p**set the output precision to

__prec__(default 3);

__prec__is the precision relative to 1 m.

See "PRECISION".

**-E**use "exact" algorithms (based on elliptic integrals) for the geodesic calculations.

These are more accurate than the (default) series expansions for |

__f__| > 0.02.

**--comment-delimiter**

set the comment delimiter to

__commentdelim__(e.g., "#" or "//"). If set, the input

lines will be scanned for this delimiter and, if found, the delimiter and the rest of

the line will be removed prior to processing and subsequently appended to the output

line (separated by a space).

**--version**

print version and exit.

**-h**print usage and exit.

**--help**

print full documentation and exit.

**--input-file**

read input from the file

__infile__instead of from standard input; a file name of "-"

stands for standard input.

**--input-string**

read input from the string

__instring__instead of from standard input. All occurrences

of the line separator character (default is a semicolon) in

__instring__are converted to

newlines before the reading begins.

**--line-separator**

set the line separator character to

__linesep__. By default this is a semicolon.

**--output-file**

write output to the file

__outfile__instead of to standard output; a file name of "-"

stands for standard output.

**INPUT**

**GeodSolve**measures all angles in degrees and all lengths (

__s12__) in meters, and all areas

(

__S12__) in meters^2. On input angles (latitude, longitude, azimuth, arc length) can be as

decimal degrees or degrees, minutes, seconds. For example, "40d30", "40d30'", "40:30",

"40.5d", and 40.5 are all equivalent. By default, latitude precedes longitude for each

point (the

**-w**flag switches this convention); however on input either may be given first

by appending (or prepending)

__N__or

__S__to the latitude and

__E__or

__W__to the longitude. Azimuths

are measured clockwise from north; however this may be overridden with

__E__or

__W__.

For details on the allowed formats for angles, see the "GEOGRAPHIC COORDINATES" section of

__GeoConvert__(1).

**AUXILIARY** **SPHERE**

Geodesics on the ellipsoid can be transferred to the

__auxiliary__

__sphere__on which the

distance is measured in terms of the arc length

__a12__(measured in degrees) instead of

__s12__.

In terms of

__a12__, 180 degrees is the distance from one equator crossing to the next or from

the minimum latitude to the maximum latitude. Geodesics with

__a12__> 180 degrees do not

correspond to shortest paths. With the

**-a**flag,

__s12__(on both input and output) is

replaced by

__a12__. The

**-a**flag does

__not__affect the full output given by the

**-f**flag (which

always includes both

__s12__and

__a12__).

**ADDITIONAL** **QUANTITIES**

The

**-f**flag reports four additional quantities.

The reduced length of the geodesic,

__m12__, is defined such that if the initial azimuth is

perturbed by d

__azi1__(radians) then the second point is displaced by

__m12__d

__azi1__in the

direction perpendicular to the geodesic.

__m12__is given in meters. On a curved surface the

reduced length obeys a symmetry relation,

__m12__+

__m21__= 0. On a flat surface, we have

__m12__=

__s12__.

__M12__and

__M21__are geodesic scales. If two geodesics are parallel at point 1 and separated

by a small distance

__dt__, then they are separated by a distance

__M12__

__dt__at point 2.

__M21__is

defined similarly (with the geodesics being parallel to one another at point 2).

__M12__and

__M21__are dimensionless quantities. On a flat surface, we have

__M12__=

__M21__= 1.

If points 1, 2, and 3 lie on a single geodesic, then the following addition 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.

Finally,

__S12__is the area between the geodesic from point 1 to point 2 and the equator;

i.e., it is the area, measured counter-clockwise, of the geodesic quadrilateral with

corners (

__lat1__,

__lon1__), (0,

__lon1__), (0,

__lon2__), and (

__lat2__,

__lon2__). It is given in meters^2.

**PRECISION**

__prec__gives precision of the output with

__prec__= 0 giving 1 m precision,

__prec__= 3 giving 1

mm precision, etc.

__prec__is the number of digits after the decimal point for lengths. For

decimal degrees, the number of digits after the decimal point is

__prec__+ 5. For DMS

(degree, minute, seconds) output, the number of digits after the decimal point in the

seconds component is

__prec__+ 1. The minimum value of

__prec__is 0 and the maximum is 10.

**ERRORS**

An illegal line of input will print an error message to standard output beginning with

"ERROR:" and causes

**GeodSolve**to return an exit code of 1. However, an error does not

cause

**GeodSolve**to terminate; following lines will be converted.

**ACCURACY**

Using the (default) series solution, GeodSolve is accurate to about 15 nm (15 nanometers)

for the WGS84 ellipsoid. The approximate maximum error (expressed as a distance) for an

ellipsoid with the same major radius as the WGS84 ellipsoid and different values of the

flattening is

|f| error

0.01 25 nm

0.02 30 nm

0.05 10 um

0.1 1.5 mm

0.2 300 mm

If

**-E**is specified, GeodSolve is accurate to about 40 nm (40 nanometers) for the WGS84

ellipsoid. The approximate maximum error (expressed as a distance) for an ellipsoid with

a quarter meridian of 10000 km and different values of the

__a/b__= 1 -

__f__is

1-f error (nm)

1/128 387

1/64 345

1/32 269

1/16 210

1/8 115

1/4 69

1/2 36

1 15

2 25

4 96

8 318

16 985

32 2352

64 6008

128 19024

**MULTIPLE** **SOLUTIONS**

The shortest distance returned for 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__.

**EXAMPLES**

Route from JFK Airport to Singapore Changi Airport:

echo 40:38:23N 073:46:44W 01:21:33N 103:59:22E |

GeodSolve -i -: -p 0

003:18:29.9 177:29:09.2 15347628

Waypoints on the route at intervals of 2000km:

for ((i = 0; i <= 16; i += 2)); do echo ${i}000000;done |

GeodSolve -l 40:38:23N 073:46:44W 003:18:29.9 -: -p 0

40:38:23.0N 073:46:44.0W 003:18:29.9

58:34:45.1N 071:49:36.7W 004:48:48.8

76:22:28.4N 065:32:17.8W 010:41:38.4

84:50:28.0N 075:04:39.2E 150:55:00.9

67:26:20.3N 098:00:51.2E 173:27:20.3

49:33:03.2N 101:06:52.6E 176:07:54.3

31:34:16.5N 102:30:46.3E 177:03:08.4

13:31:56.0N 103:26:50.7E 177:24:55.0

04:32:05.7S 104:14:48.7E 177:28:43.6

Use GeodSolve online using onworks.net services