diff --git a/doc/faq.xml b/doc/faq.xml index dabef747862..d64010e6a06 100644 --- a/doc/faq.xml +++ b/doc/faq.xml @@ -97,13 +97,13 @@ SELECT MAX(ST_NPoints(geom)) FROM sometable; - You can store point, line, polygon, multipoint, multiline, - multipolygon, and geometrycollections. In PostGIS 2.0 and above you can also store TINS and Polyhedral Surfaces in the basic geometry type. + You can store Point, LineString, Polygon, MultiPoint, MultiLineString, + MultiPolygon, and GeometryCollection geometries. In PostGIS 2.0 and above you can also store TINS and Polyhedral Surfaces in the basic geometry type. These are specified in the Open - GIS Well Known Text Format (with XYZ,XYM,XYZM extensions). There are three data types currently supported. + GIS Well Known Text Format (with Z, M, and ZM extensions). There are three data types currently supported. The standard OGC geometry data type which uses a planar coordinate system for measurement, the - geography data type which uses a geodetic coordinate system (not OGC, but you'll find a similar type in Microsoft SQL Server 2008+). Only WGS 84 long lat (SRID:4326) is supported - by the geography data type. The newest family member of the PostGIS spatial type family is raster for storing and analyzing raster data. Raster has its very own FAQ. Refer to + geography data type which uses a geodetic coordinate system, with calculations on either a sphere or spheroid. + The newest family member of the PostGIS spatial type family is raster for storing and analyzing raster data. Raster has its very own FAQ. Refer to and for more details. @@ -114,9 +114,9 @@ SELECT MAX(ST_NPoints(geom)) FROM sometable; - Short Answer: geography is a new data type that supports long range distances measurements, but most computations on it are currently + Short Answer: geography is a newer data type that supports long range distances measurements, but most computations on it are slower than they are on geometry. If - you use geography -- you don't need to learn much about planar coordinate systems. Geography is generally best + you use geography, you don't need to learn much about planar coordinate systems. Geography is generally best if all you care about is measuring distances and lengths and you have data from all over the world. Geometry data type is an older data type that has many more functions supporting it, enjoys greater support from third party tools, and operations on it are generally faster -- sometimes as much as 10 fold faster for larger geometries. diff --git a/doc/introduction.xml b/doc/introduction.xml index 3c999df3943..0a2b6f645bb 100644 --- a/doc/introduction.xml +++ b/doc/introduction.xml @@ -211,6 +211,7 @@ Markus Schaber, Maxime Guillaud, Maxime van Noppen, Michael Fuhr, +Mike Toews, Nathan Wagner, Nathaniel Clay, Nikita Shulga, diff --git a/doc/reference_constructor.xml b/doc/reference_constructor.xml index 1a775be40de..544cd38ddcc 100644 --- a/doc/reference_constructor.xml +++ b/doc/reference_constructor.xml @@ -174,8 +174,8 @@ SELECT ST_Box2dFromGeoHash('9qqj7nmxncgyy4d0dbxqz0', 0); Description - Returns a geography object from the well-known text or extended well-known representation. SRID 4326 is assumed. This - is an alias for ST_GeographyFromText. Points are always expressed in long lat form. + Returns a geography object from the well-known text or extended well-known representation. SRID 4326 is assumed if unspecified. + This is an alias for ST_GeographyFromText. Points are always expressed in long lat form. @@ -183,7 +183,10 @@ SELECT ST_Box2dFromGeoHash('9qqj7nmxncgyy4d0dbxqz0', 0); --- converting lon lat coords to geography ALTER TABLE sometable ADD COLUMN geog geography(POINT,4326); -UPDATE sometable SET geog = ST_GeogFromText('SRID=4326;POINT(' || lon || ' ' || lat || ')'); +UPDATE sometable SET geog = ST_GeogFromText('SRID=4326;POINT(' || lon || ' ' || lat || ')'); + +--- specify a geography point using EPSG:4267, NAD27 +SELECT ST_AsEWKT(ST_GeogFromText('SRID=4267;POINT(-77.0092 38.889588)')); @@ -207,7 +210,7 @@ UPDATE sometable SET geog = ST_GeogFromText('SRID=4326;POINT(' || lon || ' ' || Description - Returns a geography object from the well-known text representation. SRID 4326 is assumed. + Returns a geography object from the well-known text representation. SRID 4326 is assumed if unspecified. diff --git a/doc/reference_measure.xml b/doc/reference_measure.xml index 2b80364bb04..5f00ccfcfdd 100644 --- a/doc/reference_measure.xml +++ b/doc/reference_measure.xml @@ -643,8 +643,8 @@ SELECT ST_AsEWKT(ST_3DShortestLine(line,pt)) AS shl3d_line_pt, ST_Area - Returns the area of the surface if it is a polygon or - multi-polygon. For "geometry" type area is in SRID units. For "geography" area is in square meters. + Returns the area of the surface if it is a Polygon or + MultiPolygon. For geometry, a 2D Cartesian area is determined with units specified by the SRID. For geography, area is determined on a curved surface with units in square meters. @@ -663,12 +663,13 @@ SELECT ST_AsEWKT(ST_3DShortestLine(line,pt)) AS shl3d_line_pt, Description - Returns the area of the geometry if it is a polygon or - multi-polygon. Return the area measurement of an ST_Surface or - ST_MultiSurface value. For geometry Area is in the units of the srid. For geography area is in square meters and defaults to measuring about the spheroid of the geography (currently only WGS84). - To measure around the faster but less accurate sphere -- ST_Area(geog,false). + Returns the area of the geometry if it is a Polygon or + MultiPolygon. Return the area measurement of an ST_Surface or + ST_MultiSurface value. For geometry, a 2D Cartesian area is determined with units specified by the SRID. For geography, by default area is determined on a spheroid with units in square meters. + To measure around the faster but less accurate sphere, use ST_Area(geog,false). Enhanced: 2.0.0 - support for 2D polyhedral surfaces was introduced. + Enhanced: 2.2.0 - measurement on spheroid performed with GeographicLib for improved accuracy and robustness. &sfs_compliant; &sqlmm_compliant; SQL-MM 3: 8.1.2, 9.5.3 &P_support; @@ -680,8 +681,8 @@ SELECT ST_AsEWKT(ST_3DShortestLine(line,pt)) AS shl3d_line_pt, Examples Return area in square feet for a plot of Massachusetts land and multiply by conversion to get square meters. - Note this is in square feet because 2249 is - Mass State Plane Feet + Note this is in square feet because EPSG:2249 is + Massachusetts State Plane Feet SELECT ST_Area(the_geom) As sqft, ST_Area(the_geom)*POWER(0.3048,2) As sqm FROM (SELECT @@ -691,9 +692,9 @@ SELECT ST_Area(the_geom) As sqft, ST_Area(the_geom)*POWER(0.3048,2) As sqm ---------+------------- 928.625 | 86.27208552 -Return area square feet and transform to Massachusetts state plane meters (26986) to get square meters. +Return area square feet and transform to Massachusetts state plane meters (EPSG:26986) to get square meters. Note this is in square feet because 2249 is - Mass State Plane Feet and transformed area is in square meters since 26986 is state plane mass meters + Massachusetts State Plane Feet and transformed area is in square meters since EPSG:26986 is state plane Massachusetts meters SELECT ST_Area(the_geom) As sqft, ST_Area(ST_Transform(the_geom,26986)) As sqm @@ -705,7 +706,7 @@ SELECT ST_Area(the_geom) As sqft, ST_Area(ST_Transform(the_geom,26986)) As sqm 928.625 | 86.2724304199219 -Return area square feet and square meters using Geography data type. Note that we transform to our geometry to geography +Return area square feet and square meters using geography data type. Note that we transform to our geometry to geography (before you can do that make sure your geometry is in WGS 84 long lat 4326). Geography always measures in meters. This is just for demonstration to compare. Normally your table will be stored in geography data type already. @@ -720,9 +721,9 @@ SELECT ST_Area(the_geog)/POWER(0.3048,2) As sqft_spheroid, ST_Area(the_geog,fal ) ) ) As foo(the_geog); - sqft_spheroid | sqft_sphere | sqm_spheroid ------------------+------------------+------------------ -928.684405217197 | 927.186481558724 | 86.2776044452694 + sqft_spheroid | sqft_sphere | sqm_spheroid +------------------+------------------+------------------ + 928.684403538925 | 927.049336105925 | 86.2776042893529 --if your data is in geography already SELECT ST_Area(the_geog)/POWER(0.3048,2) As sqft, ST_Area(the_geog) As sqm @@ -758,13 +759,14 @@ SELECT ST_Area(the_geog)/POWER(0.3048,2) As sqft_spheroid, ST_Area(the_geog,fal Description Returns the azimuth in radians of the segment defined by the given - point-geometries, or NULL if the two points are coincident. The azimuth is north-based and is measured clockwise: North = 0; East = PI/2; South = PI; West = 3PI/2. - - The Azimuth is mathematical concept defined as the angle, in this case measured in radian, between a reference plane - and a point. + point geometries, or NULL if the two points are coincident. The azimuth is angle is referenced from north, and is positive clockwise: North = 0; East = π/2; South = π; West = 3π/2. + For the geography type, the forward azimuth is solved as part of the inverse geodesic problem. + The azimuth is mathematical concept defined as the angle between a reference plane and a point, with angular units in radians. + Units can be converted to degrees using a built-in PostgreSQL function degrees(), as shown in the example. Availability: 1.1.0 Enhanced: 2.0.0 support for geography was introduced. + Enhanced: 2.2.0 measurement on spheroid performed with GeographicLib for improved accuracy and robustness. Azimuth is especially useful in conjunction with ST_Translate for shifting an object along its perpendicular axis. See upgis_lineshift Plpgsqlfunctions PostGIS wiki section for example of this. @@ -773,13 +775,8 @@ SELECT ST_Area(the_geog)/POWER(0.3048,2) As sqft_spheroid, ST_Area(the_geog,fal Examples Geometry Azimuth in degrees -SELECT ST_Azimuth(ST_Point(25,45), ST_Point(75,100))/(2*pi())*360 as degA_B, - ST_Azimuth(ST_Point(75,100), ST_Point(25,45))/(2*pi())*360 As degB_A; - --- NOTE easier to remember syntax using PostgreSQL built-in degrees function -- --- Both yield same answer -- -SELECT degrees( ST_Azimuth(ST_Point(25,45), ST_Point(75,100)) ) as degA_B, - degrees( ST_Azimuth(ST_Point(75,100), ST_Point(25,45)) ) As degB_A; +SELECT degrees(ST_Azimuth(ST_Point(25, 45), ST_Point(75, 100))) AS degA_B, + degrees(ST_Azimuth(ST_Point(75, 100), ST_Point(25, 45))) AS degB_A; dega_b | degb_a ------------------+------------------ @@ -1933,8 +1930,8 @@ SELECT ST_Disjoint('POINT(0 0)'::geometry, 'LINESTRING ( 0 0, 0 2 )'::geometry); ST_Distance - For geometry type Returns the 2-dimensional cartesian minimum distance (based on spatial ref) between two geometries in - projected units. For geography type defaults to return spheroidal minimum distance between two geographies in meters. + For geometry type Returns the 2D Cartesian distance between two geometries in + projected units (based on spatial ref). For geography type defaults to return minimum geodesic distance between two geographies in meters. @@ -1975,18 +1972,19 @@ SELECT ST_Disjoint('POINT(0 0)'::geometry, 'LINESTRING ( 0 0, 0 2 )'::geometry); Description - For geometry type returns the 2-dimensional minimum cartesian distance between two geometries in - projected units (spatial ref units). For geography type defaults to return the minimum distance around WGS 84 spheroid between two geographies in meters. Pass in - false to return answer in sphere instead of spheroid. + For geometry type returns the minimum 2D Cartesian distance between two geometries in + projected units (spatial ref units). For geography type defaults to return the minimum geodesic distance between two geographies in meters. If use_spheroid is + false, a faster sphere calculation is used instead of a spheroid. &sfs_compliant; &sqlmm_compliant; SQL-MM 3: 5.1.23 &curve_support; - &sfcgal_enhanced; + &sfcgal_enhanced; Availability: 1.5.0 geography support was introduced in 1.5. Speed improvements for planar to better handle large or many vertex geometries - Enhanced: 2.1.0 improved speed for geography. See Making Geography faster for details. + Enhanced: 2.1.0 improved speed for geography. See Making Geography faster for details. Enhanced: 2.1.0 - support for curved geometries was introduced. + Enhanced: 2.2.0 - measurement on spheroid performed with GeographicLib for improved accuracy and robustness. @@ -2432,7 +2430,7 @@ FROM having the same SRID. For geography units are in meters and measurement is - defaulted to use_spheroid=true (measure around WGS 84 spheroid), for faster check, use_spheroid=false to measure along sphere. + defaulted to use_spheroid=true, for faster check, use_spheroid=false to measure along sphere. This function call will automatically include a bounding box @@ -2698,7 +2696,7 @@ t ST_Length - Returns the 2d length of the geometry if it is a linestring or multilinestring. geometry are in units of spatial reference and geography are in meters (default spheroid) + Returns the 2D length of the geometry if it is a LineString or MultiLineString. geometry are in units of spatial reference and geography are in meters (default spheroid) @@ -2716,9 +2714,12 @@ t Description - For geometry: Returns the cartesian 2D length of the geometry if it is a linestring, multilinestring, ST_Curve, ST_MultiCurve. 0 is returned for - areal geometries. For areal geometries use ST_Perimeter. Geometry: Measurements are in the units of the - spatial reference system of the geometry. Geography: Units are in meters and also acts as a Perimeter function for areal geogs. + For geometry: Returns the 2D Cartesian length of the geometry if it is a LineString, MultiLineString, ST_Curve, ST_MultiCurve. 0 is returned for + areal geometries. For areal geometries use . For geometry types, units for length measures are specified by the + spatial reference system of the geometry. + For geography types, the calculations are performed using the inverse geodesic problem, where length units are in meters. + If PostGIS is compiled with PROJ version 4.8.0 or later, the spheroid is specified by the SRID, otherwise it is exclusive to WGS84. + If use_spheroid=false, then calculations will approximate a sphere instead of a spheroid. Currently for geometry this is an alias for ST_Length2D, but this may change to support higher dimensions. Changed: 2.0.0 Breaking change -- in prior versions applying this to a MULTI/POLYGON of type geography would give you the perimeter of the POLYGON/MULTIPOLYGON. In 2.0.0 @@ -2727,13 +2728,13 @@ t &sfs_compliant; s2.1.5.1 &sqlmm_compliant; SQL-MM 3: 7.1.2, 9.3.4 Availability: 1.5.0 geography support was introduced in 1.5. - &sfcgal_enhanced; + &sfcgal_enhanced; Geometry Examples - Return length in feet for line string. Note this is in feet because 2249 is - Mass State Plane Feet + Return length in feet for line string. Note this is in feet because EPSG:2249 is + Massachusetts State Plane Feet SELECT ST_Length(ST_GeomFromText('LINESTRING(743238 2967416,743238 2967450,743265 2967450, 743265.625 2967416,743238 2967416)',2249)); @@ -2742,7 +2743,7 @@ st_length 122.630744000095 ---Transforming WGS 84 linestring to Massachusetts state plane meters +--Transforming WGS 84 LineString to Massachusetts state plane meters SELECT ST_Length( ST_Transform( ST_GeomFromEWKT('SRID=4326;LINESTRING(-72.1260 42.45, -72.1240 42.45666, -72.123 42.1546)'), @@ -2765,8 +2766,7 @@ FROM (SELECT ST_GeographyFromText( As foo; length_spheroid | length_sphere ------------------+------------------ - 34310.5703627305 | 34346.2060960742 -(1 row) + 34310.5703627288 | 34346.2060960742 @@ -2838,8 +2838,8 @@ FROM (SELECT ST_GeographyFromText( Examples - Return length in feet for a 3D cable. Note this is in feet because 2249 is - Mass State Plane Feet + Return length in feet for a 3D cable. Note this is in feet because EPSG:2249 is + Massachusetts State Plane Feet SELECT ST_3DLength(ST_GeomFromText('LINESTRING(743238 2967416 1,743238 2967450 1,743265 2967450 3, 743265.625 2967416 3,743238 2967416 3)',2249)); @@ -3342,7 +3342,7 @@ FROM (SELECT ST_Buffer(ST_GeomFromText('POINT(1 0.5)'), 3) As a, ST_Perimeter Return the length measurement of the boundary of an ST_Surface - or ST_MultiSurface geometry or geography. (Polygon, Multipolygon). geometry measurement is in units of spatial reference and geography is in meters. + or ST_MultiSurface geometry or geography. (Polygon, MultiPolygon). geometry measurement is in units of spatial reference and geography is in meters. @@ -3361,10 +3361,12 @@ FROM (SELECT ST_Buffer(ST_GeomFromText('POINT(1 0.5)'), 3) As a, Description - Returns the 2D perimeter of the geometry/geography if it is a ST_Surface, ST_MultiSurface (Polygon, Multipolygon). 0 is returned for - non-areal geometries. For linestrings use ST_Length. Measurements for geometry are in the units of the - spatial reference system of the geometry. Measurements for geography are in meters. If use_spheroid is set to false, then will - model earth as a sphere instead of a spheroid. + Returns the 2D perimeter of the geometry/geography if it is a ST_Surface, ST_MultiSurface (Polygon, MultiPolygon). 0 is returned for + non-areal geometries. For linear geometries use . For geometry types, units for perimeter measures are specified by the + spatial reference system of the geometry. + For geography types, the calculations are performed using the inverse geodesic problem, where perimeter units are in meters. + If PostGIS is compiled with PROJ version 4.8.0 or later, the spheroid is specified by the SRID, otherwise it is exclusive to WGS84. + If use_spheroid=false, then calculations will approximate a sphere instead of a spheroid. Currently this is an alias for ST_Perimeter2D, but this may change to support higher dimensions. @@ -3375,8 +3377,8 @@ FROM (SELECT ST_Buffer(ST_GeomFromText('POINT(1 0.5)'), 3) As a, Examples: Geometry - Return perimeter in feet for polygon and multipolygon. Note this is in feet because 2249 is - Mass State Plane Feet + Return perimeter in feet for Polygon and MultiPolygon. Note this is in feet because EPSG:2249 is + Massachusetts State Plane Feet SELECT ST_Perimeter(ST_GeomFromText('POLYGON((743238 2967416,743238 2967450,743265 2967450, 743265.625 2967416,743238 2967416))', 2249)); @@ -3404,7 +3406,7 @@ st_perimeter Examples: Geography - Return perimeter in meters and feet for polygon and multipolygon. Note this is geography (WGS 84 long lat) + Return perimeter in meters and feet for Polygon and MultiPolygon. Note this is geography (WGS 84 long lat) SELECT ST_Perimeter(geog) As per_meters, ST_Perimeter(geog)/0.3048 As per_ft FROM ST_GeogFromText('POLYGON((-71.1776848522251 42.3902896512902,-71.1776843766326 42.3903829478009, @@ -3415,7 +3417,7 @@ FROM ST_GeogFromText('POLYGON((-71.1776848522251 42.3902896512902,-71.1776843766 37.3790462565251 | 122.634666195949 --- Multipolygon example -- +-- MultiPolygon example -- SELECT ST_Perimeter(geog) As per_meters, ST_Perimeter(geog,false) As per_sphere_meters, ST_Perimeter(geog)/0.3048 As per_ft FROM ST_GeogFromText('MULTIPOLYGON(((-71.1044543107478 42.340674480411,-71.1044542869917 42.3406744369506, -71.1044553562977 42.340673886454,-71.1044543107478 42.340674480411)), @@ -3614,9 +3616,8 @@ SELECT ST_AsEWKT(ST_PointOnSurface(ST_GeomFromEWKT('LINESTRING(0 5 1, 0 0 1, 0 1 Description - Returns a POINT projected from a start point using an azimuth (bearing) measured in radians and distance measured in meters. - Distance, azimuth and projection are all aspects of the same operation, describing (or in the case of projection, constructing) the relationship between two points on the world. - The azimuth is sometimes called the heading or the bearing in navigation. It is measured relative to true north (azimuth zero). East is azimuth 90 (pi/2), south is azimuth 180 (pi), west is azimuth 270 (pi*1.5). + Returns a POINT projected from a start point using an azimuth (bearing) measured in radians and distance measured in meters. This is also called a direct geodesic problem. + The azimuth is sometimes called the heading or the bearing in navigation. It is measured relative to true north (azimuth zero). East is azimuth 90 (π/2), south is azimuth 180 (π), west is azimuth 270 (3π/2). The distance is given in meters. Availability: 2.0.0 @@ -3628,19 +3629,20 @@ SELECT ST_AsEWKT(ST_PointOnSurface(ST_GeomFromEWKT('LINESTRING(0 5 1, 0 0 1, 0 1 SELECT ST_AsText(ST_Project('POINT(0 0)'::geography, 100000, radians(45.0))); st_astext ------------------------------------------- - POINT(0.63523102912532 0.63947233472882) + st_astext +-------------------------------------------- + POINT(0.635231029125537 0.639472334729198) (1 row) - Example: Using radians - projected point 100,000 meters and bearing pi/4 (45 degrees) + Example: Using radians - projected point 100,000 meters and bearing pi/4 radians (45 degrees) SELECT ST_AsText(ST_Project('POINT(0 0)'::geography, 100000, pi()/4)); - st_astext ------------------------------------------- - POINT(0.63523102912532 0.63947233472882) + st_astext +-------------------------------------------- + POINT(0.635231029125537 0.639472334729198) (1 row) @@ -3852,7 +3854,7 @@ SELECT mat.name, pat.name, ST_RelateMatch(mat.val, pat.val) As satisfied If g1 and g2 are intersecting with more than one point the function will return a line with start and end in the same point but it can be any of the intersecting points. The line returned will always start in g1 and end in g2. - The length of the line this function returns will always be the same as st_distance returns for g1 and g2. + The length of the line this function returns will always be the same as ST_Distance returns for g1 and g2. Availability: 1.5.0 diff --git a/liblwgeom/Makefile.in b/liblwgeom/Makefile.in index fc27456c259..7edd3f987e9 100644 --- a/liblwgeom/Makefile.in +++ b/liblwgeom/Makefile.in @@ -92,6 +92,7 @@ SA_OBJS = \ varint.o NM_OBJS = \ + geodesic.o \ lwspheroid.o ifeq (@SFCGAL@,sfcgal) diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c index 6fc99fca10c..903e6c284d7 100644 --- a/liblwgeom/cunit/cu_geodetic.c +++ b/liblwgeom/cunit/cu_geodetic.c @@ -8,6 +8,12 @@ * This is free software; you can redistribute and/or modify it under * the terms of the GNU General Public Licence. See the COPYING file. * + * Note: Geodesic measurements have been independently verified using + * GeographicLib v 1.37 with MPFR C++ using utilities GeodSolve and + * Planimeter, with -E "exact" flag with the following env vars: + * export GEOGRAPHICLIB_DIGITS=1000 + * export WGS84_ELIPSOID="6378137 298.257223563" + * export WGS84_SPHERE="6371008.7714150598325213222 0" **********************************************************************/ #include @@ -77,13 +83,17 @@ static void test_sphere_direction(void) geographic_point_init(1, 0, &e); dist = sphere_distance(&s, &e); dir = sphere_direction(&s, &e, dist); - CU_ASSERT_DOUBLE_EQUAL(dir, M_PI / 2.0, 0.0001); + /* GeodSolve -i -E -p 16 -e 1 0 --input-string "0 0 0 1" */ + CU_ASSERT_DOUBLE_EQUAL(dir, M_PI_2, 1e-14); + CU_ASSERT_DOUBLE_EQUAL(dist, 0.0174532925199433, 1e-14); geographic_point_init(0, 0, &s); geographic_point_init(0, 1, &e); dist = sphere_distance(&s, &e); dir = sphere_direction(&s, &e, dist); - CU_ASSERT_DOUBLE_EQUAL(dir, 0.0, 0.0001); + /* GeodSolve -i -E -p 16 -e 1 0 --input-string "0 0 1 0" */ + CU_ASSERT_DOUBLE_EQUAL(dir, 0.0, 1e-14); + CU_ASSERT_DOUBLE_EQUAL(dist, 0.0174532925199433, 1e-14); } @@ -92,61 +102,79 @@ static void test_sphere_project(void) GEOGRAPHIC_POINT s, e; double dir1, dist1, dir2, dist2; - dir1 = M_PI/2; + dir1 = M_PI_2; dist1 = 0.1; - geographic_point_init(0, 0, &s); + geographic_point_init(0, 0, &s); sphere_project(&s, dist1, dir1, &e); + + CU_ASSERT_DOUBLE_EQUAL(e.lon, 0.1, 1e-14); + CU_ASSERT_DOUBLE_EQUAL(e.lat, 0.0, 1e-14); + /* Direct and inverse solutions agree */ dist2 = sphere_distance(&s, &e); dir2 = sphere_direction(&s, &e, dist1); - CU_ASSERT_DOUBLE_EQUAL(dist1, dist2, 0.0001); - CU_ASSERT_DOUBLE_EQUAL(dir1, dir2, 0.0001); + CU_ASSERT_DOUBLE_EQUAL(dist1, dist2, 1e-14); + CU_ASSERT_DOUBLE_EQUAL(dir1, dir2, 1e-14); dist1 = sphere_distance(&e, &s); dir1 = sphere_direction(&e, &s, dist1); sphere_project(&e, dist1, dir1, &s); - - CU_ASSERT_DOUBLE_EQUAL(s.lon, 0.0, 0.0001); - CU_ASSERT_DOUBLE_EQUAL(s.lat, 0.0, 0.0001); + + CU_ASSERT_DOUBLE_EQUAL(s.lon, 0.0, 1e-14); + CU_ASSERT_DOUBLE_EQUAL(s.lat, 0.0, 1e-14); geographic_point_init(0, 0.2, &e); geographic_point_init(0, 0.4, &s); dist1 = sphere_distance(&s, &e); dir1 = sphere_direction(&e, &s, dist1); - CU_ASSERT_DOUBLE_EQUAL(dir1, 0.0, 0.0001); + /* GeodSolve -i -E -p 16 -e 1 0 --input-string "0.2 0 0.4 0" */ + CU_ASSERT_DOUBLE_EQUAL(dir1, 0.0, 1e-14); + CU_ASSERT_DOUBLE_EQUAL(dist1, 0.0034906585039887, 1e-14); - geographic_point_init(0, 1, &s); + geographic_point_init(0, 1, &s); /* same start point for remainder of tests */ geographic_point_init(0, 2, &e); dist2 = sphere_distance(&s, &e); dir2 = sphere_direction(&s, &e, dist2); - CU_ASSERT_DOUBLE_EQUAL(dir2, 0.0, 0.0001); + /* GeodSolve -i -E -p 16 -e 1 0 --input-string "1 0 2 0" */ + CU_ASSERT_DOUBLE_EQUAL(dir2, 0.0, 1e-14); + CU_ASSERT_DOUBLE_EQUAL(dist2, 0.0174532925199433, 1e-14); geographic_point_init(1, 1, &e); dist2 = sphere_distance(&s, &e); dir2 = sphere_direction(&s, &e, dist2); - CU_ASSERT_DOUBLE_EQUAL(dir2, 1.57064, 0.0001); + /* GeodSolve -i -E -p 16 -e 1 0 --input-string "1 0 1 1" */ + CU_ASSERT_DOUBLE_EQUAL(dir2, 89.991273575329292895136 * M_PI / 180.0, 1e-14); + CU_ASSERT_DOUBLE_EQUAL(dist2, 0.0174506342314906, 1e-14); geographic_point_init(0, 0, &e); dist2 = sphere_distance(&s, &e); dir2 = sphere_direction(&s, &e, dist2); - CU_ASSERT_DOUBLE_EQUAL(dir2, 3.14159, 0.0001); + /* GeodSolve -i -E -p 16 -e 1 0 --input-string "1 0 0 0" */ + CU_ASSERT_DOUBLE_EQUAL(dir2, M_PI, 1e-14); + CU_ASSERT_DOUBLE_EQUAL(dist2, 0.0174532925199433, 1e-14); geographic_point_init(-1, 1, &e); dist2 = sphere_distance(&s, &e); dir2 = sphere_direction(&s, &e, dist2); - CU_ASSERT_DOUBLE_EQUAL(dir2, -1.57064, 0.0001); + /* GeodSolve -i -E -p 16 -e 1 0 --input-string "1 0 1 -1" */ + CU_ASSERT_DOUBLE_EQUAL(dir2, -89.991273575329292895136 * M_PI / 180.0, 1e-14); + CU_ASSERT_DOUBLE_EQUAL(dist2, 0.0174506342314906, 1e-14); geographic_point_init(1, 2, &e); dist2 = sphere_distance(&s, &e); dir2 = sphere_direction(&s, &e, dist2); - CU_ASSERT_DOUBLE_EQUAL(dir2, 0.785017, 0.0001); + /* GeodSolve -i -E -p 16 -e 1 0 --input-string "1 0 2 1" */ + CU_ASSERT_DOUBLE_EQUAL(dir2, 44.978182941465044354783 * M_PI / 180.0, 1e-14); + CU_ASSERT_DOUBLE_EQUAL(dist2, 0.0246782972905467, 1e-14); geographic_point_init(-1, 0, &e); dist2 = sphere_distance(&s, &e); dir2 = sphere_direction(&s, &e, dist2); - CU_ASSERT_DOUBLE_EQUAL(dir2, -2.35612, 0.0001); + /* GeodSolve -i -E -p 16 -e 1 0 --input-string "1 0 0 -1" */ + CU_ASSERT_DOUBLE_EQUAL(dir2, -134.995636455344851488216 * M_PI / 180.0, 1e-14); + CU_ASSERT_DOUBLE_EQUAL(dist2, 0.0246820563917664, 1e-14); } #if 0 @@ -247,12 +275,12 @@ static void test_gbox_from_spherical_coordinates(void) gbox_geocentric_slow = LW_FALSE; if ( - ( fabs( gbox.xmin - gbox_slow.xmin ) > gtolerance ) || - ( fabs( gbox.xmax - gbox_slow.xmax ) > gtolerance ) || - ( fabs( gbox.ymin - gbox_slow.ymin ) > gtolerance ) || - ( fabs( gbox.ymax - gbox_slow.ymax ) > gtolerance ) || - ( fabs( gbox.zmin - gbox_slow.zmin ) > gtolerance ) || - ( fabs( gbox.zmax - gbox_slow.zmax ) > gtolerance ) ) + ( fabs( gbox.xmin - gbox_slow.xmin ) > gtolerance ) || + ( fabs( gbox.xmax - gbox_slow.xmax ) > gtolerance ) || + ( fabs( gbox.ymin - gbox_slow.ymin ) > gtolerance ) || + ( fabs( gbox.ymax - gbox_slow.ymax ) > gtolerance ) || + ( fabs( gbox.zmin - gbox_slow.zmin ) > gtolerance ) || + ( fabs( gbox.zmax - gbox_slow.zmax ) > gtolerance ) ) { printf("\n-------\n"); printf("If you are seeing this, cut and paste it, it is a randomly generated test case!\n"); @@ -728,12 +756,12 @@ static void test_edge_distance_to_point(void) CU_ASSERT_DOUBLE_EQUAL(closest.lon, 0.0, 0.00001); /* Ticket #2351 */ - edge_set(149.386990599235, -26.3567415843982, 149.386990599247, -26.3567415843965, &e); + edge_set(149.386990599235, -26.3567415843982, 149.386990599247, -26.3567415843965, &e); point_set(149.386990599235, -26.3567415843982, &g); d = edge_distance_to_point(&e, &g, &closest); CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001); - // printf("CLOSE POINT(%g %g)\n", closest.lon, closest.lat); - // printf(" ORIG POINT(%g %g)\n", g.lon, g.lat); + // printf("CLOSE POINT(%g %g)\n", closest.lon, closest.lat); + // printf(" ORIG POINT(%g %g)\n", g.lon, g.lat); CU_ASSERT_DOUBLE_EQUAL(g.lat, closest.lat, 0.00001); CU_ASSERT_DOUBLE_EQUAL(g.lon, closest.lon, 0.00001); } @@ -1116,6 +1144,15 @@ static void test_lwpoly_covers_point2d(void) result = lwpoly_covers_point2d(poly, &pt_to_test); CU_ASSERT_EQUAL(result, LW_TRUE); lwgeom_free(lwg); + + /* Triangle over the antimeridian */ + lwg = lwgeom_from_wkt("POLYGON((140 52, 152.0 -6.0, -120.0 -29.0, 140 52))", LW_PARSER_CHECK_NONE); + poly = (LWPOLY*)lwg; + pt_to_test.x = -172.0; + pt_to_test.y = -13.0; + result = lwpoly_covers_point2d(poly, &pt_to_test); + CU_ASSERT_EQUAL(result, LW_TRUE); + lwgeom_free(lwg); } @@ -1212,7 +1249,7 @@ static void test_lwgeom_distance_sphere(void) lwgeom_free(lwg1); lwgeom_free(lwg2); - /* Ticket #2351 */ + /* Ticket #2351 */ lwg1 = lwgeom_from_wkt("LINESTRING(149.386990599235 -26.3567415843982,149.386990599247 -26.3567415843965)", LW_PARSER_CHECK_NONE); lwg2 = lwgeom_from_wkt("POINT(149.386990599235 -26.3567415843982)", LW_PARSER_CHECK_NONE); d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0); @@ -1241,40 +1278,65 @@ static void test_spheroid_distance(void) { GEOGRAPHIC_POINT g1, g2; double d; +#ifdef USE_PRE22GEODESIC + double epsilon; /* irregular */ +#else + const double epsilon = 1e-8; /* at least 10 nm precision */ +#endif SPHEROID s; /* Init to WGS84 */ spheroid_init(&s, 6378137.0, 6356752.314245179498); - /* One vertical degree */ + /* One vertical degree + $ GeodSolve -E -i -p 16 --input-string "0 0 1 0" */ point_set(0.0, 0.0, &g1); point_set(0.0, 1.0, &g2); d = spheroid_distance(&g1, &g2, &s); - CU_ASSERT_DOUBLE_EQUAL(d, 110574.388615329, 0.001); +#ifdef USE_PRE22GEODESIC + epsilon = 1e-6; +#endif + CU_ASSERT_DOUBLE_EQUAL(d, 110574.3885577987957342, epsilon); - /* Ten horizontal degrees */ + /* Ten horizontal degrees + $ GeodSolve -E -i -p 16 --input-string "0 -10 0 0" */ point_set(-10.0, 0.0, &g1); point_set(0.0, 0.0, &g2); d = spheroid_distance(&g1, &g2, &s); - CU_ASSERT_DOUBLE_EQUAL(d, 1113194.90793274, 0.001); +#ifdef USE_PRE22GEODESIC + epsilon = 1e-3; +#endif + CU_ASSERT_DOUBLE_EQUAL(d, 1113194.9079327357264771, epsilon); - /* One horizonal degree */ + /* One horizonal degree + $ GeodSolve -E -i -p 16 --input-string "0 -1 0 0" */ point_set(-1.0, 0.0, &g1); point_set(0.0, 0.0, &g2); d = spheroid_distance(&g1, &g2, &s); - CU_ASSERT_DOUBLE_EQUAL(d, 111319.490779, 0.001); +#ifdef USE_PRE22GEODESIC + epsilon = 1e-4; +#endif + CU_ASSERT_DOUBLE_EQUAL(d, 111319.4907932735726477, epsilon); - /* Around world w/ slight bend */ + /* Around world w/ slight bend + $ GeodSolve -E -i -p 16 --input-string "0 -180 1 0" */ point_set(-180.0, 0.0, &g1); point_set(0.0, 1.0, &g2); d = spheroid_distance(&g1, &g2, &s); - CU_ASSERT_DOUBLE_EQUAL(d, 19893357.0704483, 0.001); +#ifdef USE_PRE22GEODESIC + epsilon = 1e-5; +#endif + CU_ASSERT_DOUBLE_EQUAL(d, 19893357.0700676468277450, epsilon); - /* Up to pole */ + /* Up to pole + $ GeodSolve -E -i -p 16 --input-string "0 -180 90 0" */ point_set(-180.0, 0.0, &g1); point_set(0.0, 90.0, &g2); d = spheroid_distance(&g1, &g2, &s); - CU_ASSERT_DOUBLE_EQUAL(d, 10001965.7295318, 0.001); +#ifdef USE_PRE22GEODESIC + epsilon = 1e-6; +#endif + CU_ASSERT_DOUBLE_EQUAL(d, 10001965.7293127228117396, epsilon); } @@ -1293,51 +1355,60 @@ static void test_spheroid_area(void) /* Medford lot test polygon */ lwg = lwgeom_from_wkt("POLYGON((-122.848227067007 42.5007249610493,-122.848309475585 42.5007179884263,-122.848327688675 42.500835880696,-122.848245279942 42.5008428533324,-122.848227067007 42.5007249610493))", LW_PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg, &gbox); + /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE -r --input-string \ + "42.5007249610493 -122.848227067007;42.5007179884263 -122.848309475585;"\ + "42.500835880696 -122.848327688675;42.5008428533324 -122.848245279942" */ a1 = lwgeom_area_sphere(lwg, &s); + CU_ASSERT_DOUBLE_EQUAL(a1, 89.721147136698008, 0.1); + /* spheroid: Planimeter -E -p 20 -r --input-string \ + "42.5007249610493 -122.848227067007;42.5007179884263 -122.848309475585;"\ + "42.500835880696 -122.848327688675;42.5008428533324 -122.848245279942" */ a2 = lwgeom_area_spheroid(lwg, &s); - //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2); - CU_ASSERT_DOUBLE_EQUAL(a1, 89.7127703297, 0.1); /* sphere */ - CU_ASSERT_DOUBLE_EQUAL(a2, 89.8684316032, 0.1); /* spheroid */ + CU_ASSERT_DOUBLE_EQUAL(a2, 89.868413479309585, 0.1); lwgeom_free(lwg); /* Big-ass polygon */ lwg = lwgeom_from_wkt("POLYGON((-2 3, -2 4, -1 4, -1 3, -2 3))", LW_PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg, &gbox); + /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE -r --input-string "3 -2;4 -2;4 -1;3 -1" */ a1 = lwgeom_area_sphere(lwg, &s); + CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.106982993974659, 0.1); + /* spheroid: Planimeter -E -p 20 -r --input-string "3 -2;4 -2;4 -1;3 -1" */ a2 = lwgeom_area_spheroid(lwg, &s); - //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2); - CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.1, 10.0); /* sphere */ - CU_ASSERT_DOUBLE_EQUAL(a2, 12286574431.9, 10.0); /* spheroid */ + CU_ASSERT_DOUBLE_EQUAL(a2, 12286884908.946891319597874, 0.1); lwgeom_free(lwg); /* One-degree square */ lwg = lwgeom_from_wkt("POLYGON((8.5 2,8.5 1,9.5 1,9.5 2,8.5 2))", LW_PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg, &gbox); + /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE --input-string "2 8.5;1 8.5;1 9.5;2 9.5" */ a1 = lwgeom_area_sphere(lwg, &s); + CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.368023059138681, 0.1); + /* spheroid: Planimeter -E -p 20 --input-string "2 8.5;1 8.5;1 9.5;2 9.5" */ a2 = lwgeom_area_spheroid(lwg, &s); - //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2); - CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.1, 10.0); /* sphere */ - CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */ + CU_ASSERT_DOUBLE_EQUAL(a2, 12305128751.042900673161556, 0.1); lwgeom_free(lwg); - /* One-degree square *near* dateline */ + /* One-degree square *near* the antimeridian */ lwg = lwgeom_from_wkt("POLYGON((179.5 2,179.5 1,178.5 1,178.5 2,179.5 2))", LW_PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg, &gbox); + /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE -r --input-string "2 179.5;1 179.5;1 178.5;2 178.5" */ a1 = lwgeom_area_sphere(lwg, &s); + CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.368023059138681, 0.1); + /* spheroid: Planimeter -E -p 20 -r --input-string "2 179.5;1 179.5;1 178.5;2 178.5" */ a2 = lwgeom_area_spheroid(lwg, &s); - //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2); - CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.1, 10.0); /* sphere */ - CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */ + CU_ASSERT_DOUBLE_EQUAL(a2, 12305128751.042900673161556, 0.1); lwgeom_free(lwg); - /* One-degree square *across* dateline */ + /* One-degree square *across* the antimeridian */ lwg = lwgeom_from_wkt("POLYGON((179.5 2,179.5 1,-179.5 1,-179.5 2,179.5 2))", LW_PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg, &gbox); + /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE --input-string "2 179.5;1 179.5;1 -179.5;2 -179.5" */ a1 = lwgeom_area_sphere(lwg, &s); + CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.368023059138681, 0.1); + /* spheroid: Planimeter -E -p 20 --input-string "2 179.5;1 179.5;1 -179.5;2 -179.5" */ a2 = lwgeom_area_spheroid(lwg, &s); - //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2); - CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.3679, 10.0); /* sphere */ - CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */ + CU_ASSERT_DOUBLE_EQUAL(a2, 12305128751.042900673161556, 0.1); lwgeom_free(lwg); } @@ -1363,7 +1434,7 @@ static void test_gbox_utils(void) CU_ASSERT_DOUBLE_EQUAL(a2, 0.017764, 0.0000001); lwgeom_free(lwg); - /* One-degree square *across* dateline */ + /* One-degree square *across* antimeridian */ lwg = lwgeom_from_wkt("POLYGON((179.5 2,179.5 1,-179.5 1,-179.5 2,179.5 2))", LW_PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg, &gbox); a1 = gbox_angular_width(&gbox); @@ -1373,7 +1444,7 @@ static void test_gbox_utils(void) CU_ASSERT_DOUBLE_EQUAL(a2, 0.0174553, 0.0000001); lwgeom_free(lwg); - /* One-degree square *across* dateline */ + /* One-degree square *across* antimeridian */ lwg = lwgeom_from_wkt("POLYGON((178.5 2,178.5 1,-179.5 1,-179.5 2,178.5 2))", LW_PARSER_CHECK_NONE); lwgeom_calculate_gbox_geodetic(lwg, &gbox); gbox_centroid(&gbox, &pt); @@ -1395,18 +1466,18 @@ static void test_vector_angle(void) p1.x = 1.0; p2.y = 1.0; angle = vector_angle(&p1, &p2); - CU_ASSERT_DOUBLE_EQUAL(angle, M_PI/2, 0.00001); + CU_ASSERT_DOUBLE_EQUAL(angle, M_PI_2, 0.00001); p1.x = p2.y = 0.0; p1.y = 1.0; p2.x = 1.0; angle = vector_angle(&p1, &p2); - CU_ASSERT_DOUBLE_EQUAL(angle, M_PI/2, 0.00001); + CU_ASSERT_DOUBLE_EQUAL(angle, M_PI_2, 0.00001); p2.y = p2.x = 1.0; normalize(&p2); angle = vector_angle(&p1, &p2); - CU_ASSERT_DOUBLE_EQUAL(angle, M_PI/4, 0.00001); + CU_ASSERT_DOUBLE_EQUAL(angle, M_PI_4, 0.00001); p2.x = p2.y = p2.z = 1.0; normalize(&p2); @@ -1426,7 +1497,7 @@ static void test_vector_rotate(void) p1.x = 1.0; p2.y = 1.0; - angle = M_PI/4; + angle = M_PI_4; vector_rotate(&p1, &p2, angle, &n); //printf("%g %g %g\n\n", n.x, n.y, n.z); CU_ASSERT_DOUBLE_EQUAL(n.x, 0.707107, 0.00001); diff --git a/liblwgeom/cunit/cu_measures.c b/liblwgeom/cunit/cu_measures.c index c8929974e19..4bad0520b0c 100644 --- a/liblwgeom/cunit/cu_measures.c +++ b/liblwgeom/cunit/cu_measures.c @@ -733,9 +733,9 @@ test_lw_arc_length(void) /* Arc to right of the unit semicircle */ d = lw_arc_length(&A1, &A2, &A3); - CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI/2, 0.000001); + CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI_2, 0.000001); d = lw_arc_length(&A3, &A2, &A1); - CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI/2, 0.000001); + CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI_2, 0.000001); } static void @@ -765,7 +765,7 @@ test_lw_dist2d_pt_ptarrayarc(void) CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001); /* Point 45 degrees off arc, 2 radii from center */ - P.y = P.x = 2 * cos(M_PI/4); + P.y = P.x = 2 * cos(M_PI_4); lw_dist2d_distpts_init(&dl, DIST_MIN); rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl); CU_ASSERT_EQUAL( rv, LW_SUCCESS ); diff --git a/liblwgeom/g_box.c b/liblwgeom/g_box.c index 43ae3fcc8c9..2971381b86d 100644 --- a/liblwgeom/g_box.c +++ b/liblwgeom/g_box.c @@ -551,7 +551,7 @@ static int lwcircstring_calculate_gbox_cartesian(LWCIRCSTRING *curve, GBOX *gbox /* Initialize */ gbox->xmin = gbox->ymin = gbox->zmin = gbox->mmin = FLT_MAX; - gbox->xmax = gbox->ymax = gbox->zmax = gbox->mmax = -1 * FLT_MAX; + gbox->xmax = gbox->ymax = gbox->zmax = gbox->mmax = FLT_MIN; for ( i = 2; i < curve->points->npoints; i += 2 ) { diff --git a/liblwgeom/geodesic.c b/liblwgeom/geodesic.c new file mode 100644 index 00000000000..fd0214c7d36 --- /dev/null +++ b/liblwgeom/geodesic.c @@ -0,0 +1,1781 @@ +/** + * \file geodesic.c + * \brief Implementation of the geodesic routines in C + * + * For the full documentation see geodesic.h. + **********************************************************************/ + +/** @cond SKIP */ + +/* + * This is a C implementation of the geodesic algorithms described in + * + * C. F. F. Karney, + * Algorithms for geodesics, + * J. Geodesy 87, 43--55 (2013); + * https://dx.doi.org/10.1007/s00190-012-0578-z + * Addenda: http://geographiclib.sf.net/geod-addenda.html + * + * See the comments in geodesic.h for documentation. + * + * Copyright (c) Charles Karney (2012-2014) and licensed + * under the MIT/X11 License. For more information, see + * http://geographiclib.sourceforge.net/ + */ + +#include "geodesic.h" +#include + +#define GEOGRAPHICLIB_GEODESIC_ORDER 6 +#define nC1 GEOGRAPHICLIB_GEODESIC_ORDER +#define nC1p GEOGRAPHICLIB_GEODESIC_ORDER +#define nC2 GEOGRAPHICLIB_GEODESIC_ORDER +#define nA3 GEOGRAPHICLIB_GEODESIC_ORDER +#define nA3x nA3 +#define nC3 GEOGRAPHICLIB_GEODESIC_ORDER +#define nC3x ((nC3 * (nC3 - 1)) / 2) +#define nC4 GEOGRAPHICLIB_GEODESIC_ORDER +#define nC4x ((nC4 * (nC4 + 1)) / 2) + +typedef double real; +typedef int boolx; + +static unsigned init = 0; +static const int FALSE = 0; +static const int TRUE = 1; +static unsigned digits, maxit1, maxit2; +static real epsilon, realmin, pi, degree, NaN, + tiny, tol0, tol1, tol2, tolb, xthresh; + +static void Init() { + if (!init) { +#if defined(__DBL_MANT_DIG__) + digits = __DBL_MANT_DIG__; +#else + digits = 53; +#endif +#if defined(__DBL_EPSILON__) + epsilon = __DBL_EPSILON__; +#else + epsilon = pow(0.5, digits - 1); +#endif +#if defined(__DBL_MIN__) + realmin = __DBL_MIN__; +#else + realmin = pow(0.5, 1022); +#endif +#if defined(M_PI) + pi = M_PI; +#else + pi = atan2(0.0, -1.0); +#endif + maxit1 = 20; + maxit2 = maxit1 + digits + 10; + tiny = sqrt(realmin); + tol0 = epsilon; + /* Increase multiplier in defn of tol1 from 100 to 200 to fix inverse case + * 52.784459512564 0 -52.784459512563990912 179.634407464943777557 + * which otherwise failed for Visual Studio 10 (Release and Debug) */ + tol1 = 200 * tol0; + tol2 = sqrt(tol0); + /* Check on bisection interval */ + tolb = tol0 * tol2; + xthresh = 1000 * tol2; + degree = pi/180; + NaN = sqrt(-1.0); + init = 1; + } +} + +enum captype { + CAP_NONE = 0U, + CAP_C1 = 1U<<0, + CAP_C1p = 1U<<1, + CAP_C2 = 1U<<2, + CAP_C3 = 1U<<3, + CAP_C4 = 1U<<4, + CAP_ALL = 0x1FU, + OUT_ALL = 0x7F80U +}; + +static real sq(real x) { return x * x; } +static real log1px(real x) { + volatile real + y = 1 + x, + z = y - 1; + /* Here's the explanation for this magic: y = 1 + z, exactly, and z + * approx x, thus log(y)/z (which is nearly constant near z = 0) returns + * a good approximation to the true log(1 + x)/x. The multiplication x * + * (log(y)/z) introduces little additional error. */ + return z == 0 ? x : x * log(y) / z; +} + +static real atanhx(real x) { + real y = fabs(x); /* Enforce odd parity */ + y = log1px(2 * y/(1 - y))/2; + return x < 0 ? -y : y; +} + +static real hypotx(real x, real y) +{ return sqrt(x * x + y * y); } + +static real cbrtx(real x) { + real y = pow(fabs(x), 1/(real)(3)); /* Return the real cube root */ + return x < 0 ? -y : y; +} + +static real sumx(real u, real v, real* t) { + volatile real s = u + v; + volatile real up = s - v; + volatile real vpp = s - up; + up -= u; + vpp -= v; + *t = -(up + vpp); + /* error-free sum: + * u + v = s + t + * = round(u + v) + t */ + return s; +} + +static real minx(real x, real y) +{ return x < y ? x : y; } + +static real maxx(real x, real y) +{ return x > y ? x : y; } + +static void swapx(real* x, real* y) +{ real t = *x; *x = *y; *y = t; } + +static void SinCosNorm(real* sinx, real* cosx) { + real r = hypotx(*sinx, *cosx); + *sinx /= r; + *cosx /= r; +} + +static real AngNormalize(real x) +{ return x >= 180 ? x - 360 : (x < -180 ? x + 360 : x); } +static real AngNormalize2(real x) +{ return AngNormalize(fmod(x, (real)(360))); } + +static real AngDiff(real x, real y) { + real t, d = sumx(-x, y, &t); + if ((d - (real)(180)) + t > (real)(0)) /* y - x > 180 */ + d -= (real)(360); /* exact */ + else if ((d + (real)(180)) + t <= (real)(0)) /* y - x <= -180 */ + d += (real)(360); /* exact */ + return d + t; +} + +static real AngRound(real x) { + const real z = 1/(real)(16); + volatile real y = fabs(x); + /* The compiler mustn't "simplify" z - (z - y) to y */ + y = y < z ? z - (z - y) : y; + return x < 0 ? -y : y; +} + +static void A3coeff(struct geod_geodesic* g); +static void C3coeff(struct geod_geodesic* g); +static void C4coeff(struct geod_geodesic* g); +static real SinCosSeries(boolx sinp, + real sinx, real cosx, + const real c[], int n); +static void Lengths(const struct geod_geodesic* g, + real eps, real sig12, + real ssig1, real csig1, real dn1, + real ssig2, real csig2, real dn2, + real cbet1, real cbet2, + real* ps12b, real* pm12b, real* pm0, + boolx scalep, real* pM12, real* pM21, + /* Scratch areas of the right size */ + real C1a[], real C2a[]); +static real Astroid(real x, real y); +static real InverseStart(const struct geod_geodesic* g, + real sbet1, real cbet1, real dn1, + real sbet2, real cbet2, real dn2, + real lam12, + real* psalp1, real* pcalp1, + /* Only updated if return val >= 0 */ + real* psalp2, real* pcalp2, + /* Only updated for short lines */ + real* pdnm, + /* Scratch areas of the right size */ + real C1a[], real C2a[]); +static real Lambda12(const struct geod_geodesic* g, + real sbet1, real cbet1, real dn1, + real sbet2, real cbet2, real dn2, + real salp1, real calp1, + real* psalp2, real* pcalp2, + real* psig12, + real* pssig1, real* pcsig1, + real* pssig2, real* pcsig2, + real* peps, real* pdomg12, + boolx diffp, real* pdlam12, + /* Scratch areas of the right size */ + real C1a[], real C2a[], real C3a[]); +static real A3f(const struct geod_geodesic* g, real eps); +static void C3f(const struct geod_geodesic* g, real eps, real c[]); +static void C4f(const struct geod_geodesic* g, real eps, real c[]); +static real A1m1f(real eps); +static void C1f(real eps, real c[]); +static void C1pf(real eps, real c[]); +static real A2m1f(real eps); +static void C2f(real eps, real c[]); +static int transit(real lon1, real lon2); +static int transitdirect(real lon1, real lon2); +static void accini(real s[]); +static void acccopy(const real s[], real t[]); +static void accadd(real s[], real y); +static real accsum(const real s[], real y); +static void accneg(real s[]); + +void geod_init(struct geod_geodesic* g, real a, real f) { + if (!init) Init(); + g->a = a; + g->f = f <= 1 ? f : 1/f; + g->f1 = 1 - g->f; + g->e2 = g->f * (2 - g->f); + g->ep2 = g->e2 / sq(g->f1); /* e2 / (1 - e2) */ + g->n = g->f / ( 2 - g->f); + g->b = g->a * g->f1; + g->c2 = (sq(g->a) + sq(g->b) * + (g->e2 == 0 ? 1 : + (g->e2 > 0 ? atanhx(sqrt(g->e2)) : atan(sqrt(-g->e2))) / + sqrt(fabs(g->e2))))/2; /* authalic radius squared */ + /* The sig12 threshold for "really short". Using the auxiliary sphere + * solution with dnm computed at (bet1 + bet2) / 2, the relative error in the + * azimuth consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2. (Error + * measured for 1/100 < b/a < 100 and abs(f) >= 1/1000. For a given f and + * sig12, the max error occurs for lines near the pole. If the old rule for + * computing dnm = (dn1 + dn2)/2 is used, then the error increases by a + * factor of 2.) Setting this equal to epsilon gives sig12 = etol2. Here + * 0.1 is a safety factor (error decreased by 100) and max(0.001, abs(f)) + * stops etol2 getting too large in the nearly spherical case. */ + g->etol2 = 0.1 * tol2 / + sqrt( maxx((real)(0.001), fabs(g->f)) * minx((real)(1), 1 - g->f/2) / 2 ); + + A3coeff(g); + C3coeff(g); + C4coeff(g); +} + +void geod_lineinit(struct geod_geodesicline* l, + const struct geod_geodesic* g, + real lat1, real lon1, real azi1, unsigned caps) { + real alp1, cbet1, sbet1, phi, eps; + l->a = g->a; + l->f = g->f; + l->b = g->b; + l->c2 = g->c2; + l->f1 = g->f1; + /* If caps is 0 assume the standard direct calculation */ + l->caps = (caps ? caps : GEOD_DISTANCE_IN | GEOD_LONGITUDE) | + GEOD_LATITUDE | GEOD_AZIMUTH; /* Always allow latitude and azimuth */ + + l->lat1 = lat1; + l->lon1 = lon1; + /* Guard against underflow in salp0 */ + l->azi1 = AngRound(AngNormalize(azi1)); + /* alp1 is in [0, pi] */ + alp1 = l->azi1 * degree; + /* Enforce sin(pi) == 0 and cos(pi/2) == 0. Better to face the ensuing + * problems directly than to skirt them. */ + l->salp1 = l->azi1 == -180 ? 0 : sin(alp1); + l->calp1 = fabs(l->azi1) == 90 ? 0 : cos(alp1); + phi = lat1 * degree; + /* Ensure cbet1 = +epsilon at poles */ + sbet1 = l->f1 * sin(phi); + cbet1 = fabs(lat1) == 90 ? tiny : cos(phi); + SinCosNorm(&sbet1, &cbet1); + l->dn1 = sqrt(1 + g->ep2 * sq(sbet1)); + + /* Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0), */ + l->salp0 = l->salp1 * cbet1; /* alp0 in [0, pi/2 - |bet1|] */ + /* Alt: calp0 = hypot(sbet1, calp1 * cbet1). The following + * is slightly better (consider the case salp1 = 0). */ + l->calp0 = hypotx(l->calp1, l->salp1 * sbet1); + /* Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1). + * sig = 0 is nearest northward crossing of equator. + * With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line). + * With bet1 = pi/2, alp1 = -pi, sig1 = pi/2 + * With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2 + * Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1). + * With alp0 in (0, pi/2], quadrants for sig and omg coincide. + * No atan2(0,0) ambiguity at poles since cbet1 = +epsilon. + * With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi. */ + l->ssig1 = sbet1; l->somg1 = l->salp0 * sbet1; + l->csig1 = l->comg1 = sbet1 != 0 || l->calp1 != 0 ? cbet1 * l->calp1 : 1; + SinCosNorm(&l->ssig1, &l->csig1); /* sig1 in (-pi, pi] */ + /* SinCosNorm(somg1, comg1); -- don't need to normalize! */ + + l->k2 = sq(l->calp0) * g->ep2; + eps = l->k2 / (2 * (1 + sqrt(1 + l->k2)) + l->k2); + + if (l->caps & CAP_C1) { + real s, c; + l->A1m1 = A1m1f(eps); + C1f(eps, l->C1a); + l->B11 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C1a, nC1); + s = sin(l->B11); c = cos(l->B11); + /* tau1 = sig1 + B11 */ + l->stau1 = l->ssig1 * c + l->csig1 * s; + l->ctau1 = l->csig1 * c - l->ssig1 * s; + /* Not necessary because C1pa reverts C1a + * B11 = -SinCosSeries(TRUE, stau1, ctau1, C1pa, nC1p); */ + } + + if (l->caps & CAP_C1p) + C1pf(eps, l->C1pa); + + if (l->caps & CAP_C2) { + l->A2m1 = A2m1f(eps); + C2f(eps, l->C2a); + l->B21 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C2a, nC2); + } + + if (l->caps & CAP_C3) { + C3f(g, eps, l->C3a); + l->A3c = -l->f * l->salp0 * A3f(g, eps); + l->B31 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C3a, nC3-1); + } + + if (l->caps & CAP_C4) { + C4f(g, eps, l->C4a); + /* Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0) */ + l->A4 = sq(l->a) * l->calp0 * l->salp0 * g->e2; + l->B41 = SinCosSeries(FALSE, l->ssig1, l->csig1, l->C4a, nC4); + } +} + +real geod_genposition(const struct geod_geodesicline* l, + unsigned flags, real s12_a12, + real* plat2, real* plon2, real* pazi2, + real* ps12, real* pm12, + real* pM12, real* pM21, + real* pS12) { + real lat2 = 0, lon2 = 0, azi2 = 0, s12 = 0, + m12 = 0, M12 = 0, M21 = 0, S12 = 0; + /* Avoid warning about uninitialized B12. */ + real sig12, ssig12, csig12, B12 = 0, AB1 = 0; + real omg12, lam12, lon12; + real ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2, dn2; + unsigned outmask = + (plat2 ? GEOD_LATITUDE : 0U) | + (plon2 ? GEOD_LONGITUDE : 0U) | + (pazi2 ? GEOD_AZIMUTH : 0U) | + (ps12 ? GEOD_DISTANCE : 0U) | + (pm12 ? GEOD_REDUCEDLENGTH : 0U) | + (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) | + (pS12 ? GEOD_AREA : 0U); + + outmask &= l->caps & OUT_ALL; + if (!( TRUE /*Init()*/ && + (flags & GEOD_ARCMODE || (l->caps & GEOD_DISTANCE_IN & OUT_ALL)) )) + /* Uninitialized or impossible distance calculation requested */ + return NaN; + + if (flags & GEOD_ARCMODE) { + real s12a; + /* Interpret s12_a12 as spherical arc length */ + sig12 = s12_a12 * degree; + s12a = fabs(s12_a12); + s12a -= 180 * floor(s12a / 180); + ssig12 = s12a == 0 ? 0 : sin(sig12); + csig12 = s12a == 90 ? 0 : cos(sig12); + } else { + /* Interpret s12_a12 as distance */ + real + tau12 = s12_a12 / (l->b * (1 + l->A1m1)), + s = sin(tau12), + c = cos(tau12); + /* tau2 = tau1 + tau12 */ + B12 = - SinCosSeries(TRUE, + l->stau1 * c + l->ctau1 * s, + l->ctau1 * c - l->stau1 * s, + l->C1pa, nC1p); + sig12 = tau12 - (B12 - l->B11); + ssig12 = sin(sig12); csig12 = cos(sig12); + if (fabs(l->f) > 0.01) { + /* Reverted distance series is inaccurate for |f| > 1/100, so correct + * sig12 with 1 Newton iteration. The following table shows the + * approximate maximum error for a = WGS_a() and various f relative to + * GeodesicExact. + * erri = the error in the inverse solution (nm) + * errd = the error in the direct solution (series only) (nm) + * errda = the error in the direct solution (series + 1 Newton) (nm) + * + * f erri errd errda + * -1/5 12e6 1.2e9 69e6 + * -1/10 123e3 12e6 765e3 + * -1/20 1110 108e3 7155 + * -1/50 18.63 200.9 27.12 + * -1/100 18.63 23.78 23.37 + * -1/150 18.63 21.05 20.26 + * 1/150 22.35 24.73 25.83 + * 1/100 22.35 25.03 25.31 + * 1/50 29.80 231.9 30.44 + * 1/20 5376 146e3 10e3 + * 1/10 829e3 22e6 1.5e6 + * 1/5 157e6 3.8e9 280e6 */ + real + ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12, + csig2 = l->csig1 * csig12 - l->ssig1 * ssig12, + serr; + B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1); + serr = (1 + l->A1m1) * (sig12 + (B12 - l->B11)) - s12_a12 / l->b; + sig12 = sig12 - serr / sqrt(1 + l->k2 * sq(ssig2)); + ssig12 = sin(sig12); csig12 = cos(sig12); + /* Update B12 below */ + } + } + + /* sig2 = sig1 + sig12 */ + ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12; + csig2 = l->csig1 * csig12 - l->ssig1 * ssig12; + dn2 = sqrt(1 + l->k2 * sq(ssig2)); + if (outmask & (GEOD_DISTANCE | GEOD_REDUCEDLENGTH | GEOD_GEODESICSCALE)) { + if (flags & GEOD_ARCMODE || fabs(l->f) > 0.01) + B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1); + AB1 = (1 + l->A1m1) * (B12 - l->B11); + } + /* sin(bet2) = cos(alp0) * sin(sig2) */ + sbet2 = l->calp0 * ssig2; + /* Alt: cbet2 = hypot(csig2, salp0 * ssig2); */ + cbet2 = hypotx(l->salp0, l->calp0 * csig2); + if (cbet2 == 0) + /* I.e., salp0 = 0, csig2 = 0. Break the degeneracy in this case */ + cbet2 = csig2 = tiny; + /* tan(alp0) = cos(sig2)*tan(alp2) */ + salp2 = l->salp0; calp2 = l->calp0 * csig2; /* No need to normalize */ + + if (outmask & GEOD_DISTANCE) + s12 = flags & GEOD_ARCMODE ? l->b * ((1 + l->A1m1) * sig12 + AB1) : s12_a12; + + if (outmask & GEOD_LONGITUDE) { + /* tan(omg2) = sin(alp0) * tan(sig2) */ + somg2 = l->salp0 * ssig2; comg2 = csig2; /* No need to normalize */ + /* omg12 = omg2 - omg1 */ + omg12 = flags & GEOD_LONG_NOWRAP ? sig12 + - (atan2(ssig2, csig2) - atan2(l->ssig1, l->csig1)) + + (atan2(somg2, comg2) - atan2(l->somg1, l->comg1)) + : atan2(somg2 * l->comg1 - comg2 * l->somg1, + comg2 * l->comg1 + somg2 * l->somg1); + lam12 = omg12 + l->A3c * + ( sig12 + (SinCosSeries(TRUE, ssig2, csig2, l->C3a, nC3-1) + - l->B31)); + lon12 = lam12 / degree; + /* Use AngNormalize2 because longitude might have wrapped multiple + * times. */ + lon2 = flags & GEOD_LONG_NOWRAP ? l->lon1 + lon12 : + AngNormalize(AngNormalize(l->lon1) + AngNormalize2(lon12)); + } + + if (outmask & GEOD_LATITUDE) + lat2 = atan2(sbet2, l->f1 * cbet2) / degree; + + if (outmask & GEOD_AZIMUTH) + /* minus signs give range [-180, 180). 0- converts -0 to +0. */ + azi2 = 0 - atan2(-salp2, calp2) / degree; + + if (outmask & (GEOD_REDUCEDLENGTH | GEOD_GEODESICSCALE)) { + real + B22 = SinCosSeries(TRUE, ssig2, csig2, l->C2a, nC2), + AB2 = (1 + l->A2m1) * (B22 - l->B21), + J12 = (l->A1m1 - l->A2m1) * sig12 + (AB1 - AB2); + if (outmask & GEOD_REDUCEDLENGTH) + /* Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure + * accurate cancellation in the case of coincident points. */ + m12 = l->b * ((dn2 * (l->csig1 * ssig2) - l->dn1 * (l->ssig1 * csig2)) + - l->csig1 * csig2 * J12); + if (outmask & GEOD_GEODESICSCALE) { + real t = l->k2 * (ssig2 - l->ssig1) * (ssig2 + l->ssig1) / (l->dn1 + dn2); + M12 = csig12 + (t * ssig2 - csig2 * J12) * l->ssig1 / l->dn1; + M21 = csig12 - (t * l->ssig1 - l->csig1 * J12) * ssig2 / dn2; + } + } + + if (outmask & GEOD_AREA) { + real + B42 = SinCosSeries(FALSE, ssig2, csig2, l->C4a, nC4); + real salp12, calp12; + if (l->calp0 == 0 || l->salp0 == 0) { + /* alp12 = alp2 - alp1, used in atan2 so no need to normalized */ + salp12 = salp2 * l->calp1 - calp2 * l->salp1; + calp12 = calp2 * l->calp1 + salp2 * l->salp1; + /* The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz + * salp12 = -0 and alp12 = -180. However this depends on the sign being + * attached to 0 correctly. The following ensures the correct + * behavior. */ + if (salp12 == 0 && calp12 < 0) { + salp12 = tiny * l->calp1; + calp12 = -1; + } + } else { + /* tan(alp) = tan(alp0) * sec(sig) + * tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1) + * = calp0 * salp0 * (csig1-csig2) / (salp0^2 + calp0^2 * csig1*csig2) + * If csig12 > 0, write + * csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1) + * else + * csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1 + * No need to normalize */ + salp12 = l->calp0 * l->salp0 * + (csig12 <= 0 ? l->csig1 * (1 - csig12) + ssig12 * l->ssig1 : + ssig12 * (l->csig1 * ssig12 / (1 + csig12) + l->ssig1)); + calp12 = sq(l->salp0) + sq(l->calp0) * l->csig1 * csig2; + } + S12 = l->c2 * atan2(salp12, calp12) + l->A4 * (B42 - l->B41); + } + + if (outmask & GEOD_LATITUDE) + *plat2 = lat2; + if (outmask & GEOD_LONGITUDE) + *plon2 = lon2; + if (outmask & GEOD_AZIMUTH) + *pazi2 = azi2; + if (outmask & GEOD_DISTANCE) + *ps12 = s12; + if (outmask & GEOD_REDUCEDLENGTH) + *pm12 = m12; + if (outmask & GEOD_GEODESICSCALE) { + if (pM12) *pM12 = M12; + if (pM21) *pM21 = M21; + } + if (outmask & GEOD_AREA) + *pS12 = S12; + + return flags & GEOD_ARCMODE ? s12_a12 : sig12 / degree; +} + +void geod_position(const struct geod_geodesicline* l, real s12, + real* plat2, real* plon2, real* pazi2) { + geod_genposition(l, FALSE, s12, plat2, plon2, pazi2, 0, 0, 0, 0, 0); +} + +real geod_gendirect(const struct geod_geodesic* g, + real lat1, real lon1, real azi1, + unsigned flags, real s12_a12, + real* plat2, real* plon2, real* pazi2, + real* ps12, real* pm12, real* pM12, real* pM21, + real* pS12) { + struct geod_geodesicline l; + unsigned outmask = + (plat2 ? GEOD_LATITUDE : 0U) | + (plon2 ? GEOD_LONGITUDE : 0U) | + (pazi2 ? GEOD_AZIMUTH : 0U) | + (ps12 ? GEOD_DISTANCE : 0U) | + (pm12 ? GEOD_REDUCEDLENGTH : 0U) | + (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) | + (pS12 ? GEOD_AREA : 0U); + + geod_lineinit(&l, g, lat1, lon1, azi1, + /* Automatically supply GEOD_DISTANCE_IN if necessary */ + outmask | + (flags & GEOD_ARCMODE ? GEOD_NONE : GEOD_DISTANCE_IN)); + return geod_genposition(&l, flags, s12_a12, + plat2, plon2, pazi2, ps12, pm12, pM12, pM21, pS12); +} + +void geod_direct(const struct geod_geodesic* g, + real lat1, real lon1, real azi1, + real s12, + real* plat2, real* plon2, real* pazi2) { + geod_gendirect(g, lat1, lon1, azi1, GEOD_NOFLAGS, s12, plat2, plon2, pazi2, + 0, 0, 0, 0, 0); +} + +real geod_geninverse(const struct geod_geodesic* g, + real lat1, real lon1, real lat2, real lon2, + real* ps12, real* pazi1, real* pazi2, + real* pm12, real* pM12, real* pM21, real* pS12) { + real s12 = 0, azi1 = 0, azi2 = 0, m12 = 0, M12 = 0, M21 = 0, S12 = 0; + real lon12; + int latsign, lonsign, swapp; + real phi, sbet1, cbet1, sbet2, cbet2, s12x = 0, m12x = 0; + real dn1, dn2, lam12, slam12, clam12; + real a12 = 0, sig12, calp1 = 0, salp1 = 0, calp2 = 0, salp2 = 0; + /* index zero elements of these arrays are unused */ + real C1a[nC1 + 1], C2a[nC2 + 1], C3a[nC3]; + boolx meridian; + real omg12 = 0; + + unsigned outmask = + (ps12 ? GEOD_DISTANCE : 0U) | + (pazi1 || pazi2 ? GEOD_AZIMUTH : 0U) | + (pm12 ? GEOD_REDUCEDLENGTH : 0U) | + (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) | + (pS12 ? GEOD_AREA : 0U); + + outmask &= OUT_ALL; + /* Compute longitude difference (AngDiff does this carefully). Result is + * in [-180, 180] but -180 is only for west-going geodesics. 180 is for + * east-going and meridional geodesics. */ + lon12 = AngDiff(AngNormalize(lon1), AngNormalize(lon2)); + /* If very close to being on the same half-meridian, then make it so. */ + lon12 = AngRound(lon12); + /* Make longitude difference positive. */ + lonsign = lon12 >= 0 ? 1 : -1; + lon12 *= lonsign; + /* If really close to the equator, treat as on equator. */ + lat1 = AngRound(lat1); + lat2 = AngRound(lat2); + /* Swap points so that point with higher (abs) latitude is point 1 */ + swapp = fabs(lat1) >= fabs(lat2) ? 1 : -1; + if (swapp < 0) { + lonsign *= -1; + swapx(&lat1, &lat2); + } + /* Make lat1 <= 0 */ + latsign = lat1 < 0 ? 1 : -1; + lat1 *= latsign; + lat2 *= latsign; + /* Now we have + * + * 0 <= lon12 <= 180 + * -90 <= lat1 <= 0 + * lat1 <= lat2 <= -lat1 + * + * longsign, swapp, latsign register the transformation to bring the + * coordinates to this canonical form. In all cases, 1 means no change was + * made. We make these transformations so that there are few cases to + * check, e.g., on verifying quadrants in atan2. In addition, this + * enforces some symmetries in the results returned. */ + + phi = lat1 * degree; + /* Ensure cbet1 = +epsilon at poles */ + sbet1 = g->f1 * sin(phi); + cbet1 = lat1 == -90 ? tiny : cos(phi); + SinCosNorm(&sbet1, &cbet1); + + phi = lat2 * degree; + /* Ensure cbet2 = +epsilon at poles */ + sbet2 = g->f1 * sin(phi); + cbet2 = fabs(lat2) == 90 ? tiny : cos(phi); + SinCosNorm(&sbet2, &cbet2); + + /* If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the + * |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1 is + * a better measure. This logic is used in assigning calp2 in Lambda12. + * Sometimes these quantities vanish and in that case we force bet2 = +/- + * bet1 exactly. An example where is is necessary is the inverse problem + * 48.522876735459 0 -48.52287673545898293 179.599720456223079643 + * which failed with Visual Studio 10 (Release and Debug) */ + + if (cbet1 < -sbet1) { + if (cbet2 == cbet1) + sbet2 = sbet2 < 0 ? sbet1 : -sbet1; + } else { + if (fabs(sbet2) == -sbet1) + cbet2 = cbet1; + } + + dn1 = sqrt(1 + g->ep2 * sq(sbet1)); + dn2 = sqrt(1 + g->ep2 * sq(sbet2)); + + lam12 = lon12 * degree; + slam12 = lon12 == 180 ? 0 : sin(lam12); + clam12 = cos(lam12); /* lon12 == 90 isn't interesting */ + + meridian = lat1 == -90 || slam12 == 0; + + if (meridian) { + + /* Endpoints are on a single full meridian, so the geodesic might lie on + * a meridian. */ + + real ssig1, csig1, ssig2, csig2; + calp1 = clam12; salp1 = slam12; /* Head to the target longitude */ + calp2 = 1; salp2 = 0; /* At the target we're heading north */ + + /* tan(bet) = tan(sig) * cos(alp) */ + ssig1 = sbet1; csig1 = calp1 * cbet1; + ssig2 = sbet2; csig2 = calp2 * cbet2; + + /* sig12 = sig2 - sig1 */ + sig12 = atan2(maxx(csig1 * ssig2 - ssig1 * csig2, (real)(0)), + csig1 * csig2 + ssig1 * ssig2); + { + real dummy; + Lengths(g, g->n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, + cbet1, cbet2, &s12x, &m12x, &dummy, + (outmask & GEOD_GEODESICSCALE) != 0U, &M12, &M21, C1a, C2a); + } + /* Add the check for sig12 since zero length geodesics might yield m12 < + * 0. Test case was + * + * echo 20.001 0 20.001 0 | GeodSolve -i + * + * In fact, we will have sig12 > pi/2 for meridional geodesic which is + * not a shortest path. */ + if (sig12 < 1 || m12x >= 0) { + m12x *= g->b; + s12x *= g->b; + a12 = sig12 / degree; + } else + /* m12 < 0, i.e., prolate and too close to anti-podal */ + meridian = FALSE; + } + + if (!meridian && + sbet1 == 0 && /* and sbet2 == 0 */ + /* Mimic the way Lambda12 works with calp1 = 0 */ + (g->f <= 0 || lam12 <= pi - g->f * pi)) { + + /* Geodesic runs along equator */ + calp1 = calp2 = 0; salp1 = salp2 = 1; + s12x = g->a * lam12; + sig12 = omg12 = lam12 / g->f1; + m12x = g->b * sin(sig12); + if (outmask & GEOD_GEODESICSCALE) + M12 = M21 = cos(sig12); + a12 = lon12 / g->f1; + + } else if (!meridian) { + + /* Now point1 and point2 belong within a hemisphere bounded by a + * meridian and geodesic is neither meridional or equatorial. */ + + /* Figure a starting point for Newton's method */ + real dnm = 0; + sig12 = InverseStart(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, + lam12, + &salp1, &calp1, &salp2, &calp2, &dnm, + C1a, C2a); + + if (sig12 >= 0) { + /* Short lines (InverseStart sets salp2, calp2, dnm) */ + s12x = sig12 * g->b * dnm; + m12x = sq(dnm) * g->b * sin(sig12 / dnm); + if (outmask & GEOD_GEODESICSCALE) + M12 = M21 = cos(sig12 / dnm); + a12 = sig12 / degree; + omg12 = lam12 / (g->f1 * dnm); + } else { + + /* Newton's method. This is a straightforward solution of f(alp1) = + * lambda12(alp1) - lam12 = 0 with one wrinkle. f(alp) has exactly one + * root in the interval (0, pi) and its derivative is positive at the + * root. Thus f(alp) is positive for alp > alp1 and negative for alp < + * alp1. During the course of the iteration, a range (alp1a, alp1b) is + * maintained which brackets the root and with each evaluation of + * f(alp) the range is shrunk, if possible. Newton's method is + * restarted whenever the derivative of f is negative (because the new + * value of alp1 is then further from the solution) or if the new + * estimate of alp1 lies outside (0,pi); in this case, the new starting + * guess is taken to be (alp1a + alp1b) / 2. */ + real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0; + unsigned numit = 0; + /* Bracketing range */ + real salp1a = tiny, calp1a = 1, salp1b = tiny, calp1b = -1; + boolx tripn, tripb; + for (tripn = FALSE, tripb = FALSE; numit < maxit2; ++numit) { + /* the WGS84 test set: mean = 1.47, sd = 1.25, max = 16 + * WGS84 and random input: mean = 2.85, sd = 0.60 */ + real dv, + v = (Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1, + &salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2, + &eps, &omg12, numit < maxit1, &dv, C1a, C2a, C3a) + - lam12); + /* 2 * tol0 is approximately 1 ulp for a number in [0, pi]. */ + /* Reversed test to allow escape with NaNs */ + if (tripb || !(fabs(v) >= (tripn ? 8 : 2) * tol0)) break; + /* Update bracketing values */ + if (v > 0 && (numit > maxit1 || calp1/salp1 > calp1b/salp1b)) + { salp1b = salp1; calp1b = calp1; } + else if (v < 0 && (numit > maxit1 || calp1/salp1 < calp1a/salp1a)) + { salp1a = salp1; calp1a = calp1; } + if (numit < maxit1 && dv > 0) { + real + dalp1 = -v/dv; + real + sdalp1 = sin(dalp1), cdalp1 = cos(dalp1), + nsalp1 = salp1 * cdalp1 + calp1 * sdalp1; + if (nsalp1 > 0 && fabs(dalp1) < pi) { + calp1 = calp1 * cdalp1 - salp1 * sdalp1; + salp1 = nsalp1; + SinCosNorm(&salp1, &calp1); + /* In some regimes we don't get quadratic convergence because + * slope -> 0. So use convergence conditions based on epsilon + * instead of sqrt(epsilon). */ + tripn = fabs(v) <= 16 * tol0; + continue; + } + } + /* Either dv was not postive or updated value was outside legal + * range. Use the midpoint of the bracket as the next estimate. + * This mechanism is not needed for the WGS84 ellipsoid, but it does + * catch problems with more eccentric ellipsoids. Its efficacy is + * such for the WGS84 test set with the starting guess set to alp1 = + * 90deg: + * the WGS84 test set: mean = 5.21, sd = 3.93, max = 24 + * WGS84 and random input: mean = 4.74, sd = 0.99 */ + salp1 = (salp1a + salp1b)/2; + calp1 = (calp1a + calp1b)/2; + SinCosNorm(&salp1, &calp1); + tripn = FALSE; + tripb = (fabs(salp1a - salp1) + (calp1a - calp1) < tolb || + fabs(salp1 - salp1b) + (calp1 - calp1b) < tolb); + } + { + real dummy; + Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, + cbet1, cbet2, &s12x, &m12x, &dummy, + (outmask & GEOD_GEODESICSCALE) != 0U, &M12, &M21, C1a, C2a); + } + m12x *= g->b; + s12x *= g->b; + a12 = sig12 / degree; + omg12 = lam12 - omg12; + } + } + + if (outmask & GEOD_DISTANCE) + s12 = 0 + s12x; /* Convert -0 to 0 */ + + if (outmask & GEOD_REDUCEDLENGTH) + m12 = 0 + m12x; /* Convert -0 to 0 */ + + if (outmask & GEOD_AREA) { + real + /* From Lambda12: sin(alp1) * cos(bet1) = sin(alp0) */ + salp0 = salp1 * cbet1, + calp0 = hypotx(calp1, salp1 * sbet1); /* calp0 > 0 */ + real alp12; + if (calp0 != 0 && salp0 != 0) { + real + /* From Lambda12: tan(bet) = tan(sig) * cos(alp) */ + ssig1 = sbet1, csig1 = calp1 * cbet1, + ssig2 = sbet2, csig2 = calp2 * cbet2, + k2 = sq(calp0) * g->ep2, + eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2), + /* Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0). */ + A4 = sq(g->a) * calp0 * salp0 * g->e2; + real C4a[nC4]; + real B41, B42; + SinCosNorm(&ssig1, &csig1); + SinCosNorm(&ssig2, &csig2); + C4f(g, eps, C4a); + B41 = SinCosSeries(FALSE, ssig1, csig1, C4a, nC4); + B42 = SinCosSeries(FALSE, ssig2, csig2, C4a, nC4); + S12 = A4 * (B42 - B41); + } else + /* Avoid problems with indeterminate sig1, sig2 on equator */ + S12 = 0; + + if (!meridian && + omg12 < (real)(0.75) * pi && /* Long difference too big */ + sbet2 - sbet1 < (real)(1.75)) { /* Lat difference too big */ + /* Use tan(Gamma/2) = tan(omg12/2) + * * (tan(bet1/2)+tan(bet2/2))/(1+tan(bet1/2)*tan(bet2/2)) + * with tan(x/2) = sin(x)/(1+cos(x)) */ + real + somg12 = sin(omg12), domg12 = 1 + cos(omg12), + dbet1 = 1 + cbet1, dbet2 = 1 + cbet2; + alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ), + domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) ); + } else { + /* alp12 = alp2 - alp1, used in atan2 so no need to normalize */ + real + salp12 = salp2 * calp1 - calp2 * salp1, + calp12 = calp2 * calp1 + salp2 * salp1; + /* The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz + * salp12 = -0 and alp12 = -180. However this depends on the sign + * being attached to 0 correctly. The following ensures the correct + * behavior. */ + if (salp12 == 0 && calp12 < 0) { + salp12 = tiny * calp1; + calp12 = -1; + } + alp12 = atan2(salp12, calp12); + } + S12 += g->c2 * alp12; + S12 *= swapp * lonsign * latsign; + /* Convert -0 to 0 */ + S12 += 0; + } + + /* Convert calp, salp to azimuth accounting for lonsign, swapp, latsign. */ + if (swapp < 0) { + swapx(&salp1, &salp2); + swapx(&calp1, &calp2); + if (outmask & GEOD_GEODESICSCALE) + swapx(&M12, &M21); + } + + salp1 *= swapp * lonsign; calp1 *= swapp * latsign; + salp2 *= swapp * lonsign; calp2 *= swapp * latsign; + + if (outmask & GEOD_AZIMUTH) { + /* minus signs give range [-180, 180). 0- converts -0 to +0. */ + azi1 = 0 - atan2(-salp1, calp1) / degree; + azi2 = 0 - atan2(-salp2, calp2) / degree; + } + + if (outmask & GEOD_DISTANCE) + *ps12 = s12; + if (outmask & GEOD_AZIMUTH) { + if (pazi1) *pazi1 = azi1; + if (pazi2) *pazi2 = azi2; + } + if (outmask & GEOD_REDUCEDLENGTH) + *pm12 = m12; + if (outmask & GEOD_GEODESICSCALE) { + if (pM12) *pM12 = M12; + if (pM21) *pM21 = M21; + } + if (outmask & GEOD_AREA) + *pS12 = S12; + + /* Returned value in [0, 180] */ + return a12; +} + +void geod_inverse(const struct geod_geodesic* g, + real lat1, real lon1, real lat2, real lon2, + real* ps12, real* pazi1, real* pazi2) { + geod_geninverse(g, lat1, lon1, lat2, lon2, ps12, pazi1, pazi2, 0, 0, 0, 0); +} + +real SinCosSeries(boolx sinp, real sinx, real cosx, const real c[], int n) { + /* Evaluate + * y = sinp ? sum(c[i] * sin( 2*i * x), i, 1, n) : + * sum(c[i] * cos((2*i+1) * x), i, 0, n-1) + * using Clenshaw summation. N.B. c[0] is unused for sin series + * Approx operation count = (n + 5) mult and (2 * n + 2) add */ + real ar, y0, y1; + c += (n + sinp); /* Point to one beyond last element */ + ar = 2 * (cosx - sinx) * (cosx + sinx); /* 2 * cos(2 * x) */ + y0 = n & 1 ? *--c : 0; y1 = 0; /* accumulators for sum */ + /* Now n is even */ + n /= 2; + while (n--) { + /* Unroll loop x 2, so accumulators return to their original role */ + y1 = ar * y0 - y1 + *--c; + y0 = ar * y1 - y0 + *--c; + } + return sinp + ? 2 * sinx * cosx * y0 /* sin(2 * x) * y0 */ + : cosx * (y0 - y1); /* cos(x) * (y0 - y1) */ +} + +void Lengths(const struct geod_geodesic* g, + real eps, real sig12, + real ssig1, real csig1, real dn1, + real ssig2, real csig2, real dn2, + real cbet1, real cbet2, + real* ps12b, real* pm12b, real* pm0, + boolx scalep, real* pM12, real* pM21, + /* Scratch areas of the right size */ + real C1a[], real C2a[]) { + real s12b = 0, m12b = 0, m0 = 0, M12 = 0, M21 = 0; + real A1m1, AB1, A2m1, AB2, J12; + + /* Return m12b = (reduced length)/b; also calculate s12b = distance/b, + * and m0 = coefficient of secular term in expression for reduced length. */ + C1f(eps, C1a); + C2f(eps, C2a); + A1m1 = A1m1f(eps); + AB1 = (1 + A1m1) * (SinCosSeries(TRUE, ssig2, csig2, C1a, nC1) - + SinCosSeries(TRUE, ssig1, csig1, C1a, nC1)); + A2m1 = A2m1f(eps); + AB2 = (1 + A2m1) * (SinCosSeries(TRUE, ssig2, csig2, C2a, nC2) - + SinCosSeries(TRUE, ssig1, csig1, C2a, nC2)); + m0 = A1m1 - A2m1; + J12 = m0 * sig12 + (AB1 - AB2); + /* Missing a factor of b. + * Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure accurate + * cancellation in the case of coincident points. */ + m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) - csig1 * csig2 * J12; + /* Missing a factor of b */ + s12b = (1 + A1m1) * sig12 + AB1; + if (scalep) { + real csig12 = csig1 * csig2 + ssig1 * ssig2; + real t = g->ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2); + M12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1; + M21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2; + } + *ps12b = s12b; + *pm12b = m12b; + *pm0 = m0; + if (scalep) { + *pM12 = M12; + *pM21 = M21; + } +} + +real Astroid(real x, real y) { + /* Solve k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive root k. + * This solution is adapted from Geocentric::Reverse. */ + real k; + real + p = sq(x), + q = sq(y), + r = (p + q - 1) / 6; + if ( !(q == 0 && r <= 0) ) { + real + /* Avoid possible division by zero when r = 0 by multiplying equations + * for s and t by r^3 and r, resp. */ + S = p * q / 4, /* S = r^3 * s */ + r2 = sq(r), + r3 = r * r2, + /* The discrimant of the quadratic equation for T3. This is zero on + * the evolute curve p^(1/3)+q^(1/3) = 1 */ + disc = S * (S + 2 * r3); + real u = r; + real v, uv, w; + if (disc >= 0) { + real T3 = S + r3, T; + /* Pick the sign on the sqrt to maximize abs(T3). This minimizes loss + * of precision due to cancellation. The result is unchanged because + * of the way the T is used in definition of u. */ + T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc); /* T3 = (r * t)^3 */ + /* N.B. cbrtx always returns the real root. cbrtx(-8) = -2. */ + T = cbrtx(T3); /* T = r * t */ + /* T can be zero; but then r2 / T -> 0. */ + u += T + (T != 0 ? r2 / T : 0); + } else { + /* T is complex, but the way u is defined the result is real. */ + real ang = atan2(sqrt(-disc), -(S + r3)); + /* There are three possible cube roots. We choose the root which + * avoids cancellation. Note that disc < 0 implies that r < 0. */ + u += 2 * r * cos(ang / 3); + } + v = sqrt(sq(u) + q); /* guaranteed positive */ + /* Avoid loss of accuracy when u < 0. */ + uv = u < 0 ? q / (v - u) : u + v; /* u+v, guaranteed positive */ + w = (uv - q) / (2 * v); /* positive? */ + /* Rearrange expression for k to avoid loss of accuracy due to + * subtraction. Division by 0 not possible because uv > 0, w >= 0. */ + k = uv / (sqrt(uv + sq(w)) + w); /* guaranteed positive */ + } else { /* q == 0 && r <= 0 */ + /* y = 0 with |x| <= 1. Handle this case directly. + * for y small, positive root is k = abs(y)/sqrt(1-x^2) */ + k = 0; + } + return k; +} + +real InverseStart(const struct geod_geodesic* g, + real sbet1, real cbet1, real dn1, + real sbet2, real cbet2, real dn2, + real lam12, + real* psalp1, real* pcalp1, + /* Only updated if return val >= 0 */ + real* psalp2, real* pcalp2, + /* Only updated for short lines */ + real* pdnm, + /* Scratch areas of the right size */ + real C1a[], real C2a[]) { + real salp1 = 0, calp1 = 0, salp2 = 0, calp2 = 0, dnm = 0; + + /* Return a starting point for Newton's method in salp1 and calp1 (function + * value is -1). If Newton's method doesn't need to be used, return also + * salp2 and calp2 and function value is sig12. */ + real + sig12 = -1, /* Return value */ + /* bet12 = bet2 - bet1 in [0, pi); bet12a = bet2 + bet1 in (-pi, 0] */ + sbet12 = sbet2 * cbet1 - cbet2 * sbet1, + cbet12 = cbet2 * cbet1 + sbet2 * sbet1; +#if defined(__GNUC__) && __GNUC__ == 4 && \ + (__GNUC_MINOR__ < 6 || defined(__MINGW32__)) + /* Volatile declaration needed to fix inverse cases + * 88.202499451857 0 -88.202499451857 179.981022032992859592 + * 89.262080389218 0 -89.262080389218 179.992207982775375662 + * 89.333123580033 0 -89.333123580032997687 179.99295812360148422 + * which otherwise fail with g++ 4.4.4 x86 -O3 (Linux) + * and g++ 4.4.0 (mingw) and g++ 4.6.1 (tdm mingw). */ + real sbet12a; + { + volatile real xx1 = sbet2 * cbet1; + volatile real xx2 = cbet2 * sbet1; + sbet12a = xx1 + xx2; + } +#else + real sbet12a = sbet2 * cbet1 + cbet2 * sbet1; +#endif + boolx shortline = cbet12 >= 0 && sbet12 < (real)(0.5) && + cbet2 * lam12 < (real)(0.5); + real omg12 = lam12, somg12, comg12, ssig12, csig12; + if (shortline) { + real sbetm2 = sq(sbet1 + sbet2); + /* sin((bet1+bet2)/2)^2 + * = (sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2) */ + sbetm2 /= sbetm2 + sq(cbet1 + cbet2); + dnm = sqrt(1 + g->ep2 * sbetm2); + omg12 /= g->f1 * dnm; + } + somg12 = sin(omg12); comg12 = cos(omg12); + + salp1 = cbet2 * somg12; + calp1 = comg12 >= 0 ? + sbet12 + cbet2 * sbet1 * sq(somg12) / (1 + comg12) : + sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12); + + ssig12 = hypotx(salp1, calp1); + csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12; + + if (shortline && ssig12 < g->etol2) { + /* really short lines */ + salp2 = cbet1 * somg12; + calp2 = sbet12 - cbet1 * sbet2 * + (comg12 >= 0 ? sq(somg12) / (1 + comg12) : 1 - comg12); + SinCosNorm(&salp2, &calp2); + /* Set return value */ + sig12 = atan2(ssig12, csig12); + } else if (fabs(g->n) > (real)(0.1) || /* No astroid calc if too eccentric */ + csig12 >= 0 || + ssig12 >= 6 * fabs(g->n) * pi * sq(cbet1)) { + /* Nothing to do, zeroth order spherical approximation is OK */ + } else { + /* Scale lam12 and bet2 to x, y coordinate system where antipodal point + * is at origin and singular point is at y = 0, x = -1. */ + real y, lamscale, betscale; + /* Volatile declaration needed to fix inverse case + * 56.320923501171 0 -56.320923501171 179.664747671772880215 + * which otherwise fails with g++ 4.4.4 x86 -O3 */ + volatile real x; + if (g->f >= 0) { /* In fact f == 0 does not get here */ + /* x = dlong, y = dlat */ + { + real + k2 = sq(sbet1) * g->ep2, + eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2); + lamscale = g->f * cbet1 * A3f(g, eps) * pi; + } + betscale = lamscale * cbet1; + + x = (lam12 - pi) / lamscale; + y = sbet12a / betscale; + } else { /* f < 0 */ + /* x = dlat, y = dlong */ + real + cbet12a = cbet2 * cbet1 - sbet2 * sbet1, + bet12a = atan2(sbet12a, cbet12a); + real m12b, m0, dummy; + /* In the case of lon12 = 180, this repeats a calculation made in + * Inverse. */ + Lengths(g, g->n, pi + bet12a, + sbet1, -cbet1, dn1, sbet2, cbet2, dn2, + cbet1, cbet2, &dummy, &m12b, &m0, FALSE, + &dummy, &dummy, C1a, C2a); + x = -1 + m12b / (cbet1 * cbet2 * m0 * pi); + betscale = x < -(real)(0.01) ? sbet12a / x : + -g->f * sq(cbet1) * pi; + lamscale = betscale / cbet1; + y = (lam12 - pi) / lamscale; + } + + if (y > -tol1 && x > -1 - xthresh) { + /* strip near cut */ + if (g->f >= 0) { + salp1 = minx((real)(1), -(real)(x)); calp1 = - sqrt(1 - sq(salp1)); + } else { + calp1 = maxx((real)(x > -tol1 ? 0 : -1), (real)(x)); + salp1 = sqrt(1 - sq(calp1)); + } + } else { + /* Estimate alp1, by solving the astroid problem. + * + * Could estimate alpha1 = theta + pi/2, directly, i.e., + * calp1 = y/k; salp1 = -x/(1+k); for f >= 0 + * calp1 = x/(1+k); salp1 = -y/k; for f < 0 (need to check) + * + * However, it's better to estimate omg12 from astroid and use + * spherical formula to compute alp1. This reduces the mean number of + * Newton iterations for astroid cases from 2.24 (min 0, max 6) to 2.12 + * (min 0 max 5). The changes in the number of iterations are as + * follows: + * + * change percent + * 1 5 + * 0 78 + * -1 16 + * -2 0.6 + * -3 0.04 + * -4 0.002 + * + * The histogram of iterations is (m = number of iterations estimating + * alp1 directly, n = number of iterations estimating via omg12, total + * number of trials = 148605): + * + * iter m n + * 0 148 186 + * 1 13046 13845 + * 2 93315 102225 + * 3 36189 32341 + * 4 5396 7 + * 5 455 1 + * 6 56 0 + * + * Because omg12 is near pi, estimate work with omg12a = pi - omg12 */ + real k = Astroid(x, y); + real + omg12a = lamscale * ( g->f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k ); + somg12 = sin(omg12a); comg12 = -cos(omg12a); + /* Update spherical estimate of alp1 using omg12 instead of lam12 */ + salp1 = cbet2 * somg12; + calp1 = sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12); + } + } + /* Sanity check on starting guess. Backwards check allows NaN through. */ + if (!(salp1 <= 0)) + SinCosNorm(&salp1, &calp1); + else { + salp1 = 1; calp1 = 0; + } + + *psalp1 = salp1; + *pcalp1 = calp1; + if (shortline) + *pdnm = dnm; + if (sig12 >= 0) { + *psalp2 = salp2; + *pcalp2 = calp2; + } + return sig12; +} + +real Lambda12(const struct geod_geodesic* g, + real sbet1, real cbet1, real dn1, + real sbet2, real cbet2, real dn2, + real salp1, real calp1, + real* psalp2, real* pcalp2, + real* psig12, + real* pssig1, real* pcsig1, + real* pssig2, real* pcsig2, + real* peps, real* pdomg12, + boolx diffp, real* pdlam12, + /* Scratch areas of the right size */ + real C1a[], real C2a[], real C3a[]) { + real salp2 = 0, calp2 = 0, sig12 = 0, + ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0, domg12 = 0, dlam12 = 0; + real salp0, calp0; + real somg1, comg1, somg2, comg2, omg12, lam12; + real B312, h0, k2; + + if (sbet1 == 0 && calp1 == 0) + /* Break degeneracy of equatorial line. This case has already been + * handled. */ + calp1 = -tiny; + + /* sin(alp1) * cos(bet1) = sin(alp0) */ + salp0 = salp1 * cbet1; + calp0 = hypotx(calp1, salp1 * sbet1); /* calp0 > 0 */ + + /* tan(bet1) = tan(sig1) * cos(alp1) + * tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1) */ + ssig1 = sbet1; somg1 = salp0 * sbet1; + csig1 = comg1 = calp1 * cbet1; + SinCosNorm(&ssig1, &csig1); + /* SinCosNorm(&somg1, &comg1); -- don't need to normalize! */ + + /* Enforce symmetries in the case abs(bet2) = -bet1. Need to be careful + * about this case, since this can yield singularities in the Newton + * iteration. + * sin(alp2) * cos(bet2) = sin(alp0) */ + salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1; + /* calp2 = sqrt(1 - sq(salp2)) + * = sqrt(sq(calp0) - sq(sbet2)) / cbet2 + * and subst for calp0 and rearrange to give (choose positive sqrt + * to give alp2 in [0, pi/2]). */ + calp2 = cbet2 != cbet1 || fabs(sbet2) != -sbet1 ? + sqrt(sq(calp1 * cbet1) + + (cbet1 < -sbet1 ? + (cbet2 - cbet1) * (cbet1 + cbet2) : + (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 : + fabs(calp1); + /* tan(bet2) = tan(sig2) * cos(alp2) + * tan(omg2) = sin(alp0) * tan(sig2). */ + ssig2 = sbet2; somg2 = salp0 * sbet2; + csig2 = comg2 = calp2 * cbet2; + SinCosNorm(&ssig2, &csig2); + /* SinCosNorm(&somg2, &comg2); -- don't need to normalize! */ + + /* sig12 = sig2 - sig1, limit to [0, pi] */ + sig12 = atan2(maxx(csig1 * ssig2 - ssig1 * csig2, (real)(0)), + csig1 * csig2 + ssig1 * ssig2); + + /* omg12 = omg2 - omg1, limit to [0, pi] */ + omg12 = atan2(maxx(comg1 * somg2 - somg1 * comg2, (real)(0)), + comg1 * comg2 + somg1 * somg2); + k2 = sq(calp0) * g->ep2; + eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2); + C3f(g, eps, C3a); + B312 = (SinCosSeries(TRUE, ssig2, csig2, C3a, nC3-1) - + SinCosSeries(TRUE, ssig1, csig1, C3a, nC3-1)); + h0 = -g->f * A3f(g, eps); + domg12 = salp0 * h0 * (sig12 + B312); + lam12 = omg12 + domg12; + + if (diffp) { + if (calp2 == 0) + dlam12 = - 2 * g->f1 * dn1 / sbet1; + else { + real dummy; + Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, + cbet1, cbet2, &dummy, &dlam12, &dummy, + FALSE, &dummy, &dummy, C1a, C2a); + dlam12 *= g->f1 / (calp2 * cbet2); + } + } + + *psalp2 = salp2; + *pcalp2 = calp2; + *psig12 = sig12; + *pssig1 = ssig1; + *pcsig1 = csig1; + *pssig2 = ssig2; + *pcsig2 = csig2; + *peps = eps; + *pdomg12 = domg12; + if (diffp) + *pdlam12 = dlam12; + + return lam12; +} + +real A3f(const struct geod_geodesic* g, real eps) { + /* Evaluate sum(A3x[k] * eps^k, k, 0, nA3x-1) by Horner's method */ + real v = 0; + int i; + for (i = nA3x; i; ) + v = eps * v + g->A3x[--i]; + return v; +} + +void C3f(const struct geod_geodesic* g, real eps, real c[]) { + /* Evaluate C3 coeffs by Horner's method + * Elements c[1] thru c[nC3 - 1] are set */ + int i, j, k; + real mult = 1; + for (j = nC3x, k = nC3 - 1; k; ) { + real t = 0; + for (i = nC3 - k; i; --i) + t = eps * t + g->C3x[--j]; + c[k--] = t; + } + + for (k = 1; k < nC3; ) { + mult *= eps; + c[k++] *= mult; + } +} + +void C4f(const struct geod_geodesic* g, real eps, real c[]) { + /* Evaluate C4 coeffs by Horner's method + * Elements c[0] thru c[nC4 - 1] are set */ + int i, j, k; + real mult = 1; + for (j = nC4x, k = nC4; k; ) { + real t = 0; + for (i = nC4 - k + 1; i; --i) + t = eps * t + g->C4x[--j]; + c[--k] = t; + } + + for (k = 1; k < nC4; ) { + mult *= eps; + c[k++] *= mult; + } +} + +/* Generated by Maxima on 2010-09-04 10:26:17-04:00 */ + +/* The scale factor A1-1 = mean value of (d/dsigma)I1 - 1 */ +real A1m1f(real eps) { + real + eps2 = sq(eps), + t = eps2*(eps2*(eps2+4)+64)/256; + return (t + eps) / (1 - eps); +} + +/* The coefficients C1[l] in the Fourier expansion of B1 */ +void C1f(real eps, real c[]) { + real + eps2 = sq(eps), + d = eps; + c[1] = d*((6-eps2)*eps2-16)/32; + d *= eps; + c[2] = d*((64-9*eps2)*eps2-128)/2048; + d *= eps; + c[3] = d*(9*eps2-16)/768; + d *= eps; + c[4] = d*(3*eps2-5)/512; + d *= eps; + c[5] = -7*d/1280; + d *= eps; + c[6] = -7*d/2048; +} + +/* The coefficients C1p[l] in the Fourier expansion of B1p */ +void C1pf(real eps, real c[]) { + real + eps2 = sq(eps), + d = eps; + c[1] = d*(eps2*(205*eps2-432)+768)/1536; + d *= eps; + c[2] = d*(eps2*(4005*eps2-4736)+3840)/12288; + d *= eps; + c[3] = d*(116-225*eps2)/384; + d *= eps; + c[4] = d*(2695-7173*eps2)/7680; + d *= eps; + c[5] = 3467*d/7680; + d *= eps; + c[6] = 38081*d/61440; +} + +/* The scale factor A2-1 = mean value of (d/dsigma)I2 - 1 */ +real A2m1f(real eps) { + real + eps2 = sq(eps), + t = eps2*(eps2*(25*eps2+36)+64)/256; + return t * (1 - eps) - eps; +} + +/* The coefficients C2[l] in the Fourier expansion of B2 */ +void C2f(real eps, real c[]) { + real + eps2 = sq(eps), + d = eps; + c[1] = d*(eps2*(eps2+2)+16)/32; + d *= eps; + c[2] = d*(eps2*(35*eps2+64)+384)/2048; + d *= eps; + c[3] = d*(15*eps2+80)/768; + d *= eps; + c[4] = d*(7*eps2+35)/512; + d *= eps; + c[5] = 63*d/1280; + d *= eps; + c[6] = 77*d/2048; +} + +/* The scale factor A3 = mean value of (d/dsigma)I3 */ +void A3coeff(struct geod_geodesic* g) { + g->A3x[0] = 1; + g->A3x[1] = (g->n-1)/2; + g->A3x[2] = (g->n*(3*g->n-1)-2)/8; + g->A3x[3] = ((-g->n-3)*g->n-1)/16; + g->A3x[4] = (-2*g->n-3)/64; + g->A3x[5] = -3/(real)(128); +} + +/* The coefficients C3[l] in the Fourier expansion of B3 */ +void C3coeff(struct geod_geodesic* g) { + g->C3x[0] = (1-g->n)/4; + g->C3x[1] = (1-g->n*g->n)/8; + g->C3x[2] = ((3-g->n)*g->n+3)/64; + g->C3x[3] = (2*g->n+5)/128; + g->C3x[4] = 3/(real)(128); + g->C3x[5] = ((g->n-3)*g->n+2)/32; + g->C3x[6] = ((-3*g->n-2)*g->n+3)/64; + g->C3x[7] = (g->n+3)/128; + g->C3x[8] = 5/(real)(256); + g->C3x[9] = (g->n*(5*g->n-9)+5)/192; + g->C3x[10] = (9-10*g->n)/384; + g->C3x[11] = 7/(real)(512); + g->C3x[12] = (7-14*g->n)/512; + g->C3x[13] = 7/(real)(512); + g->C3x[14] = 21/(real)(2560); +} + +/* Generated by Maxima on 2012-10-19 08:02:34-04:00 */ + +/* The coefficients C4[l] in the Fourier expansion of I4 */ +void C4coeff(struct geod_geodesic* g) { + g->C4x[0] = (g->n*(g->n*(g->n*(g->n*(100*g->n+208)+572)+3432)-12012)+30030)/ + 45045; + g->C4x[1] = (g->n*(g->n*(g->n*(64*g->n+624)-4576)+6864)-3003)/15015; + g->C4x[2] = (g->n*((14144-10656*g->n)*g->n-4576)-858)/45045; + g->C4x[3] = ((-224*g->n-4784)*g->n+1573)/45045; + g->C4x[4] = (1088*g->n+156)/45045; + g->C4x[5] = 97/(real)(15015); + g->C4x[6] = (g->n*(g->n*((-64*g->n-624)*g->n+4576)-6864)+3003)/135135; + g->C4x[7] = (g->n*(g->n*(5952*g->n-11648)+9152)-2574)/135135; + g->C4x[8] = (g->n*(5792*g->n+1040)-1287)/135135; + g->C4x[9] = (468-2944*g->n)/135135; + g->C4x[10] = 1/(real)(9009); + g->C4x[11] = (g->n*((4160-1440*g->n)*g->n-4576)+1716)/225225; + g->C4x[12] = ((4992-8448*g->n)*g->n-1144)/225225; + g->C4x[13] = (1856*g->n-936)/225225; + g->C4x[14] = 8/(real)(10725); + g->C4x[15] = (g->n*(3584*g->n-3328)+1144)/315315; + g->C4x[16] = (1024*g->n-208)/105105; + g->C4x[17] = -136/(real)(63063); + g->C4x[18] = (832-2560*g->n)/405405; + g->C4x[19] = -128/(real)(135135); + g->C4x[20] = 128/(real)(99099); +} + +int transit(real lon1, real lon2) { + real lon12; + /* Return 1 or -1 if crossing prime meridian in east or west direction. + * Otherwise return zero. */ + /* Compute lon12 the same way as Geodesic::Inverse. */ + lon1 = AngNormalize(lon1); + lon2 = AngNormalize(lon2); + lon12 = AngDiff(lon1, lon2); + return lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 : + (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0); +} + +int transitdirect(real lon1, real lon2) { + lon1 = fmod(lon1, (real)(720)); + lon2 = fmod(lon2, (real)(720)); + return ( ((lon2 >= 0 && lon2 < 360) || lon2 < -360 ? 0 : 1) - + ((lon1 >= 0 && lon1 < 360) || lon1 < -360 ? 0 : 1) ); +} + +void accini(real s[]) { + /* Initialize an accumulator; this is an array with two elements. */ + s[0] = s[1] = 0; +} + +void acccopy(const real s[], real t[]) { + /* Copy an accumulator; t = s. */ + t[0] = s[0]; t[1] = s[1]; +} + +void accadd(real s[], real y) { + /* Add y to an accumulator. */ + real u, z = sumx(y, s[1], &u); + s[0] = sumx(z, s[0], &s[1]); + if (s[0] == 0) + s[0] = u; + else + s[1] = s[1] + u; +} + +real accsum(const real s[], real y) { + /* Return accumulator + y (but don't add to accumulator). */ + real t[2]; + acccopy(s, t); + accadd(t, y); + return t[0]; +} + +void accneg(real s[]) { + /* Negate an accumulator. */ + s[0] = -s[0]; s[1] = -s[1]; +} + +void geod_polygon_init(struct geod_polygon* p, boolx polylinep) { + p->lat0 = p->lon0 = p->lat = p->lon = NaN; + p->polyline = (polylinep != 0); + accini(p->P); + accini(p->A); + p->num = p->crossings = 0; +} + +void geod_polygon_addpoint(const struct geod_geodesic* g, + struct geod_polygon* p, + real lat, real lon) { + lon = AngNormalize(lon); + if (p->num == 0) { + p->lat0 = p->lat = lat; + p->lon0 = p->lon = lon; + } else { + real s12, S12; + geod_geninverse(g, p->lat, p->lon, lat, lon, + &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12); + accadd(p->P, s12); + if (!p->polyline) { + accadd(p->A, S12); + p->crossings += transit(p->lon, lon); + } + p->lat = lat; p->lon = lon; + } + ++p->num; +} + +void geod_polygon_addedge(const struct geod_geodesic* g, + struct geod_polygon* p, + real azi, real s) { + if (p->num) { /* Do nothing is num is zero */ + real lat, lon, S12; + geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_NOWRAP, s, + &lat, &lon, 0, + 0, 0, 0, 0, p->polyline ? 0 : &S12); + accadd(p->P, s); + if (!p->polyline) { + accadd(p->A, S12); + p->crossings += transitdirect(p->lon, lon); + } + p->lat = lat; p->lon = lon; + ++p->num; + } +} + +unsigned geod_polygon_compute(const struct geod_geodesic* g, + const struct geod_polygon* p, + boolx reverse, boolx sign, + real* pA, real* pP) { + real s12, S12, t[2], area0; + int crossings; + if (p->num < 2) { + if (pP) *pP = 0; + if (!p->polyline && pA) *pA = 0; + return p->num; + } + if (p->polyline) { + if (pP) *pP = p->P[0]; + return p->num; + } + geod_geninverse(g, p->lat, p->lon, p->lat0, p->lon0, + &s12, 0, 0, 0, 0, 0, &S12); + if (pP) *pP = accsum(p->P, s12); + acccopy(p->A, t); + accadd(t, S12); + crossings = p->crossings + transit(p->lon, p->lon0); + area0 = 4 * pi * g->c2; + if (crossings & 1) + accadd(t, (t[0] < 0 ? 1 : -1) * area0/2); + /* area is with the clockwise sense. If !reverse convert to + * counter-clockwise convention. */ + if (!reverse) + accneg(t); + /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */ + if (sign) { + if (t[0] > area0/2) + accadd(t, -area0); + else if (t[0] <= -area0/2) + accadd(t, +area0); + } else { + if (t[0] >= area0) + accadd(t, -area0); + else if (t[0] < 0) + accadd(t, +area0); + } + if (pA) *pA = 0 + t[0]; + return p->num; +} + +unsigned geod_polygon_testpoint(const struct geod_geodesic* g, + const struct geod_polygon* p, + real lat, real lon, + boolx reverse, boolx sign, + real* pA, real* pP) { + real perimeter, tempsum, area0; + int crossings, i; + unsigned num = p->num + 1; + if (num == 1) { + if (pP) *pP = 0; + if (!p->polyline && pA) *pA = 0; + return num; + } + perimeter = p->P[0]; + tempsum = p->polyline ? 0 : p->A[0]; + crossings = p->crossings; + for (i = 0; i < (p->polyline ? 1 : 2); ++i) { + real s12, S12; + geod_geninverse(g, + i == 0 ? p->lat : lat, i == 0 ? p->lon : lon, + i != 0 ? p->lat0 : lat, i != 0 ? p->lon0 : lon, + &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12); + perimeter += s12; + if (!p->polyline) { + tempsum += S12; + crossings += transit(i == 0 ? p->lon : lon, + i != 0 ? p->lon0 : lon); + } + } + + if (pP) *pP = perimeter; + if (p->polyline) + return num; + + area0 = 4 * pi * g->c2; + if (crossings & 1) + tempsum += (tempsum < 0 ? 1 : -1) * area0/2; + /* area is with the clockwise sense. If !reverse convert to + * counter-clockwise convention. */ + if (!reverse) + tempsum *= -1; + /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */ + if (sign) { + if (tempsum > area0/2) + tempsum -= area0; + else if (tempsum <= -area0/2) + tempsum += area0; + } else { + if (tempsum >= area0) + tempsum -= area0; + else if (tempsum < 0) + tempsum += area0; + } + if (pA) *pA = 0 + tempsum; + return num; +} + +unsigned geod_polygon_testedge(const struct geod_geodesic* g, + const struct geod_polygon* p, + real azi, real s, + boolx reverse, boolx sign, + real* pA, real* pP) { + real perimeter, tempsum, area0; + int crossings; + unsigned num = p->num + 1; + if (num == 1) { /* we don't have a starting point! */ + if (pP) *pP = NaN; + if (!p->polyline && pA) *pA = NaN; + return 0; + } + perimeter = p->P[0] + s; + if (p->polyline) { + if (pP) *pP = perimeter; + return num; + } + + tempsum = p->A[0]; + crossings = p->crossings; + { + real lat, lon, s12, S12; + geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_NOWRAP, s, + &lat, &lon, 0, + 0, 0, 0, 0, &S12); + tempsum += S12; + crossings += transitdirect(p->lon, lon); + geod_geninverse(g, lat, lon, p->lat0, p->lon0, + &s12, 0, 0, 0, 0, 0, &S12); + perimeter += s12; + tempsum += S12; + crossings += transit(lon, p->lon0); + } + + area0 = 4 * pi * g->c2; + if (crossings & 1) + tempsum += (tempsum < 0 ? 1 : -1) * area0/2; + /* area is with the clockwise sense. If !reverse convert to + * counter-clockwise convention. */ + if (!reverse) + tempsum *= -1; + /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */ + if (sign) { + if (tempsum > area0/2) + tempsum -= area0; + else if (tempsum <= -area0/2) + tempsum += area0; + } else { + if (tempsum >= area0) + tempsum -= area0; + else if (tempsum < 0) + tempsum += area0; + } + if (pP) *pP = perimeter; + if (pA) *pA = 0 + tempsum; + return num; +} + +void geod_polygonarea(const struct geod_geodesic* g, + real lats[], real lons[], int n, + real* pA, real* pP) { + int i; + struct geod_polygon p; + geod_polygon_init(&p, FALSE); + for (i = 0; i < n; ++i) + geod_polygon_addpoint(g, &p, lats[i], lons[i]); + geod_polygon_compute(g, &p, FALSE, TRUE, pA, pP); +} + +/** @endcond */ diff --git a/liblwgeom/geodesic.h b/liblwgeom/geodesic.h new file mode 100644 index 00000000000..6b2afc5a620 --- /dev/null +++ b/liblwgeom/geodesic.h @@ -0,0 +1,781 @@ +/** + * \file geodesic.h + * \brief Header for the geodesic routines in C + * + * This an implementation in C of the geodesic algorithms described in + * - C. F. F. Karney, + * + * Algorithms for geodesics, + * J. Geodesy 87, 43--55 (2013); + * DOI: + * 10.1007/s00190-012-0578-z; + * addenda: + * geod-addenda.html. + * . + * The principal advantages of these algorithms over previous ones (e.g., + * Vincenty, 1975) are + * - accurate to round off for |f| < 1/50; + * - the solution of the inverse problem is always found; + * - differential and integral properties of geodesics are computed. + * + * The shortest path between two points on the ellipsoid at (\e lat1, \e + * lon1) and (\e lat2, \e lon2) is called the geodesic. Its length is + * \e s12 and the geodesic from point 1 to point 2 has forward azimuths + * \e azi1 and \e azi2 at the two end points. + * + * Traditionally two geodesic problems are considered: + * - the direct problem -- given \e lat1, \e lon1, \e s12, and \e azi1, + * determine \e lat2, \e lon2, and \e azi2. This is solved by the function + * geod_direct(). + * - the inverse problem -- given \e lat1, \e lon1, and \e lat2, \e lon2, + * determine \e s12, \e azi1, and \e azi2. This is solved by the function + * geod_inverse(). + * + * The ellipsoid is specified by its equatorial radius \e a (typically in + * meters) and flattening \e f. The routines are accurate to round off with + * double precision arithmetic provided that |f| < 1/50; for the + * WGS84 ellipsoid, the errors are less than 15 nanometers. (Reasonably + * accurate results are obtained for |f| < 1/5.) For a prolate + * ellipsoid, specify \e f < 0. + * + * The routines also calculate several other quantities of interest + * - \e 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 + * quadrilateral with corners (\e lat1,\e lon1), (0,\e lon1), (0,\e lon2), + * and (\e lat2,\e lon2). + * - \e m12, the reduced length of the geodesic is defined such that if + * the initial azimuth is perturbed by \e dazi1 (radians) then the + * second point is displaced by \e m12 \e dazi1 in the direction + * perpendicular to the geodesic. On a curved surface the reduced + * length obeys a symmetry relation, \e m12 + \e m21 = 0. On a flat + * surface, we have \e m12 = \e s12. + * - \e M12 and \e M21 are geodesic scales. If two geodesics are + * parallel at point 1 and separated by a small distance \e dt, then + * they are separated by a distance \e M12 \e dt at point 2. \e M21 + * is defined similarly (with the geodesics being parallel to one + * another at point 2). On a flat surface, we have \e M12 = \e M21 + * = 1. + * - \e a12 is the arc length on the auxiliary sphere. This is a + * construct for converting the problem to one in spherical + * trigonometry. \e a12 is measured in degrees. The spherical arc + * length from one equator crossing to the next is always 180°. + * + * If points 1, 2, and 3 lie on a single geodesic, then the following + * addition rules hold: + * - \e s13 = \e s12 + \e s23 + * - \e a13 = \e a12 + \e a23 + * - \e S13 = \e S12 + \e S23 + * - \e m13 = \e m12 \e M23 + \e m23 \e M21 + * - \e M13 = \e M12 \e M23 − (1 − \e M12 \e M21) \e + * m23 / \e m12 + * - \e M31 = \e M32 \e M21 − (1 − \e M23 \e M32) \e + * m12 / \e m23 + * + * 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: + * - \e lat1 = −\e lat2 (with neither point at a pole). If \e azi1 = + * \e azi2, the geodesic is unique. Otherwise there are two geodesics + * and the second one is obtained by setting [\e azi1, \e azi2] = [\e + * azi2, \e azi1], [\e M12, \e M21] = [\e M21, \e M12], \e S12 = + * −\e S12. (This occurs when the longitude difference is near + * ±180° for oblate ellipsoids.) + * - \e lon2 = \e lon1 ± 180° (with neither point at a pole). + * If \e azi1 = 0° or ±180°, the geodesic is unique. + * Otherwise there are two geodesics and the second one is obtained by + * setting [\e azi1, \e azi2] = [−\e azi1, −\e azi2], \e S12 + * = −\e S12. (This occurs when \e lat2 is near −\e lat1 for + * prolate ellipsoids.) + * - Points 1 and 2 at opposite poles. There are infinitely many + * geodesics which can be generated by setting [\e azi1, \e azi2] = + * [\e azi1, \e azi2] + [\e d, −\e d], for arbitrary \e d. (For + * spheres, this prescription applies when points 1 and 2 are + * antipodal.) + * - \e s12 = 0 (coincident points). There are infinitely many geodesics + * which can be generated by setting [\e azi1, \e azi2] = [\e azi1, \e + * azi2] + [\e d, \e d], for arbitrary \e d. + * + * These routines are a simple transcription of the corresponding C++ classes + * in GeographicLib. The "class + * data" is represented by the structs geod_geodesic, geod_geodesicline, + * geod_polygon and pointers to these objects are passed as initial arguments + * to the member functions. Most of the internal comments have been retained. + * However, in the process of transcription some documentation has been lost + * and the documentation for the C++ classes, GeographicLib::Geodesic, + * GeographicLib::GeodesicLine, and GeographicLib::PolygonAreaT, should be + * consulted. The C++ code remains the "reference implementation". Think + * twice about restructuring the internals of the C code since this may make + * porting fixes from the C++ code more difficult. + * + * Copyright (c) Charles Karney (2012-2014) and licensed + * under the MIT/X11 License. For more information, see + * http://geographiclib.sourceforge.net/ + * + * This library was distributed with + * GeographicLib 1.40. + **********************************************************************/ + +#if !defined(GEODESIC_H) +#define GEODESIC_H 1 + +/** + * The major version of the geodesic library. (This tracks the version of + * GeographicLib.) + **********************************************************************/ +#define GEODESIC_VERSION_MAJOR 1 +/** + * The minor version of the geodesic library. (This tracks the version of + * GeographicLib.) + **********************************************************************/ +#define GEODESIC_VERSION_MINOR 40 +/** + * The patch level of the geodesic library. (This tracks the version of + * GeographicLib.) + **********************************************************************/ +#define GEODESIC_VERSION_PATCH 0 + +#if defined(__cplusplus) +extern "C" { +#endif + + /** + * The struct containing information about the ellipsoid. This must be + * initialized by geod_init() before use. + **********************************************************************/ + struct geod_geodesic { + double a; /**< the equatorial radius */ + double f; /**< the flattening */ + /**< @cond SKIP */ + double f1, e2, ep2, n, b, c2, etol2; + double A3x[6], C3x[15], C4x[21]; + /**< @endcond */ + }; + + /** + * The struct containing information about a single geodesic. This must be + * initialized by geod_lineinit() before use. + **********************************************************************/ + struct geod_geodesicline { + double lat1; /**< the starting latitude */ + double lon1; /**< the starting longitude */ + double azi1; /**< the starting azimuth */ + double a; /**< the equatorial radius */ + double f; /**< the flattening */ + /**< @cond SKIP */ + double b, c2, f1, salp0, calp0, k2, + salp1, calp1, ssig1, csig1, dn1, stau1, ctau1, somg1, comg1, + A1m1, A2m1, A3c, B11, B21, B31, A4, B41; + double C1a[6+1], C1pa[6+1], C2a[6+1], C3a[6], C4a[6]; + /**< @endcond */ + unsigned caps; /**< the capabilities */ + }; + + /** + * The struct for accumulating information about a geodesic polygon. This is + * used for computing the perimeter and area of a polygon. This must be + * initialized by geod_polygon_init() before use. + **********************************************************************/ + struct geod_polygon { + double lat; /**< the current latitude */ + double lon; /**< the current longitude */ + /**< @cond SKIP */ + double lat0; + double lon0; + double A[2]; + double P[2]; + int polyline; + int crossings; + /**< @endcond */ + unsigned num; /**< the number of points so far */ + }; + + /** + * Initialize a geod_geodesic object. + * + * @param[out] g a pointer to the object to be initialized. + * @param[in] a the equatorial radius (meters). + * @param[in] f the flattening. + **********************************************************************/ + void geod_init(struct geod_geodesic* g, double a, double f); + + /** + * Initialize a geod_geodesicline object. + * + * @param[out] l a pointer to the object to be initialized. + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] lat1 latitude of point 1 (degrees). + * @param[in] lon1 longitude of point 1 (degrees). + * @param[in] azi1 azimuth at point 1 (degrees). + * @param[in] caps bitor'ed combination of geod_mask() values specifying the + * capabilities the geod_geodesicline object should possess, i.e., which + * quantities can be returned in calls to geod_position() and + * geod_genposition(). + * + * \e g must have been initialized with a call to geod_init(). \e lat1 + * should be in the range [−90°, 90°]; \e lon1 and \e azi1 + * should be in the range [−540°, 540°). + * + * The geod_mask values are [see geod_mask()]: + * - \e caps |= GEOD_LATITUDE for the latitude \e lat2; this is + * added automatically, + * - \e caps |= GEOD_LONGITUDE for the latitude \e lon2, + * - \e caps |= GEOD_AZIMUTH for the latitude \e azi2; this is + * added automatically, + * - \e caps |= GEOD_DISTANCE for the distance \e s12, + * - \e caps |= GEOD_REDUCEDLENGTH for the reduced length \e m12, + * - \e caps |= GEOD_GEODESICSCALE for the geodesic scales \e M12 + * and \e M21, + * - \e caps |= GEOD_AREA for the area \e S12, + * - \e caps |= GEOD_DISTANCE_IN permits the length of the + * geodesic to be given in terms of \e s12; without this capability the + * length can only be specified in terms of arc length. + * . + * A value of \e caps = 0 is treated as GEOD_LATITUDE | GEOD_LONGITUDE | + * GEOD_AZIMUTH | GEOD_DISTANCE_IN (to support the solution of the "standard" + * direct problem). + **********************************************************************/ + void geod_lineinit(struct geod_geodesicline* l, + const struct geod_geodesic* g, + double lat1, double lon1, double azi1, unsigned caps); + + /** + * Solve the direct geodesic problem. + * + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] lat1 latitude of point 1 (degrees). + * @param[in] lon1 longitude of point 1 (degrees). + * @param[in] azi1 azimuth at point 1 (degrees). + * @param[in] s12 distance between point 1 and point 2 (meters); it can be + * negative. + * @param[out] plat2 pointer to the latitude of point 2 (degrees). + * @param[out] plon2 pointer to the longitude of point 2 (degrees). + * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). + * + * \e g must have been initialized with a call to geod_init(). \e lat1 + * should be in the range [−90°, 90°]; \e lon1 and \e azi1 + * should be in the range [−540°, 540°). The values of \e lon2 + * and \e azi2 returned are in the range [−180°, 180°). Any of + * the "return" arguments \e plat2, etc., may be replaced by 0, if you do not + * need some quantities computed. + * + * If either point is at a pole, the azimuth is defined by keeping the + * longitude fixed, writing \e 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°.) + * + * Example, determine the point 10000 km NE of JFK: + @code + struct geod_geodesic g; + double lat, lon; + geod_init(&g, 6378137, 1/298.257223563); + geod_direct(&g, 40.64, -73.78, 45.0, 10e6, &lat, &lon, 0); + printf("%.5f %.5f\n", lat, lon); + @endcode + **********************************************************************/ + void geod_direct(const struct geod_geodesic* g, + double lat1, double lon1, double azi1, double s12, + double* plat2, double* plon2, double* pazi2); + + /** + * Solve the inverse geodesic problem. + * + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] lat1 latitude of point 1 (degrees). + * @param[in] lon1 longitude of point 1 (degrees). + * @param[in] lat2 latitude of point 2 (degrees). + * @param[in] lon2 longitude of point 2 (degrees). + * @param[out] ps12 pointer to the distance between point 1 and point 2 + * (meters). + * @param[out] pazi1 pointer to the azimuth at point 1 (degrees). + * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). + * + * \e g must have been initialized with a call to geod_init(). \e lat1 + * and \e lat2 should be in the range [−90°, 90°]; \e lon1 and + * \e lon2 should be in the range [−540°, 540°). The values of + * \e azi1 and \e azi2 returned are in the range [−180°, 180°). + * Any of the "return" arguments \e ps12, etc., may be replaced by 0, if you + * do not need some quantities computed. + * + * If either point is at a pole, the azimuth is defined by keeping the + * longitude fixed, writing \e lat = ±(90° − ε), and + * 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. + * + * Example, determine the distance between JFK and Singapore Changi Airport: + @code + struct geod_geodesic g; + double s12; + geod_init(&g, 6378137, 1/298.257223563); + geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, 0, 0); + printf("%.3f\n", s12); + @endcode + **********************************************************************/ + void geod_inverse(const struct geod_geodesic* g, + double lat1, double lon1, double lat2, double lon2, + double* ps12, double* pazi1, double* pazi2); + + /** + * Compute the position along a geod_geodesicline. + * + * @param[in] l a pointer to the geod_geodesicline object specifying the + * geodesic line. + * @param[in] s12 distance between point 1 and point 2 (meters); it can be + * negative. + * @param[out] plat2 pointer to the latitude of point 2 (degrees). + * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires + * that \e l was initialized with \e caps |= GEOD_LONGITUDE. + * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). + * + * \e l must have been initialized with a call to geod_lineinit() with \e + * caps |= GEOD_DISTANCE_IN. The values of \e lon2 and \e azi2 returned are + * in the range [−180°, 180°). Any of the "return" arguments + * \e plat2, etc., may be replaced by 0, if you do not need some quantities + * computed. + * + * Example, compute way points between JFK and Singapore Changi Airport + * the "obvious" way using geod_direct(): + @code + struct geod_geodesic g; + double s12, azi1, lat[101],lon[101]; + int i; + geod_init(&g, 6378137, 1/298.257223563); + geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0); + for (i = 0; i < 101; ++i) { + geod_direct(&g, 40.64, -73.78, azi1, i * s12 * 0.01, lat + i, lon + i, 0); + printf("%.5f %.5f\n", lat[i], lon[i]); + } + @endcode + * A faster way using geod_position(): + @code + struct geod_geodesic g; + struct geod_geodesicline l; + double s12, azi1, lat[101],lon[101]; + int i; + geod_init(&g, 6378137, 1/298.257223563); + geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0); + geod_lineinit(&l, &g, 40.64, -73.78, azi1, 0); + for (i = 0; i < 101; ++i) { + geod_position(&l, i * s12 * 0.01, lat + i, lon + i, 0); + printf("%.5f %.5f\n", lat[i], lon[i]); + } + @endcode + **********************************************************************/ + void geod_position(const struct geod_geodesicline* l, double s12, + double* plat2, double* plon2, double* pazi2); + + /** + * The general direct geodesic problem. + * + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] lat1 latitude of point 1 (degrees). + * @param[in] lon1 longitude of point 1 (degrees). + * @param[in] azi1 azimuth at point 1 (degrees). + * @param[in] flags bitor'ed combination of geod_flags(); \e flags & + * GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags & + * GEOD_LONG_NOWRAP prevents the value of \e lon2 being wrapped into + * the range [−180°, 180°). + * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, 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. + * @param[out] plat2 pointer to the latitude of point 2 (degrees). + * @param[out] plon2 pointer to the longitude of point 2 (degrees). + * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). + * @param[out] ps12 pointer to the distance between point 1 and point 2 + * (meters). + * @param[out] pm12 pointer to the reduced length of geodesic (meters). + * @param[out] pM12 pointer to the geodesic scale of point 2 relative to + * point 1 (dimensionless). + * @param[out] pM21 pointer to the geodesic scale of point 1 relative to + * point 2 (dimensionless). + * @param[out] pS12 pointer to the area under the geodesic + * (meters2). + * @return \e a12 arc length of between point 1 and point 2 (degrees). + * + * \e g must have been initialized with a call to geod_init(). \e lat1 + * should be in the range [−90°, 90°]; \e lon1 and \e azi1 + * should be in the range [−540°, 540°). The function + * value \e a12 equals \e s12_a12 if \e flags & GEOD_ARCMODE. Any of the + * "return" arguments \e plat2, etc., may be replaced by 0, if you do not + * need some quantities computed. + * + * With \e flags & GEOD_LONG_NOWRAP bit set, the quantity \e lon2 − + * \e lon1 indicates how many times the geodesic wrapped around the + * ellipsoid. Because \e lon2 might be outside the normal allowed range + * for longitudes, [−540°, 540°), be sure to normalize it, + * e.g., with fmod(\e lon2, 360.0) before using it in subsequent + * calculations + **********************************************************************/ + double geod_gendirect(const struct geod_geodesic* g, + double lat1, double lon1, double azi1, + unsigned flags, double s12_a12, + double* plat2, double* plon2, double* pazi2, + double* ps12, double* pm12, double* pM12, double* pM21, + double* pS12); + + /** + * The general inverse geodesic calculation. + * + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] lat1 latitude of point 1 (degrees). + * @param[in] lon1 longitude of point 1 (degrees). + * @param[in] lat2 latitude of point 2 (degrees). + * @param[in] lon2 longitude of point 2 (degrees). + * @param[out] ps12 pointer to the distance between point 1 and point 2 + * (meters). + * @param[out] pazi1 pointer to the azimuth at point 1 (degrees). + * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). + * @param[out] pm12 pointer to the reduced length of geodesic (meters). + * @param[out] pM12 pointer to the geodesic scale of point 2 relative to + * point 1 (dimensionless). + * @param[out] pM21 pointer to the geodesic scale of point 1 relative to + * point 2 (dimensionless). + * @param[out] pS12 pointer to the area under the geodesic + * (meters2). + * @return \e a12 arc length of between point 1 and point 2 (degrees). + * + * \e g must have been initialized with a call to geod_init(). \e lat1 + * and \e lat2 should be in the range [−90°, 90°]; \e lon1 and + * \e lon2 should be in the range [−540°, 540°). Any of the + * "return" arguments \e ps12, etc., may be replaced by 0, if you do not need + * some quantities computed. + **********************************************************************/ + double geod_geninverse(const struct geod_geodesic* g, + double lat1, double lon1, double lat2, double lon2, + double* ps12, double* pazi1, double* pazi2, + double* pm12, double* pM12, double* pM21, + double* pS12); + + /** + * The general position function. + * + * @param[in] l a pointer to the geod_geodesicline object specifying the + * geodesic line. + * @param[in] flags bitor'ed combination of geod_flags(); \e flags & + * GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags & + * GEOD_LONG_NOWRAP prevents the value of \e lon2 being wrapped into + * the range [−180°, 180°); if \e flags & GEOD_ARCMODE is + * 0, then \e l must have been initialized with \e caps |= + * GEOD_DISTANCE_IN. + * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, 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. + * @param[out] plat2 pointer to the latitude of point 2 (degrees). + * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires + * that \e l was initialized with \e caps |= GEOD_LONGITUDE. + * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). + * @param[out] ps12 pointer to the distance between point 1 and point 2 + * (meters); requires that \e l was initialized with \e caps |= + * GEOD_DISTANCE. + * @param[out] pm12 pointer to the reduced length of geodesic (meters); + * requires that \e l was initialized with \e caps |= GEOD_REDUCEDLENGTH. + * @param[out] pM12 pointer to the geodesic scale of point 2 relative to + * point 1 (dimensionless); requires that \e l was initialized with \e caps + * |= GEOD_GEODESICSCALE. + * @param[out] pM21 pointer to the geodesic scale of point 1 relative to + * point 2 (dimensionless); requires that \e l was initialized with \e caps + * |= GEOD_GEODESICSCALE. + * @param[out] pS12 pointer to the area under the geodesic + * (meters2); requires that \e l was initialized with \e caps |= + * GEOD_AREA. + * @return \e a12 arc length of between point 1 and point 2 (degrees). + * + * \e l must have been initialized with a call to geod_lineinit() with \e + * caps |= GEOD_DISTANCE_IN. The value \e azi2 returned is in the range + * [−180°, 180°). Any of the "return" arguments \e plat2, + * etc., may be replaced by 0, if you do not need some quantities + * computed. Requesting a value which \e l is not capable of computing + * is not an error; the corresponding argument will not be altered. + * + * With \e flags & GEOD_LONG_NOWRAP bit set, the quantity \e lon2 − + * \e lon1 indicates how many times the geodesic wrapped around the + * ellipsoid. Because \e lon2 might be outside the normal allowed range + * for longitudes, [−540°, 540°), be sure to normalize it, + * e.g., with fmod(\e lon2, 360.0) before using it in subsequent + * calculations + * + * Example, compute way points between JFK and Singapore Changi Airport + * using geod_genposition(). In this example, the points are evenly space in + * arc length (and so only approximately equally space in distance). This is + * faster than using geod_position() would be appropriate if drawing the path + * on a map. + @code + struct geod_geodesic g; + struct geod_geodesicline l; + double a12, azi1, lat[101], lon[101]; + int i; + geod_init(&g, 6378137, 1/298.257223563); + a12 = geod_geninverse(&g, 40.64, -73.78, 1.36, 103.99, + 0, &azi1, 0, 0, 0, 0, 0); + geod_lineinit(&l, &g, 40.64, -73.78, azi1, GEOD_LATITUDE | GEOD_LONGITUDE); + for (i = 0; i < 101; ++i) { + geod_genposition(&l, 1, i * a12 * 0.01, + lat + i, lon + i, 0, 0, 0, 0, 0, 0); + printf("%.5f %.5f\n", lat[i], lon[i]); + } + @endcode + **********************************************************************/ + double geod_genposition(const struct geod_geodesicline* l, + unsigned flags, double s12_a12, + double* plat2, double* plon2, double* pazi2, + double* ps12, double* pm12, + double* pM12, double* pM21, + double* pS12); + + /** + * Initialize a geod_polygon object. + * + * @param[out] p a pointer to the object to be initialized. + * @param[in] polylinep non-zero if a polyline instead of a polygon. + * + * If \e polylinep is zero, then the sequence of vertices and edges added by + * geod_polygon_addpoint() and geod_polygon_addedge() define a polygon and + * the perimeter and area are returned by geod_polygon_compute(). If \e + * polylinep is non-zero, then the vertices and edges define a polyline and + * only the perimeter is returned by geod_polygon_compute(). + * + * The area and perimeter are accumulated at two times the standard floating + * point precision to guard against the loss of accuracy with many-sided + * polygons. At any point you can ask for the perimeter and area so far. + * + * An example of the use of this function is given in the documentation for + * geod_polygon_compute(). + **********************************************************************/ + void geod_polygon_init(struct geod_polygon* p, int polylinep); + + /** + * Add a point to the polygon or polyline. + * + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in,out] p a pointer to the geod_polygon object specifying the + * polygon. + * @param[in] lat the latitude of the point (degrees). + * @param[in] lon the longitude of the point (degrees). + * + * \e g and \e p must have been initialized with calls to geod_init() and + * geod_polygon_init(), respectively. The same \e g must be used for all the + * points and edges in a polygon. \e lat should be in the range + * [−90°, 90°] and \e lon should be in the range + * [−540°, 540°). + * + * An example of the use of this function is given in the documentation for + * geod_polygon_compute(). + **********************************************************************/ + void geod_polygon_addpoint(const struct geod_geodesic* g, + struct geod_polygon* p, + double lat, double lon); + + /** + * Add an edge to the polygon or polyline. + * + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in,out] p a pointer to the geod_polygon object specifying the + * polygon. + * @param[in] azi azimuth at current point (degrees). + * @param[in] s distance from current point to next point (meters). + * + * \e g and \e p must have been initialized with calls to geod_init() and + * geod_polygon_init(), respectively. The same \e g must be used for all the + * points and edges in a polygon. \e azi should be in the range + * [−540°, 540°). This does nothing if no points have been + * added yet. The \e lat and \e lon fields of \e p give the location of + * the new vertex. + **********************************************************************/ + void geod_polygon_addedge(const struct geod_geodesic* g, + struct geod_polygon* p, + double azi, double s); + + /** + * Return the results for a polygon. + * + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] p a pointer to the geod_polygon object specifying the polygon. + * @param[in] reverse if non-zero then clockwise (instead of + * counter-clockwise) traversal counts as a positive area. + * @param[in] sign if non-zero then return a signed result for the area if + * the polygon is traversed in the "wrong" direction instead of returning + * the area for the rest of the earth. + * @param[out] pA pointer to the area of the polygon (meters2); + * only set if \e polyline is non-zero in the call to geod_polygon_init(). + * @param[out] pP pointer to the perimeter of the polygon or length of the + * polyline (meters). + * @return the number of points. + * + * The area and perimeter are accumulated at two times the standard floating + * point precision to guard against the loss of accuracy with many-sided + * polygons. Only simple polygons (which are not self-intersecting) are + * allowed. There's no need to "close" the polygon by repeating the first + * vertex. Set \e pA or \e pP to zero, if you do not want the corresponding + * quantity returned. + * + * Example, compute the perimeter and area of the geodesic triangle with + * vertices (0°N,0°E), (0°N,90°E), (90°N,0°E). + @code + double A, P; + int n; + struct geod_geodesic g; + struct geod_polygon p; + geod_init(&g, 6378137, 1/298.257223563); + geod_polygon_init(&p, 0); + + geod_polygon_addpoint(&g, &p, 0, 0); + geod_polygon_addpoint(&g, &p, 0, 90); + geod_polygon_addpoint(&g, &p, 90, 0); + n = geod_polygon_compute(&g, &p, 0, 1, &A, &P); + printf("%d %.8f %.3f\n", n, P, A); + @endcode + **********************************************************************/ + unsigned geod_polygon_compute(const struct geod_geodesic* g, + const struct geod_polygon* p, + int reverse, int sign, + double* pA, double* pP); + + /** + * Return the results assuming a tentative final test point is added; + * however, the data for the test point is not saved. This lets you report a + * running result for the perimeter and area as the user moves the mouse + * cursor. Ordinary floating point arithmetic is used to accumulate the data + * for the test point; thus the area and perimeter returned are less accurate + * than if geod_polygon_addpoint() and geod_polygon_compute() are used. + * + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] p a pointer to the geod_polygon object specifying the polygon. + * @param[in] lat the latitude of the test point (degrees). + * @param[in] lon the longitude of the test point (degrees). + * @param[in] reverse if non-zero then clockwise (instead of + * counter-clockwise) traversal counts as a positive area. + * @param[in] sign if non-zero then return a signed result for the area if + * the polygon is traversed in the "wrong" direction instead of returning + * the area for the rest of the earth. + * @param[out] pA pointer to the area of the polygon (meters2); + * only set if \e polyline is non-zero in the call to geod_polygon_init(). + * @param[out] pP pointer to the perimeter of the polygon or length of the + * polyline (meters). + * @return the number of points. + * + * \e lat should be in the range [−90°, 90°] and \e + * lon should be in the range [−540°, 540°). + **********************************************************************/ + unsigned geod_polygon_testpoint(const struct geod_geodesic* g, + const struct geod_polygon* p, + double lat, double lon, + int reverse, int sign, + double* pA, double* pP); + + /** + * Return the results assuming a tentative final test point is added via an + * azimuth and distance; however, the data for the test point is not saved. + * This lets you report a running result for the perimeter and area as the + * user moves the mouse cursor. Ordinary floating point arithmetic is used + * to accumulate the data for the test point; thus the area and perimeter + * returned are less accurate than if geod_polygon_addedge() and + * geod_polygon_compute() are used. + * + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] p a pointer to the geod_polygon object specifying the polygon. + * @param[in] azi azimuth at current point (degrees). + * @param[in] s distance from current point to final test point (meters). + * @param[in] reverse if non-zero then clockwise (instead of + * counter-clockwise) traversal counts as a positive area. + * @param[in] sign if non-zero then return a signed result for the area if + * the polygon is traversed in the "wrong" direction instead of returning + * the area for the rest of the earth. + * @param[out] pA pointer to the area of the polygon (meters2); + * only set if \e polyline is non-zero in the call to geod_polygon_init(). + * @param[out] pP pointer to the perimeter of the polygon or length of the + * polyline (meters). + * @return the number of points. + * + * \e azi should be in the range [−540°, 540°). + **********************************************************************/ + unsigned geod_polygon_testedge(const struct geod_geodesic* g, + const struct geod_polygon* p, + double azi, double s, + int reverse, int sign, + double* pA, double* pP); + + /** + * A simple interface for computing the area of a geodesic polygon. + * + * @param[in] g a pointer to the geod_geodesic object specifying the + * ellipsoid. + * @param[in] lats an array of latitudes of the polygon vertices (degrees). + * @param[in] lons an array of longitudes of the polygon vertices (degrees). + * @param[in] n the number of vertices. + * @param[out] pA pointer to the area of the polygon (meters2). + * @param[out] pP pointer to the perimeter of the polygon (meters). + * + * \e lats should be in the range [−90°, 90°]; \e lons should + * be in the range [−540°, 540°). + * + * Only simple polygons (which are not self-intersecting) are allowed. + * There's no need to "close" the polygon by repeating the first vertex. The + * area returned is signed with counter-clockwise traversal being treated as + * positive. + * + * Example, compute the area of Antarctica: + @code + double + lats[] = {-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7, + -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7}, + lons[] = {-74, -102, -102, -131, -163, 163, 172, 140, 113, + 88, 59, 25, -4, -14, -33, -46, -61}; + struct geod_geodesic g; + double A, P; + geod_init(&g, 6378137, 1/298.257223563); + geod_polygonarea(&g, lats, lons, (sizeof lats) / (sizeof lats[0]), &A, &P); + printf("%.0f %.2f\n", A, P); + @endcode + **********************************************************************/ + void geod_polygonarea(const struct geod_geodesic* g, + double lats[], double lons[], int n, + double* pA, double* pP); + + /** + * mask values for the \e caps argument to geod_lineinit(). + **********************************************************************/ + enum geod_mask { + GEOD_NONE = 0U, /**< Calculate nothing */ + GEOD_LATITUDE = 1U<<7 | 0U, /**< Calculate latitude */ + GEOD_LONGITUDE = 1U<<8 | 1U<<3, /**< Calculate longitude */ + GEOD_AZIMUTH = 1U<<9 | 0U, /**< Calculate azimuth */ + GEOD_DISTANCE = 1U<<10 | 1U<<0, /**< Calculate distance */ + GEOD_DISTANCE_IN = 1U<<11 | 1U<<0 | 1U<<1, /**< Allow distance as input */ + GEOD_REDUCEDLENGTH= 1U<<12 | 1U<<0 | 1U<<2, /**< Calculate reduced length */ + GEOD_GEODESICSCALE= 1U<<13 | 1U<<0 | 1U<<2, /**< Calculate geodesic scale */ + GEOD_AREA = 1U<<14 | 1U<<4, /**< Calculate reduced length */ + GEOD_ALL = 0x7F80U| 0x1FU /**< Calculate everything */ + }; + + /** + * flag values for the \e flags argument to geod_gendirect() and + * geod_genposition() + **********************************************************************/ + enum geod_flags { + GEOD_NOFLAGS = 0U, /**< No flags */ + GEOD_ARCMODE = 1U<<0, /**< Position given in terms of arc distance */ + GEOD_LONG_NOWRAP = 1U<<15 /**< Don't wrap longitude */ + }; + +#if defined(__cplusplus) +} +#endif + +#endif diff --git a/liblwgeom/liblwgeom_internal.h b/liblwgeom/liblwgeom_internal.h index fff65278c36..49cbf401421 100644 --- a/liblwgeom/liblwgeom_internal.h +++ b/liblwgeom/liblwgeom_internal.h @@ -32,15 +32,13 @@ #include "liblwgeom.h" +/* Define to enable pre-version 2.2 geodesic functions for geography types + (e.g. Vincenty for ST_Distance); otherwise use GeographicLib */ +#define USE_PRE22GEODESIC +#undef USE_PRE22GEODESIC /** -* PI -*/ -#define PI 3.1415926535897932384626433832795 - - -/** -* Floating point comparitors. +* Floating point comparators. */ #define FP_TOLERANCE 1e-12 #define FP_IS_ZERO(A) (fabs(A) <= FP_TOLERANCE) diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h index 74c976d01d7..ae75ca26061 100644 --- a/liblwgeom/lwgeodetic.h +++ b/liblwgeom/lwgeodetic.h @@ -57,8 +57,8 @@ typedef struct /** * Conversion functions */ -#define deg2rad(d) (PI * (d) / 180.0) -#define rad2deg(r) (180.0 * (r) / PI) +#define deg2rad(d) (M_PI * (d) / 180.0) +#define rad2deg(r) (180.0 * (r) / M_PI) /** * Ape a java function @@ -143,4 +143,4 @@ int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double * Maintain consistent units (radians?) throughout all calculations * Put an index pointer onto LWGEOM itself, and cache the indexed LWGEOM instead of a bare tree * only primitive objects should get a tree -*/ \ No newline at end of file +*/ diff --git a/liblwgeom/lwspheroid.c b/liblwgeom/lwspheroid.c index a16f9172120..f8b68232845 100644 --- a/liblwgeom/lwspheroid.c +++ b/liblwgeom/lwspheroid.c @@ -14,6 +14,9 @@ #include "lwgeodetic.h" #include "lwgeom_log.h" +/* GeographicLib */ +#include "geodesic.h" + /** * Initialize spheroid object based on major and minor axis */ @@ -26,6 +29,7 @@ void spheroid_init(SPHEROID *s, double a, double b) s->radius = (2.0 * a + b ) / 3.0; } +#ifdef USE_PRE22GEODESIC static double spheroid_mu2(double alpha, const SPHEROID *s) { double b2 = POW2(s->b); @@ -41,6 +45,113 @@ static double spheroid_big_b(double u2) { return (u2 / 1024.0) * (256.0 + u2 * (-128.0 + u2 * (74.0 - 47.0 * u2))); } +#endif /* def USE_PRE22GEODESIC */ + + +#ifndef USE_PRE22GEODESIC + +/** +* Computes the shortest distance along the surface of the spheroid +* between two points, using the inverse geodesic problem from +* GeographicLib (Karney 2013). +* +* @param a - location of first point +* @param b - location of second point +* @param s - spheroid to calculate on +* @return spheroidal distance between a and b in spheroid units +*/ +double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid) +{ + struct geod_geodesic gd; + geod_init(&gd, spheroid->a, spheroid->f); + double lat1 = a->lat * 180.0 / M_PI; + double lon1 = a->lon * 180.0 / M_PI; + double lat2 = b->lat * 180.0 / M_PI; + double lon2 = b->lon * 180.0 / M_PI; + double s12; /* return distance */ + geod_inverse(&gd, lat1, lon1, lat2, lon2, &s12, 0, 0); + return s12; +} + +/** +* Computes the forward azimuth of the geodesic joining two points on +* the spheroid, using the inverse geodesic problem (Karney 2013). +* +* @param r - location of first point +* @param s - location of second point +* @return azimuth of line joining r to s (but not reverse) +*/ +double spheroid_direction(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid) +{ + struct geod_geodesic gd; + geod_init(&gd, spheroid->a, spheroid->f); + double lat1 = a->lat * 180.0 / M_PI; + double lon1 = a->lon * 180.0 / M_PI; + double lat2 = b->lat * 180.0 / M_PI; + double lon2 = b->lon * 180.0 / M_PI; + double azi1; /* return azimuth */ + geod_inverse(&gd, lat1, lon1, lat2, lon2, 0, &azi1, 0); + return azi1 * M_PI / 180.0; +} + +/** +* Given a location, an azimuth and a distance, computes the location of +* the projected point. Using the direct geodesic problem from +* GeographicLib (Karney 2013). +* +* @param r - location of first point +* @param distance - distance in meters +* @param azimuth - azimuth in radians +* @return g - location of projected point +*/ +int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g) +{ + struct geod_geodesic gd; + geod_init(&gd, spheroid->a, spheroid->f); + double lat1 = r->lat * 180.0 / M_PI; + double lon1 = r->lon * 180.0 / M_PI; + double lat2, lon2; /* return projected position */ + geod_direct(&gd, lat1, lon1, azimuth * 180.0 / M_PI, distance, &lat2, &lon2, 0); + g->lat = lat2 * M_PI / 180.0; + g->lon = lon2 * M_PI / 180.0; + return LW_SUCCESS; +} + + +static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *spheroid) +{ + /* Return zero on non-sensical inputs */ + if ( ! pa || pa->npoints < 4 ) + return 0.0; + + struct geod_geodesic gd; + geod_init(&gd, spheroid->a, spheroid->f); + struct geod_polygon poly; + geod_polygon_init(&poly, 0); + int i; + double area; /* returned polygon area */ + POINT2D p; /* long/lat units are degrees */ + + /* Pass points from point array; don't close the linearring */ + for ( i = 0; i < pa->npoints - 1; i++ ) + { + getPoint2d_p(pa, i, &p); + geod_polygon_addpoint(&gd, &poly, p.y, p.x); + LWDEBUGF(4, "geod_polygon_addpoint %d: %.12g %.12g", i, p.y, p.x); + } + i = geod_polygon_compute(&gd, &poly, 0, 1, &area, 0); + if ( i != pa->npoints - 1 ) + { + lwerror("ptarray_area_spheroid: different number of points %d vs %d", + i, pa->npoints - 1); + } + LWDEBUGF(4, "geod_polygon_compute area: %.12g", area); + return fabs(area); +} + +/* Above use GeographicLib */ +#else /* ndef USE_PRE22GEODESIC */ +/* Below use pre-version 2.2 geodesic functions */ /** * Computes the shortest distance along the surface of the spheroid @@ -244,7 +355,7 @@ int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double { azimuth = azimuth + M_PI * 2.0; } - if (azimuth > (PI * 2.0)) + if (azimuth > (M_PI * 2.0)) { azimuth = azimuth - M_PI * 2.0; } @@ -503,6 +614,8 @@ static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *sphero } return fabs(area); } +#endif /* else USE_PRE22GEODESIC */ + /** * Calculate the area of an LWGEOM. Anything except POLYGON, MULTIPOLYGON * and GEOMETRYCOLLECTION return zero immediately. Multi's recurse, polygons diff --git a/liblwgeom/measures.c b/liblwgeom/measures.c index 26567a179ac..9b1efef2c82 100644 --- a/liblwgeom/measures.c +++ b/liblwgeom/measures.c @@ -2259,16 +2259,16 @@ azimuth_pt_pt(const POINT2D *A, const POINT2D *B, double *d) { if ( A->x == B->x ) { - if ( A->y < B->y ) *d=0.0; - else if ( A->y > B->y ) *d=M_PI; + if ( A->y < B->y ) *d = 0.0; + else if ( A->y > B->y ) *d = M_PI; /* 180 degrees */ else return 0; return 1; } if ( A->y == B->y ) { - if ( A->x < B->x ) *d=M_PI/2; - else if ( A->x > B->x ) *d=M_PI+(M_PI/2); + if ( A->x < B->x ) *d = M_PI_2; /* 90 degrees */ + else if ( A->x > B->x ) *d = M_PI + M_PI_2; /* 270 degrees */ else return 0; return 1; } @@ -2282,7 +2282,7 @@ azimuth_pt_pt(const POINT2D *A, const POINT2D *B, double *d) else /* ( A->y > B->y ) - equality case handled above */ { *d=atan(fabs(A->y - B->y) / fabs(A->x - B->x) ) - + (M_PI/2); + + M_PI_2; } } @@ -2296,7 +2296,7 @@ azimuth_pt_pt(const POINT2D *A, const POINT2D *B, double *d) else /* ( A->y < B->y ) - equality case handled above */ { *d=atan(fabs(A->y - B->y) / fabs(A->x - B->x) ) - + (M_PI+(M_PI/2)); + + (M_PI + M_PI_2); } } diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c index b848ead2afe..11cffd0977b 100644 --- a/postgis/geography_measurement.c +++ b/postgis/geography_measurement.c @@ -27,6 +27,14 @@ #include "geography_measurement_trees.h" /* For circ_tree caching */ #include "lwgeom_transform.h" /* For SRID functions */ +#ifdef USE_PRE22GEODESIC +/* round to 100 nm precision */ +#define INVMINDIST 1.0e9 +#else +/* round to 10 nm precision */ +#define INVMINDIST 1.0e8 +#endif + Datum geography_distance(PG_FUNCTION_ARGS); Datum geography_distance_uncached(PG_FUNCTION_ARGS); Datum geography_distance_tree(PG_FUNCTION_ARGS); @@ -162,6 +170,9 @@ Datum geography_distance(PG_FUNCTION_ARGS) PG_FREE_IF_COPY(g1, 0); PG_FREE_IF_COPY(g2, 1); + /* Knock off any funny business at the nanometer level, ticket #2168 */ + distance = round(distance * INVMINDIST) / INVMINDIST; + /* Something went wrong, negative return... should already be eloged, return NULL */ if ( distance < 0.0 ) { @@ -169,9 +180,6 @@ Datum geography_distance(PG_FUNCTION_ARGS) PG_RETURN_NULL(); } - /* Knock off any funny business at the micrometer level, ticket #2168 */ - distance = round(distance * 10e8) / 10e8; - PG_RETURN_FLOAT8(distance); } @@ -424,6 +432,7 @@ Datum geography_area(PG_FUNCTION_ARGS) else lwgeom_calculate_gbox_geodetic(lwgeom, &gbox); +#ifdef USE_PRE22GEODESIC /* Test for cases that are currently not handled by spheroid code */ if ( use_spheroid ) { @@ -434,6 +443,7 @@ Datum geography_area(PG_FUNCTION_ARGS) if ( gbox.zmax > 0.0 && gbox.zmin < 0.0 ) use_spheroid = LW_FALSE; } +#endif /* ifdef USE_PRE22GEODESIC */ /* User requests spherical calculation, turn our spheroid into a sphere */ if ( ! use_spheroid ) diff --git a/postgis/lwgeom_spheroid.c b/postgis/lwgeom_spheroid.c index 0a830330421..96504de0a9c 100644 --- a/postgis/lwgeom_spheroid.c +++ b/postgis/lwgeom_spheroid.c @@ -449,7 +449,7 @@ double distance_sphere_method(double lat1, double long1,double lat2,double long2 R = Geocent_a / (sqrt(1.0e0 - Geocent_e2 * sin2_lat)); /* 90 - lat1, but in radians */ - S = R * sin(M_PI/2.0-lat1) ; + S = R * sin(M_PI_2 - lat1) ; deltaX = long2 - long1; /* in rads */ deltaY = lat2 - lat1; /* in rads */