diff options
Diffstat (limited to '')
34 files changed, 763 insertions, 157 deletions
@@ -54,7 +54,9 @@ MAN0 = libtellurian.h.0 MAN3 = $(MAN3_FUNC) $(MAN3_CONST) MAN7 = libtellurian.7 -MAN3_FUNC = $(OBJ_PUBLIC:.o=.3) +MAN3_FUNC =\ + $(OBJ_PUBLIC:.o=.3)\ + libtellurian_coarse_distance_radius.3 MAN3_CONST =\ LIBTELLURIAN_EQUATORIAL_RADIUS.3\ @@ -1,7 +1,7 @@ NAME libtellurian - Geodesy library -SYNPOSIS +SYNOPSIS #include <libtellurian.h> Link with -ltellurian -lm. @@ -16,8 +16,8 @@ DESCRIPTION function name uses the suffix "_radians", in which case the function uses radians instead of degrees. All functions that input or output angles have a - version that uses degress and a version that uses + version that uses degrees and a version that uses radians. Unless otherwise specified, functions that - use degress convert the input and output to and from + use degrees convert the input and output to and from radians and call the version of the function that uses radians. @@ -1,5 +1,52 @@ -Finish libtellurian_end_point_radians Finish libtellurian_vincenty_inverse__ +Document the details of NaN returns in libtellurian_distance.3 +Document what happens to the azimuths when the starting point is a pole. -Add a method for creating a locally optimised spherical model of the earth, -and add method for calculating distances on this model +libtellurian_coarse_distance can be simplified to (d/2R)² +which would be an good alternative for simply sorting distances, +and it would still be easy to afterwards convert the values for +a selection of them into approximate distances. + +Add inverse function of libtellurian_effective_gravity_radians + +Add inverse function of libtellurian_elevated_gravity_radians + +And functions for converting between coordinate systems + Geodetic (h, φ, λ) : Angle between normal and equatorial plane + Inverse of Cartesian (see Cartesian for definition of N) + Iteration is required unless φ or h is known + λ = atan2(Y, X) + h = p / cos φ - N + φ = tan⁻¹(Zp⁻¹/(1 - e²N/(N + h))) + where p = √(X² + Y²) + https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_ECEF_to_geodetic_coordinates + When h=0 + λ = atan2(Y, X) + φ = tan⁻¹(Zp⁻¹/(1 - e²)) + where p = √(X² + Y²) + wouldn't it easier to convert to geocentric intermittently + Geocentric ϕ : Angle from centre of earth + ϕ = tan⁻¹((1 - e²) tan φ) + Geometric : Arc length = geodetic + Parametric/reduced β : spherical angle resulting in same distance from polar axis as + β = tan⁻¹((1 - f) tan φ) + Ellipsoidal-harmonic + TODO + Rectifying μ + μ = πm(φ) / 2m(½π), where m(u) = a(1 - e²) ∫{0→u} √(1 - e² sin² v)⁻³ dv + Authalic ξ + ξ = sin⁻¹(q(φ) / q(½π)), where q(u) = ((1 - e²) sin u)/(1 - e² sin u) + (1 - e²) e⁻¹ tanh⁻¹ (e sin u) + q(½π) = 1 + (1 - e²) e⁻¹ tanh⁻¹ e + Conformal χ + χ = tan⁻¹ (sinh [sinh⁻¹ tan φ - e tanh⁻¹ (e sin φ)]) + Isometric ψ + ψ = sinh⁻¹ tan φ - e than⁻¹ (e sin φ) + Astronomical Φ : angle between equatorial plane and the true vertical direction + the true vertical direction, is the direction of gravity which is affect + byu the centrifugal acceleration in addition to the gravitational acceleration. + Cartesian (geocentric Cartesian) + X = (N + h) cos φ cos λ + Y = (N + h) cos φ sin λ + Z = (Nb²/a² + h) sin φ + where + N = a²/√(a² cos² φ + b² sin² φ) @@ -14,6 +14,10 @@ #define degrees(RAD) (180.0 / (double)M_PI * (RAD)) #define haversin(X) (fma(-0.5, cos(X), 0.5)) +#define geodetic(GEOCENTRIC)\ + atan(tan(GEOCENTRIC) * (LIBTELLURIAN_EQUATORIAL_RADIUS * LIBTELLURIAN_EQUATORIAL_RADIUS /\ + (LIBTELLURIAN_POLAR_RADIUS * LIBTELLURIAN_POLAR_RADIUS))) + void libtellurian_vincenty_inverse__(double latitude1, double longitude1, double latitude2, double longitude2, double *distance_out, double *azimuth1_out, double *azimuth2_out); diff --git a/libtellurian.7 b/libtellurian.7 index 0933aa9..0ab50be 100644 --- a/libtellurian.7 +++ b/libtellurian.7 @@ -1,8 +1,8 @@ -.TH LIBTELLURIAN 7 libtellurian +.TH LIBTELLURIAN 7 LIBTELLURIAN .SH NAME libtellurian \- Geodesy library -.SH SYNPOSIS +.SH SYNOPSIS .nf #include <libtellurian.h> .fi @@ -23,9 +23,9 @@ and degrees for the angle unit, except when the function name uses the suffix \(dq_radians\(dq, in which case the function uses radians instead of degrees. All functions that input or output angles have a -version that uses degress and a version that uses +version that uses degrees and a version that uses radians. Unless otherwise specified, functions that -use degress convert the input and output to and from +use degrees convert the input and output to and from radians and call the version of the function that uses radians. .PP @@ -95,3 +95,16 @@ radius of curvature for some latitude. See .BR libtellurian.h (0) for a listing of available constants. +.SH NOTES +Users, and developers in particular, shall be aware that +this library uses GPS coordinates which are geocentric. +This means that the latitude is base one the angle from +the centre on the earth. There are a number of other +systems, most famously as a linear function of the +arc length, or more commonly in geodesy, the geodesic +latitude which is the enable between the local normal +from the surface and equatorial plane. Users should +take great care ensure the correct correct system (GPS) +is used, and developers should be careful designing +functions for the GPS coordinate system and be aware +the geodesic formulae usually use geodesic latitudes. diff --git a/libtellurian.h b/libtellurian.h index 609f61d..e96cfec 100644 --- a/libtellurian.h +++ b/libtellurian.h @@ -64,7 +64,7 @@ /** * The circumference, in meters, of a sphere inscribed in the Earth's - * spheroid and intersecting with it's pole (the circumference of + * spheroid and intersecting with its pole (the circumference of * a circle with Earth's polar radius) */ #define LIBTELLURIAN_POLAR_CIRCUMFERENCE 39940652.74224401 /* 2bπ */ @@ -123,15 +123,27 @@ * The geocentric gravitational constant, in cubic meters per square second * * This is the (universal) gravitational constant (6.67430e-11) multiplied - * by the mass of the Earth, however this value more reliable that the - * gravitation constant and the mass of th Earth, and should thus be used - * instead of multiplying the universial gravitational constant with the - * the msas of the Earth + * by the mass of the Earth; however, this value is more reliable than the + * gravitational constant and the mass of the Earth. It should therefore be used + * instead of multiplying the universal gravitational constant by the mass of + * the Earth. */ #define LIBTELLURIAN_GEOCENTRIC_GRAVITATIONAL_CONSTANT 3.986004418e14 /** + * The radius, in meters, used by the `libtellurian_coarse_distance` + * and `libtellurian_coarse_distance_radians` + * + * If you divide the return value of those function by this constant + * and then multiply the result with the local radius + * (`libtellurian_sea_level`), you will acquire or better local + * approximation of the actual distance. + */ +extern const double libtellurian_coarse_distance_radius; + + +/** * Calculate the distance of the nominal sea level (geocentric * radius), for some point on the Earth's surface, from the * centre of the Earth @@ -156,7 +168,7 @@ double libtellurian_sea_level_radians(double latitude); /** * Calculate the distance between two points on the Earth's surface * - * This function is gives an approximate value using + * This function gives an approximate value using * a sphere as a model for the Earth * * @param latitude1 GPS latitude coordinate for the first point, in degrees @@ -172,7 +184,7 @@ double libtellurian_coarse_distance(double latitude1, double longitude1, /** * Calculate the distance between two points on the Earth's surface * - * This function is gives an approximate value using + * This function gives an approximate value using * a sphere as a model for the Earth * * @param latitude1 GPS latitude coordinate for the first point, in radians @@ -188,25 +200,34 @@ double libtellurian_coarse_distance_radians(double latitude1, double longitude1, /** * Calculate the distance and azimuths between two points on the Earth's surface * - * This function is gives good approximate values + * This function gives good approximate values * using an oblate spheroid as a model for the Earth * - * @param latitude1 GPS latitude coordinate for the first point, in degrees - * @param longitude1 GPS longitude coordinate for the first point, in degrees - * @param latitude2 GPS latitude coordinate for the second point, in degrees - * @param longitude2 GPS longitude coordinate for the second point, in degrees + * @param latitude1 GPS latitude coordinate for the starting point, in degrees + * @param longitude1 GPS longitude coordinate for the starting point, in degrees + * @param latitude2 GPS latitude coordinate for the end point, in degrees + * @param longitude2 GPS longitude coordinate for the end point, in degrees * @param azimuth1_out Output parameter for the forward azimuth, in degrees, - * at the first point; or `NULL` + * at the starting point; or `NULL` * @param azimuth2_out Output parameter for the forward azimuth, in degrees, - * at the second point; or `NULL` + * at the end point; or `NULL` * @return Approximate distance between the two points * - * Calculating the the forward azimuths is not free, but it + * Calculating the forward azimuths is not free, but it * is cheap to compute them (especially the first one) when * most of the computations for the distance have been * performed. If you have no need for an azimuth you can set * the corresponding output parameter to `NULL`, and the * function will not compute it. + * + * If the two points are the same, the distance will be 0 + * and the azimuths will be NaN, denoting that any direction + * can be taken. If the two points are antipodal, the distance + * will be `LIBTELLURIAN_MERIDIONAL_CIRCUMFERENCE` and the + * azimuths will be NaN denoting that there are multiple directions + * to take, specifying either any (if the two points are the + * poles) or either north or south (and the two azimuths shall + * be the opposite of each other). */ LIBTELLURIAN_WUR__ double libtellurian_distance(double latitude1, double longitude1, @@ -216,25 +237,34 @@ double libtellurian_distance(double latitude1, double longitude1, /** * Calculate the distance and azimuths between two points on the Earth's surface * - * This function is gives good approximate values + * This function gives good approximate values * using an oblate spheroid as a model for the Earth * - * @param latitude1 GPS latitude coordinate for the first point, in radians - * @param longitude1 GPS longitude coordinate for the first point, in radians - * @param latitude2 GPS latitude coordinate for the second point, in radians - * @param longitude2 GPS longitude coordinate for the second point, in radians + * @param latitude1 GPS latitude coordinate for the starting point, in radians + * @param longitude1 GPS longitude coordinate for the starting point, in radians + * @param latitude2 GPS latitude coordinate for the end point, in radians + * @param longitude2 GPS longitude coordinate for the end point, in radians * @param azimuth1_out Output parameter for the forward azimuth, in radians, - * at the first point; or `NULL` + * at the starting point; or `NULL` * @param azimuth2_out Output parameter for the forward azimuth, in radians, - * at the second point; or `NULL` + * at the end point; or `NULL` * @return Approximate distance between the two points * - * Calculating the the forward azimuths is not free, but it + * Calculating the forward azimuths is not free, but it * is cheap to compute them (especially the first one) when * most of the computations for the distance have been * performed. If you have no need for an azimuth you can set * the corresponding output parameter to `NULL`, and the * function will not compute it. + * + * If the two points are the same, the distance will be 0 + * and the azimuths will be NaN, denoting that any direction + * can be taken. If the two points are antipodal, the distance + * will be `LIBTELLURIAN_MERIDIONAL_CIRCUMFERENCE` and the + * azimuths will be NaN denoting that there are multiple directions + * to take, specifying either any (if the two points are the + * poles) or either north or south (and the two azimuths shall + * be the opposite of each other). */ LIBTELLURIAN_WUR__ double libtellurian_distance_radians(double latitude1, double longitude1, @@ -244,21 +274,28 @@ double libtellurian_distance_radians(double latitude1, double longitude1, /** * Calculate the azimuths between two points on the Earth's surface * - * This function is gives good approximate values + * This function gives good approximate values * using an oblate spheroid as a model for the Earth * - * @param latitude1 GPS latitude coordinate for the first point, in degrees - * @param longitude1 GPS longitude coordinate for the first point, in degrees - * @param latitude2 GPS latitude coordinate for the second point, in degrees - * @param longitude2 GPS longitude coordinate for the second point, in degrees + * @param latitude1 GPS latitude coordinate for the starting point, in degrees + * @param longitude1 GPS longitude coordinate for the starting point, in degrees + * @param latitude2 GPS latitude coordinate for the end point, in degrees + * @param longitude2 GPS longitude coordinate for the end point, in degrees * @param azimuth1_out Output parameter for the forward azimuth, in degrees, - * at the first point; or `NULL` + * at the starting point; or `NULL` * @param azimuth2_out Output parameter for the forward azimuth, in degrees, - * at the second point; or `NULL` + * at the end point; or `NULL` * - * This function is identical to libtellurian_distance` + * This function is identical to `libtellurian_distance`. * except it does not compute the distance between the * points + * + * If the two points are the same, the azimuths will be NaN, + * denoting that any direction can be taken. If the two points + * are antipodal, the azimuths will be NaN denoting that there are + * multiple directions to take, specifying either any (if the + * two points are the poles) or either north or south (and + * the two azimuths shall be the opposite of each other). */ void libtellurian_azimuth(double latitude1, double longitude1, double latitude2, double longitude2, @@ -267,21 +304,28 @@ void libtellurian_azimuth(double latitude1, double longitude1, /** * Calculate the azimuths between two points on the Earth's surface * - * This function is gives good approximate values + * This function gives good approximate values * using an oblate spheroid as a model for the Earth * - * @param latitude1 GPS latitude coordinate for the first point, in radians - * @param longitude1 GPS longitude coordinate for the first point, in radians - * @param latitude2 GPS latitude coordinate for the second point, in radians - * @param longitude2 GPS longitude coordinate for the second point, in radians + * @param latitude1 GPS latitude coordinate for the starting point, in radians + * @param longitude1 GPS longitude coordinate for the starting point, in radians + * @param latitude2 GPS latitude coordinate for the end point, in radians + * @param longitude2 GPS longitude coordinate for the end point, in radians * @param azimuth1_out Output parameter for the forward azimuth, in radians, - * at the first point; or `NULL` + * at the starting point; or `NULL` * @param azimuth2_out Output parameter for the forward azimuth, in radians, - * at the second point; or `NULL` + * at the end point; or `NULL` * * - * This function is identical to libtellurian_distance_radians` + * This function is identical to `libtellurian_distance_radians`. * except it does not compute the distance between the points + * + * If the two points are the same, the azimuths will be NaN, + * denoting that any direction can be taken. If the two points + * are antipodal, the azimuths will be NaN denoting that there are + * multiple directions to take, specifying either any (if the + * two points are the poles) or either north or south (and + * the two azimuths shall be the opposite of each other). */ void libtellurian_azimuth_radians(double latitude1, double longitude1, double latitude2, double longitude2, @@ -299,8 +343,8 @@ void libtellurian_azimuth_radians(double latitude1, double longitude1, * coordinate, in degrees; or `NULL` * @param longitude2_out Output parameter for the end point's GPS longitude * coordinate, in degrees; or `NULL` - * @param azimuth2_out Output parameter for the direction from the end point - * to the starting point, in degrees; or `NULL` + * @param azimuth2_out Output parameter for the forward azimuth at the endpoint, + * in degrees; or `NULL` */ void libtellurian_end_point(double latitude1, double longitude1, double azimuth1, double distance, double *latitude2_out, double *longitude2_out, double *azimuth2_out); @@ -317,15 +361,15 @@ void libtellurian_end_point(double latitude1, double longitude1, double azimuth1 * coordinate, in radians; or `NULL` * @param longitude2_out Output parameter for the end point's GPS longitude * coordinate, in radians; or `NULL` - * @param azimuth2_out Output parameter for the direction from the end point - * to the starting point, in radians; or `NULL` + * @param azimuth2_out Output parameter for the forward azimuth at the endpoint, + * in radians; or `NULL` */ void libtellurian_end_point_radians(double latitude1, double longitude1, double azimuth1, double distance, double *latitude2_out, double *longitude2_out, double *azimuth2_out); /** * Calculate the normal gravity for some point on - * the Earth's surface, where the Earth's is model + * the Earth's surface, where the Earth is modelled * as an oblate spheroid * * @param latitude GPS latitude coordinate, in degrees @@ -336,7 +380,7 @@ double libtellurian_normal_gravity(double latitude); /** * Calculate the normal gravity for some point on - * the Earth's surface, where the Earth's is model + * the Earth's surface, where the Earth is modelled * as an oblate spheroid * * @param latitude GPS latitude coordinate, in radians @@ -373,7 +417,7 @@ double libtellurian_effective_gravity_radians(double gravity, double latitude); /** * Calculate the gravity adjusted for the elevation - * above the altitude where the gravity is measure + * above the altitude where the gravity is measured * * Altitudes above circa 100000 meters is out of range * for this function (that would be in outer space) @@ -388,7 +432,7 @@ double libtellurian_elevated_gravity(double gravity, double latitude, double alt /** * Calculate the gravity adjusted for the elevation - * above the altitude where the gravity is measure + * above the altitude where the gravity is measured * * Altitudes above circa 100000 meters is out of range * for this function (that would be in outer space) diff --git a/libtellurian.h.0 b/libtellurian.h.0 index 472d35c..5907df9 100644 --- a/libtellurian.h.0 +++ b/libtellurian.h.0 @@ -1,8 +1,8 @@ -.TH LIBTELLURIAN.H 0 libtellurian +.TH LIBTELLURIAN.H 0 LIBTELLURIAN .SH NAME libtellurian.h \- Geodesy library header -.SH SYNPOSIS +.SH SYNOPSIS .nf #include <libtellurian.h> diff --git a/libtellurian_azimuth.c b/libtellurian_azimuth.c index 4bfc743..c39596a 100644 --- a/libtellurian_azimuth.c +++ b/libtellurian_azimuth.c @@ -21,5 +21,40 @@ libtellurian_azimuth(double latitude1, double longitude1, #else -TODO_TEST + + +static int +check(double da, double db, double dc, double dd, + double ra, double rb, double rc, double rd) +{ + double az1_ref = 1, az2_ref = 2, az1 = 3, az2 = 4, az1b = 5, az2b = 6; + + libtellurian_azimuth_radians(ra, rb, rc, rd, &az1_ref, &az2_ref); + az1_ref *= 180 / M_PI; + az2_ref *= 180 / M_PI; + + libtellurian_azimuth(da, db, dc, dd, &az1, &az2); + libtellurian_azimuth(da, db, dc, dd, &az1b, NULL); + libtellurian_azimuth(da, db, dc, dd, NULL, &az2b); + libtellurian_azimuth(da, db, dc, dd, NULL, NULL); + return az1 == az1_ref && az2 == az2_ref && az1 == az1b && az2 == az2b; +} + +int +main(void) +{ + double a1 = 1, a2 = 2; + + ASSERT(check(30, 45, 60, 90, D30, D45, D60, D90)); + ASSERT(check(30, 45, 60, 45, D30, D45, D60, D45)); + ASSERT(check(45, 30, 60, 45, D45, D30, D60, D45)); + ASSERT(check(45, 30, 30, 45, D45, D30, D30, D45)); + + libtellurian_azimuth(0, 0, 0, 90, &a1, &a2); + ASSERT(a1 == a2); + ASSERT(a1 == 90); + return 0; +} + + #endif diff --git a/libtellurian_azimuth_radians.c b/libtellurian_azimuth_radians.c index 18c954e..469eb45 100644 --- a/libtellurian_azimuth_radians.c +++ b/libtellurian_azimuth_radians.c @@ -17,9 +17,16 @@ libtellurian_azimuth_radians(double latitude1, double longitude1, #else -#if 1 -TODO_TEST -#else +static int +eq(double a, double b) +{ + if (isfinite(a) && isfinite(b)) + return a == b; + if (isnan(a) && isnan(b)) + return 1; + printf("Unexpected floating-point class\n"); + return 0; +} static int check(double a, double b, double c, double d) @@ -34,7 +41,7 @@ check(double a, double b, double c, double d) libtellurian_azimuth_radians(a, b, c, d, &az1b, NULL); libtellurian_azimuth_radians(a, b, c, d, NULL, &az2b); libtellurian_azimuth_radians(a, b, c, d, NULL, NULL); - return az1 == az1_ref && az2 == az2_ref && az1 == az1b && az2 == az2b; + return eq(az1, az1_ref) && eq(az2, az2_ref) && eq(az1, az1b) && eq(az2, az2b); } int @@ -49,7 +56,5 @@ main(void) return 0; } -#endif - #endif diff --git a/libtellurian_azimuthal_radius.3 b/libtellurian_azimuthal_radius.3 index a1be8aa..81c97fe 100644 --- a/libtellurian_azimuthal_radius.3 +++ b/libtellurian_azimuthal_radius.3 @@ -1,8 +1,8 @@ -.TH LIBTELLURIAN_AZIMUTHAL_RADIUS 3 libtellurian +.TH LIBTELLURIAN_AZIMUTHAL_RADIUS 3 LIBTELLURIAN .SH NAME -libtellurian_azimuthal_radius \- Calculate a azimuthal radius of curvature +libtellurian_azimuthal_radius \- Calculate an azimuthal radius of curvature -.SH SYNPOSIS +.SH SYNOPSIS .nf #include <libtellurian.h> diff --git a/libtellurian_azimuthal_radius_radians.c b/libtellurian_azimuthal_radius_radians.c index d58a903..e9433e3 100644 --- a/libtellurian_azimuthal_radius_radians.c +++ b/libtellurian_azimuthal_radius_radians.c @@ -6,8 +6,9 @@ double libtellurian_azimuthal_radius_radians(double latitude, double azimuth) { - double m = libtellurian_meridian_radius_radians(latitude); - double n = libtellurian_transverse_radius_radians(latitude); + double lat = geodetic(latitude); + double m = libtellurian_meridian_radius_radians(lat); + double n = libtellurian_transverse_radius_radians(lat); double c2 = cos(azimuth) * cos(azimuth); double s2 = sin(azimuth) * sin(azimuth); double x = c2 / m; diff --git a/libtellurian_coarse_distance.3 b/libtellurian_coarse_distance.3 index 9a3973b..c45f242 100644 --- a/libtellurian_coarse_distance.3 +++ b/libtellurian_coarse_distance.3 @@ -1,8 +1,8 @@ -.TH LIBTELLURIAN_COARSE_DISTANCE 3 libtellurian +.TH LIBTELLURIAN_COARSE_DISTANCE 3 LIBTELLURIAN .SH NAME libtellurian_coarse_distance \- Calculate distance between two locations -.SH SYNPOSIS +.SH SYNOPSIS .nf #include <libtellurian.h> @@ -11,6 +11,8 @@ double libtellurian_coarse_distance(double \fIlatitude1\fP, double \fIlongitude1 double libtellurian_coarse_distance_radians(double \fIlatitude1\fP, double \fIlongitude1\fP, double \fIlatitude2\fP, double \fIlongitude2\fP); + +extern const double libtellurian_coarse_distance_radius; .fi .PP Link with @@ -42,6 +44,17 @@ function except that and .I longitude2 shall be specified in radians. +.PP +The constant variable +.B libtellurian_coarse_distance_radius +is set to the radius used in the model of the +Earth that these functions use. If you want a +better local approximation of distances, you +can multiple the return values by the local radius +(you can use the +.BR libtellurian_sea_level (3) +function) divided by +.BR libtellurian_coarse_distance_radius . .SH RETURN VALUE The diff --git a/libtellurian_coarse_distance_radians.c b/libtellurian_coarse_distance_radians.c index f22ad28..65f1cac 100644 --- a/libtellurian_coarse_distance_radians.c +++ b/libtellurian_coarse_distance_radians.c @@ -6,6 +6,9 @@ #ifndef TEST +const double libtellurian_coarse_distance_radius = R; + + double libtellurian_coarse_distance_radians(double latitude1, double longitude1, double latitude2, double longitude2) @@ -48,6 +51,7 @@ main(void) ASSERT(check(-D45, -D45, D45, D45, M_PI / 1.5)); ASSERT(check(D45, -D45, D45, D45, M_PI / 3.0)); ASSERT(check(-D45, D45, D45, D45, M_PI / 2.0)); + ASSERT(libtellurian_coarse_distance_radius == R); return 0; } diff --git a/libtellurian_coarse_distance_radius.3 b/libtellurian_coarse_distance_radius.3 new file mode 120000 index 0000000..48e4642 --- /dev/null +++ b/libtellurian_coarse_distance_radius.3 @@ -0,0 +1 @@ +libtellurian_coarse_distance.3
\ No newline at end of file diff --git a/libtellurian_distance.3 b/libtellurian_distance.3 index 49181b3..f6f9135 100644 --- a/libtellurian_distance.3 +++ b/libtellurian_distance.3 @@ -1,8 +1,8 @@ -.TH LIBTELLURIAN_DISTANCE 3 libtellurian +.TH LIBTELLURIAN_DISTANCE 3 LIBTELLURIAN .SH NAME libtellurian_distance \- Calculate distance between two locations -.SH SYNPOSIS +.SH SYNOPSIS .nf #include <libtellurian.h> @@ -51,7 +51,7 @@ unless is .IR NULL , and the forward azimuth for the second point -will be calculated and stored, in degress, in +will be calculated and stored, in degrees, in .I *azimuth2_out unless .I azimuth2_out @@ -104,6 +104,44 @@ functions do not return any value. .SH ERRORS None. +.SH APPLICATION USAGE +The application should use +.BR libtellurian_azimuth () +or +.BR libtellurian_azimuth_radians () +whenever the distance will not be used by the +application, likewise, if +.I *azimuth1_out +will not be used +.I azimuth1_out +should be set to +.I NULL +and if +.I *azimuth2_out +will not be used +.I azimuth2_out +should be set to +.IR NULL . +Although this currently has almost no impact on +the application, future versions of the library +may select different methods for performing the +calculations depending on what output is enabled. + +.SH NOTES +The (forward) azimuths are defined as the +direction you travelling when you travel +along a great ellipsoid (which is the +shortest distance) from +.RI ( latitude1 ", " longitude1 ) +in the direction of +.RI ( latitude2 ", " longitude2 ). +Hence, +.I *azimuth2_out +will be the direction you would continue +along the great ellipsoid, \m[red]not\m[] +the direction for returning to +.RI ( latitude1 ", " longitude1 ). + .SH SEE ALSO .BR libtellurian (7), .BR libtellurian_coarse_distance (3) diff --git a/libtellurian_distance.c b/libtellurian_distance.c index 250dddc..dec0edd 100644 --- a/libtellurian_distance.c +++ b/libtellurian_distance.c @@ -23,5 +23,48 @@ libtellurian_distance(double latitude1, double longitude1, #else -TODO_TEST + + +static int +check(double da, double db, double dc, double dd, + double ra, double rb, double rc, double rd) +{ + double az1_ref = 1, az2_ref = 2, az1 = 3, az2 = 4, az1b = 5, az2b = 6; + double s_ref, s1, s2, s3, s4; + + s_ref = libtellurian_distance_radians(ra, rb, rc, rd, &az1_ref, &az2_ref); + az1_ref *= 180 / M_PI; + az2_ref *= 180 / M_PI; + + s1 = libtellurian_distance(da, db, dc, dd, &az1, &az2); + s2 = libtellurian_distance(da, db, dc, dd, &az1b, NULL); + s3 = libtellurian_distance(da, db, dc, dd, NULL, &az2b); + s4 = libtellurian_distance(da, db, dc, dd, NULL, NULL); + ASSERT(s1 == s_ref); + ASSERT(s2 == s_ref); + ASSERT(s3 == s_ref); + ASSERT(s4 == s_ref); + return az1 == az1_ref && az2 == az2_ref && az1 == az1b && az2 == az2b; +} + +int +main(void) +{ + double a1 = 1, a2 = 2, s; + + ASSERT(check(30, 45, 60, 90, D30, D45, D60, D90)); + ASSERT(check(30, 45, 60, 45, D30, D45, D60, D45)); + ASSERT(check(45, 30, 60, 45, D45, D30, D60, D45)); + ASSERT(check(45, 30, 30, 45, D45, D30, D30, D45)); + + s = libtellurian_distance(0, 0, 0, 90, &a1, &a2); + fprintf(stderr, "%lg\n", s); + fprintf(stderr, "%lg\n", s / LIBTELLURIAN_MERIDIONAL_CIRCUMFERENCE); + (void) s; + ASSERT(a1 == a2); + ASSERT(a1 == 90); + return 0; +} + + #endif diff --git a/libtellurian_distance_radians.c b/libtellurian_distance_radians.c index d7c5547..7401edd 100644 --- a/libtellurian_distance_radians.c +++ b/libtellurian_distance_radians.c @@ -17,5 +17,34 @@ libtellurian_distance_radians(double latitude1, double longitude1, #else -TODO_TEST + + +int +main(void) +{ + double rs, ra1, ra2; + double s, a1, a2; + + libtellurian_vincenty_inverse__(1, 2, 3, 4, &rs, &ra1, &ra2); + + s = libtellurian_distance_radians(1, 2, 3, 4, &a1, &a2); + ASSERT(s == rs); + ASSERT(a1 == ra1); + ASSERT(a2 == ra2); + + s = libtellurian_distance_radians(1, 2, 3, 4, &a1, NULL); + ASSERT(s == rs); + ASSERT(a1 == ra1); + + s = libtellurian_distance_radians(1, 2, 3, 4, NULL, &a2); + ASSERT(s == rs); + ASSERT(a2 == ra2); + + s = libtellurian_distance_radians(1, 2, 3, 4, NULL, NULL); + ASSERT(s == rs); + + return 0; +} + + #endif diff --git a/libtellurian_effective_gravity.3 b/libtellurian_effective_gravity.3 index d462964..43299a6 100644 --- a/libtellurian_effective_gravity.3 +++ b/libtellurian_effective_gravity.3 @@ -1,8 +1,8 @@ -.TH LIBTELLURIAN_EFFECTIVE_GRAVITY 3 libtellurian +.TH LIBTELLURIAN_EFFECTIVE_GRAVITY 3 LIBTELLURIAN .SH NAME libtellurian_effective_gravity \- Calculate effective gravity -.SH SYNPOSIS +.SH SYNOPSIS .nf #include <libtellurian.h> diff --git a/libtellurian_effective_gravity_radians.c b/libtellurian_effective_gravity_radians.c index 54e2b8e..cb23246 100644 --- a/libtellurian_effective_gravity_radians.c +++ b/libtellurian_effective_gravity_radians.c @@ -9,7 +9,8 @@ double libtellurian_effective_gravity_radians(double gravity, double latitude) { - return fma(-K / gravity, cos(latitude) * cos(latitude), gravity); + double lat = geodetic(latitude); /* actual one is geodetic and one is geographical, there is no actual difference */ + return fma(-K / gravity, cos(lat) * cos(lat), gravity); } diff --git a/libtellurian_elevated_gravity.3 b/libtellurian_elevated_gravity.3 index fb728b9..4b4dcd0 100644 --- a/libtellurian_elevated_gravity.3 +++ b/libtellurian_elevated_gravity.3 @@ -1,8 +1,8 @@ -.TH LIBTELLURIAN_ELEVATED_GRAVITY 3 libtellurian +.TH LIBTELLURIAN_ELEVATED_GRAVITY 3 LIBTELLURIAN .SH NAME libtellurian_elevated_gravity \- Calculate gravity at some altitude -.SH SYNPOSIS +.SH SYNOPSIS .nf #include <libtellurian.h> diff --git a/libtellurian_elevated_gravity_radians.c b/libtellurian_elevated_gravity_radians.c index 1c97f1a..edf9678 100644 --- a/libtellurian_elevated_gravity_radians.c +++ b/libtellurian_elevated_gravity_radians.c @@ -6,6 +6,7 @@ double libtellurian_elevated_gravity_radians(double gravity, double latitude, double altitude) { + double lat = geodetic(latitude); double a = LIBTELLURIAN_EQUATORIAL_RADIUS; double b = LIBTELLURIAN_POLAR_RADIUS; double omega = LIBTELLURIAN_ANGULAR_VELOCITY; @@ -15,7 +16,7 @@ libtellurian_elevated_gravity_radians(double gravity, double latitude, double al double neg_k1 = fma(-2.0, f + m, -2.0) / a; double k2 = 4.0 * f / a; double neg_k3 = -3.0 / (a * a); - double sin2_phi = sin(latitude) * sin(latitude); + double sin2_phi = sin(lat) * sin(lat); double s = fma(neg_k3, altitude, fma(k2, sin2_phi, neg_k1)); return fma(s * altitude, gravity, gravity); } diff --git a/libtellurian_end_point.3 b/libtellurian_end_point.3 index 15fe232..d657893 100644 --- a/libtellurian_end_point.3 +++ b/libtellurian_end_point.3 @@ -1,19 +1,19 @@ -.TH LIBTELLURIAN_END_POINT 3 libtellurian +.TH LIBTELLURIAN_END_POINT 3 LIBTELLURIAN .SH NAME libtellurian_end_point \- Calculate travel end-point -.SH SYNPOSIS +.SH SYNOPSIS .nf #include <libtellurian.h> void libtellurian_end_point(double \fIlatitude1\fP, double \fIlongitude1\fP, double \fIazimuth1\fP, double \fIdistance\fP, - double \fIlatitude2_out\fP, double \fIlongitude2_out\fP, + double *\fIlatitude2_out\fP, double *\fIlongitude2_out\fP, double *\fIazimuth2_out\fP); void libtellurian_end_point_radians(double \fIlatitude1\fP, double \fIlongitude1\fP, double \fIazimuth1\fP, double \fIdistance\fP, - double \fIlatitude2_out\fP, double \fIlongitude2_out\fP, + double *\fIlatitude2_out\fP, double *\fIlongitude2_out\fP, double *\fIazimuth2_out\fP); .fi .PP @@ -38,7 +38,10 @@ The location the traveller will end at will be stored at .RI ( *latitude2_out ", " *longitude2_out ). along the ellipsoid. The azimuth from this -point to the starting point will be stored in +point to continue the on the extended, circular +path (eventually returning to the starting +pointer after having traveled the azimuthal +circumference of the Earth) in .IR *azimuth2_out . .PP However @@ -58,14 +61,14 @@ Likewise .I azimuth1 shall be in degrees, and .I *azimuth2_out -will be in degress. +will be in degrees. .PP The .BR libtellurian_end_point_radians () function is identical to the .BR libtellurian_end_point () function except that radians are -used instead of degress for +used instead of degrees for .IR latitude1 , .IR longitude1 , .IR azimuth1 , @@ -80,5 +83,37 @@ None. .SH ERRORS None. +.SH APPLICATION USAGE +If +.I *azimuth2_out +will not be used by the application +.I azimuth2_out +should be set to +.I NULL +and likwise for +.I latitude2_out +and +.IR longitude2_out . +Although this currently has almost no impact on +the application, future versions of the library +may select different methods for performing the +calculations depending on what output is enabled. + +.SH NOTES +The (forward) azimuths are defined as the +direction you travelling when you travel +along a great ellipsoid (which is the +shortest distance) from +.RI ( latitude1 ", " longitude1 ) +in the direction specified by +.I azimuth1 +at this point. Hence, +.I *azimuth2_out +will be the direction you would continue +along the great ellipsoid, at +.RI ( *latitude2_out ", " *longitude2_out ), +\m[red]not\m[] the direction for returning to +.RI ( latitude1 ", " longitude1 ). + .SH SEE ALSO .BR libtellurian (7) diff --git a/libtellurian_end_point.c b/libtellurian_end_point.c index 57af198..df0046e 100644 --- a/libtellurian_end_point.c +++ b/libtellurian_end_point.c @@ -22,5 +22,76 @@ libtellurian_end_point(double latitude1, double longitude1, double azimuth1, dou #else -TODO_TEST + + +# define PM LIBTELLURIAN_MERIDIONAL_CIRCUMFERENCE +# define PE LIBTELLURIAN_EQUATORIAL_CIRCUMFERENCE + +static int +check(double dlat, double dlon, double daz, double rlat, double rlon, double raz, double s) +{ + double lat1, lat2, lon1, lon2, az1, az2; + libtellurian_end_point_radians(rlat, rlon, raz, s, &lat1, &lon1, &az1); + lat1 *= 180 / M_PI; + lon1 *= 180 / M_PI; + az1 *= 180 / M_PI; + libtellurian_end_point(dlat, dlon, daz, s, &lat2, &lon2, &az2); + if (lat2 != lat1 || lon2 != lon1 || az2 != az1) + return 0; + libtellurian_end_point(dlat, dlon, daz, s, &lat2, NULL, NULL); + libtellurian_end_point(dlat, dlon, daz, s, NULL, &lon2, NULL); + libtellurian_end_point(dlat, dlon, daz, s, NULL, NULL, &az2); + if (lat2 != lat1 || lon2 != lon1 || az2 != az1) + return 0; + libtellurian_end_point(dlat, dlon, daz, s, NULL, &lon2, &az2); + if (lon2 != lon1 || az2 != az1) + return 0; + libtellurian_end_point(dlat, dlon, daz, s, &lat2, NULL, &az2); + if (lat2 != lat1 || az2 != az1) + return 0; + libtellurian_end_point(dlat, dlon, daz, s, &lat2, &lon2, NULL); + if (lat2 != lat1 || lon2 != lon1) + return 0; + libtellurian_end_point(dlat, dlon, daz, s, NULL, NULL, NULL); + return 1; +} + +int +main(void) +{ + ASSERT(check(0, 0, 0, 0, 0, 0, 0)); + ASSERT(check(0, 0, 0, 0, 0, 0, PM / 4)); + ASSERT(check(0, 0, 0, 0, 0, 0, PM / 2)); + ASSERT(check(0, 0, 0, 0, 0, 0, PM)); + ASSERT(check(0, 0, 0, 0, 0, 0, 2 * PM)); + ASSERT(check(0, 0, 0, 0, 0, 0, 3 * PM / 2)); + ASSERT(check(0, 0, 90, 0, 0, D90, PE)); + ASSERT(check(0, 0, 90, 0, 0, D90, PE / 8)); + ASSERT(check(0, 0, 90, 0, 0, D90, PE / 4)); + ASSERT(check(0, 0, 90, 0, 0, D90, PE / 2)); + + ASSERT(check(0, 0, 0, 0, 0, 0, 1000.0)); + ASSERT(check(0, 0, 30, 0, 0, D30, 1000.0)); + ASSERT(check(0, 0, 45, 0, 0, D45, 1000.0)); + ASSERT(check(0, 0, 60, 0, 0, D60, 1000.0)); + ASSERT(check(0, 0, 90, 0, 0, D90, 1000.0)); + ASSERT(check(0, 0, 180, 0, 0, D180, 1000.0)); + + ASSERT(check(0, 30, 0, 0, D30, 0, 1000.0)); + ASSERT(check(0, 45, 30, 0, D45, D30, 1000.0)); + ASSERT(check(0, 60, 45, 0, D60, D45, 1000.0)); + ASSERT(check(0, 90, 60, 0, D90, D60, 1000.0)); + ASSERT(check(0, 180, 90, 0, D180, D90, 1000.0)); + ASSERT(check(0, 45, 180, 0, D45, D180, 1000.0)); + + ASSERT(check(45, 30, 0, D45, D30, 0, 1000.0)); + ASSERT(check(60, 45, 30, D60, D45, D30, 1000.0)); + ASSERT(check(90, 60, 45, D90, D60, D45, 1000.0)); + ASSERT(check(-45, 90, 60, -D45, D90, D60, 1000.0)); + ASSERT(check(-60, 180, 90, -D60, D180, D90, 1000.0)); + ASSERT(check(-90, 45, 180, -D90, D45, D180, 1000.0)); + return 0; +} + + #endif diff --git a/libtellurian_end_point_radians.c b/libtellurian_end_point_radians.c index c512b59..ff00f5c 100644 --- a/libtellurian_end_point_radians.c +++ b/libtellurian_end_point_radians.c @@ -3,8 +3,6 @@ #ifndef TEST -/* TODO complete the implementation */ - void libtellurian_end_point_radians(double latitude1, double longitude1, double azimuth1, double distance, double *latitude2_out, double *longitude2_out, double *azimuth2_out) @@ -12,6 +10,7 @@ libtellurian_end_point_radians(double latitude1, double longitude1, double azimu /* * Vincenty's formula for the "direct problem" */ + double t, sigma, cos_2sigma_m, sin_sigma, cos_sigma; double x, x0, y, C, L, v, lambda; @@ -38,23 +37,26 @@ libtellurian_end_point_radians(double latitude1, double longitude1, double azimu double B = fma(fma(fma(-47.0, uu, 74.0), uu, -128.0), uu, 256.0) * (uu / 1024.0); double sigma_0 = distance / (b * A); + double sigma_prev; + int max_interations = 20; sigma = sigma_0; - /* { */ + do { + sigma_prev = sigma; - cos_2sigma_m = cos(fma(2.0, sigma1, sigma)); - sin_sigma = sin(sigma); - cos_sigma = cos(sigma); + cos_2sigma_m = cos(fma(2.0, sigma1, sigma)); + sin_sigma = sin(sigma); + cos_sigma = cos(sigma); - v = cos_sigma * fma(2 * cos_2sigma_m, cos_2sigma_m, -1.0); - t = fma(2.0 * cos_2sigma_m, 2.0 * cos_2sigma_m, -3.0); - t *= fma(2.0 * sin_sigma, 2.0 * sin_sigma, -3.0); - t = fma(B * cos_2sigma_m / -6.0, t, v); - t = fma(0.25 * B, t, cos_2sigma_m); - sigma = fma(B * sin_sigma, t, sigma_0); + v = cos_sigma * fma(2 * cos_2sigma_m, cos_2sigma_m, -1.0); + t = fma(2.0 * cos_2sigma_m, 2.0 * cos_2sigma_m, -3.0); + t *= fma(2.0 * sin_sigma, 2.0 * sin_sigma, -3.0); + t = fma(B * cos_2sigma_m / -6.0, t, v); + t = fma(0.25 * B, t, cos_2sigma_m); + sigma = fma(B * sin_sigma, t, sigma_0); - /* } repeat until sigma converges */ + } while (fabs(sigma - sigma_prev) > 1e-14 && --max_interations); if (latitude2_out || azimuth2_out) { x0 = fma(cos_sigma * cos_alpha1, cos_u1, -sin_u1 * sin_alpha); @@ -82,5 +84,137 @@ libtellurian_end_point_radians(double latitude1, double longitude1, double azimu #else -TODO_TEST + + +# define PM LIBTELLURIAN_MERIDIONAL_CIRCUMFERENCE +# define PE LIBTELLURIAN_EQUATORIAL_CIRCUMFERENCE + +static int +approx(double a, double b) +{ + return fabs(a - b) <= 1e-8; +} + +static double +normal_angle(double u) +{ + u = fmod(u + 4.0 * M_PI, 2.0 * M_PI); + if (u > M_PI) + u -= 2.0 * M_PI; + return u; +} + +static void +end_point(double a, double b, double c, double d, + double *lat_out, double *lon_out, double *az_out) +{ + libtellurian_end_point_radians(a, b, c, d, lat_out, lon_out, az_out); + if (lat_out) { + *lat_out = normal_angle(*lat_out); + if (*lat_out > M_PI / 2) { + *lat_out = M_PI - *lat_out; + if (*lon_out) + *lon_out += M_PI; + } else if (*lat_out < -M_PI / 2) { + *lat_out = -M_PI - *lat_out; + if (*lon_out) + *lon_out += M_PI; + } + } + if (lon_out) + *lon_out = normal_angle(*lon_out); + if (az_out) + *az_out = normal_angle(*az_out); +} + +int +main(void) +{ + double lat, lon, az; + + end_point(0, 0, 0, 0, &lat, &lon, &az); + ASSERT(lat == 0); + ASSERT(lon == 0); + ASSERT(az == 0); + + end_point(0, 0, 3, 0, &lat, &lon, &az); + ASSERT(lat == 0); + ASSERT(lon == 0); + ASSERT(az == 3); + + end_point(0, 0, 6, 0, &lat, &lon, &az); + ASSERT(lat == 0); + ASSERT(lon == 0); + ASSERT(az == normal_angle(6)); + + end_point(0, 0, 0, PM / 4, &lat, &lon, &az); + ASSERT(approx(lat, D90)); + ASSERT(approx(lon, 0)); + + end_point(0, 0, 0, PM / 4, NULL, &lon, &az); + ASSERT(approx(lon, 0)); + + end_point(0, 0, 0, PM / 2, &lat, &lon, &az); + ASSERT(approx(lat, 0)); + ASSERT(approx(lon, D180)); + ASSERT(approx(az, D180)); + + end_point(0, 0, 0, PM / 2, &lat, NULL, &az); + ASSERT(approx(lat, 0)); + ASSERT(approx(az, D180)); + + end_point(0, 0, 0, PM / 2, &lat, &lon, NULL); + ASSERT(approx(lat, 0)); + ASSERT(approx(lon, D180)); + + end_point(0, 0, 0, PM / 2, &lat, NULL, NULL); + ASSERT(approx(lat, 0)); + + end_point(0, 0, 0, PM / 2, NULL, NULL, NULL); + + end_point(0, 0, 0, PM / 2, NULL, NULL, &az); + ASSERT(approx(az, D180)); + + end_point(0, 0, 0, PM / 2, NULL, &lon, NULL); + ASSERT(approx(lon, D180)); + + end_point(0, 0, 0, PM, &lat, &lon, &az); + ASSERT(approx(lat, 0)); + ASSERT(approx(lon, 0)); + ASSERT(approx(az, 0)); + + end_point(0, 0, 0, 2 * PM, &lat, &lon, &az); + ASSERT(approx(lat, 0)); + ASSERT(approx(lon, 0)); + ASSERT(approx(az, 0)); + + end_point(0, 0, 0, 3 * PM / 2, &lat, &lon, &az); + ASSERT(approx(lat, 0)); + ASSERT(approx(lon, D180)); + ASSERT(approx(az, D180)); + + end_point(0, 0, D90, PE, &lat, &lon, &az); + ASSERT(approx(lat, 0)); + ASSERT(approx(lon, 0)); + ASSERT(approx(az, D90)); + + end_point(0, 0, D90, PE / 8, &lat, &lon, &az); + ASSERT(approx(lat, 0)); + ASSERT(approx(lon, D180 / 4)); + ASSERT(approx(az, D90)); + + end_point(0, 0, D90, PE / 4, &lat, &lon, &az); + ASSERT(approx(lat, 0)); + ASSERT(approx(lon, D90)); + ASSERT(approx(az, D90)); + + end_point(0, 0, D90, PE / 2, &lat, &lon, &az); + ASSERT(approx(lat, 0)); + ASSERT(approx(lon, D180)); + ASSERT(approx(az, D90)); + + return 0; +} + + #endif diff --git a/libtellurian_gaussian_radius.3 b/libtellurian_gaussian_radius.3 index aaa84dc..cd6033a 100644 --- a/libtellurian_gaussian_radius.3 +++ b/libtellurian_gaussian_radius.3 @@ -1,8 +1,8 @@ -.TH LIBTELLURIAN_GAUSSIAN_RADIUS 3 libtellurian +.TH LIBTELLURIAN_GAUSSIAN_RADIUS 3 LIBTELLURIAN .SH NAME libtellurian_gaussian_radius \- Calculate a Gaussian radius of curvature -.SH SYNPOSIS +.SH SYNOPSIS .nf #include <libtellurian.h> diff --git a/libtellurian_gaussian_radius_radians.c b/libtellurian_gaussian_radius_radians.c index 02d0586..7b14ed5 100644 --- a/libtellurian_gaussian_radius_radians.c +++ b/libtellurian_gaussian_radius_radians.c @@ -11,7 +11,7 @@ libtellurian_gaussian_radius_radians(double latitude) double c = b / a; double neg_e2 = fma(c, c, -1.0); double neg_e2_plus_1 = c * c; - double s = sin(latitude); + double s = sin(geodetic(latitude)); double denom = fma(neg_e2, s * s, 1.0); return a * sqrt(neg_e2_plus_1) / denom; } diff --git a/libtellurian_meridian_radius.3 b/libtellurian_meridian_radius.3 index a13b12b..433d4d4 100644 --- a/libtellurian_meridian_radius.3 +++ b/libtellurian_meridian_radius.3 @@ -1,8 +1,8 @@ -.TH LIBTELLURIAN_MERIDIAN_RADIUS 3 libtellurian +.TH LIBTELLURIAN_MERIDIAN_RADIUS 3 LIBTELLURIAN .SH NAME libtellurian_meridian_radius \- Calculate a meridian radius of curvature -.SH SYNPOSIS +.SH SYNOPSIS .nf #include <libtellurian.h> diff --git a/libtellurian_meridian_radius_radians.c b/libtellurian_meridian_radius_radians.c index f550289..d5213c1 100644 --- a/libtellurian_meridian_radius_radians.c +++ b/libtellurian_meridian_radius_radians.c @@ -9,7 +9,7 @@ libtellurian_meridian_radius_radians(double latitude) double a = LIBTELLURIAN_EQUATORIAL_RADIUS; double b = LIBTELLURIAN_POLAR_RADIUS; double neg_e2 = fma(b / a, b / a, -1.0); - double s = sin(latitude); + double s = sin(geodetic(latitude)); return fma(a, neg_e2, a) / pow(fma(neg_e2, s * s, 1.0), 1.5); } diff --git a/libtellurian_normal_gravity.3 b/libtellurian_normal_gravity.3 index cc3a6fa..44ca3b5 100644 --- a/libtellurian_normal_gravity.3 +++ b/libtellurian_normal_gravity.3 @@ -1,8 +1,8 @@ -.TH LIBTELLURIAN_NORMAL_GRAVITY 3 libtellurian +.TH LIBTELLURIAN_NORMAL_GRAVITY 3 LIBTELLURIAN .SH NAME libtellurian_normal_gravity \- Calculate normal gravity -.SH SYNPOSIS +.SH SYNOPSIS .nf #include <libtellurian.h> diff --git a/libtellurian_sea_level.3 b/libtellurian_sea_level.3 index 9401641..e6394b3 100644 --- a/libtellurian_sea_level.3 +++ b/libtellurian_sea_level.3 @@ -1,8 +1,8 @@ -.TH LIBTELLURIAN_SEA_LEVEL 3 libtellurian +.TH LIBTELLURIAN_SEA_LEVEL 3 LIBTELLURIAN .SH NAME libtellurian_sea_level \- Calculate a geocentric radius -.SH SYNPOSIS +.SH SYNOPSIS .nf #include <libtellurian.h> diff --git a/libtellurian_sea_level_radians.c b/libtellurian_sea_level_radians.c index e536025..ec03985 100644 --- a/libtellurian_sea_level_radians.c +++ b/libtellurian_sea_level_radians.c @@ -6,15 +6,9 @@ double libtellurian_sea_level_radians(double latitude) { - double a = LIBTELLURIAN_EQUATORIAL_RADIUS; - double b = LIBTELLURIAN_POLAR_RADIUS; - double c = cos(latitude); - double s = sin(latitude); - double x = a * c * a; - double y = b * s * b; - double num = fma(x, x, y * y); - double denom = fma(x, c, y * s); - return sqrt(num / denom); + double x = cos(latitude) * LIBTELLURIAN_EQUATORIAL_RADIUS; + double y = sin(latitude) * LIBTELLURIAN_POLAR_RADIUS; + return sqrt(fma(x, x, y * y)); } diff --git a/libtellurian_transverse_radius.3 b/libtellurian_transverse_radius.3 index aa24be4..27ceece 100644 --- a/libtellurian_transverse_radius.3 +++ b/libtellurian_transverse_radius.3 @@ -1,8 +1,8 @@ -.TH LIBTELLURIAN_TRANSVERSE_RADIUS 3 libtellurian +.TH LIBTELLURIAN_TRANSVERSE_RADIUS 3 LIBTELLURIAN .SH NAME libtellurian_transverse_radius \- Calculate a transverse radius of curvature -.SH SYNPOSIS +.SH SYNOPSIS .nf #include <libtellurian.h> diff --git a/libtellurian_transverse_radius_radians.c b/libtellurian_transverse_radius_radians.c index c4f80aa..3e4950e 100644 --- a/libtellurian_transverse_radius_radians.c +++ b/libtellurian_transverse_radius_radians.c @@ -9,7 +9,7 @@ libtellurian_transverse_radius_radians(double latitude) double a = LIBTELLURIAN_EQUATORIAL_RADIUS; double b = LIBTELLURIAN_POLAR_RADIUS; double neg_e2 = fma(b / a, b / a, -1.0); - double s = sin(latitude); + double s = sin(geodetic(latitude)); return a * pow(fma(neg_e2, s * s, 1.0), -0.5); } diff --git a/libtellurian_vincenty_inverse__.c b/libtellurian_vincenty_inverse__.c index c50752e..4d8d6ea 100644 --- a/libtellurian_vincenty_inverse__.c +++ b/libtellurian_vincenty_inverse__.c @@ -1,10 +1,9 @@ /* See LICENSE file for copyright and license details. */ #include "common.h" #ifndef TEST +#include <stdio.h> // TODO remove -/* TODO complete the implementation */ - void libtellurian_vincenty_inverse__(double latitude1, double longitude1, double latitude2, double longitude2, double *distance_out, double *azimuth1_out, double *azimuth2_out) @@ -34,42 +33,62 @@ libtellurian_vincenty_inverse__(double latitude1, double longitude1, double lati double sin_u1_sin_u2 = sin_u1 * sin_u2; double L = longitude2 - longitude1; + double lambda_prev; + int max_interations = 20; lambda = L; - /* { */ + do { + lambda_prev = lambda; - cos_lambda = cos(lambda); - sin_lambda = sin(lambda); + cos_lambda = cos(lambda); + sin_lambda = sin(lambda); - y = cos_u2 * sin_lambda; - x = fma(-sin_u1_cos_u2, cos_lambda, cos_u1_sin_u2); + y = cos_u2 * sin_lambda; + x = fma(-sin_u1_cos_u2, cos_lambda, cos_u1_sin_u2); - sin_sigma = sqrt(fma(y, y, x * x)); - cos_sigma = fma(cos_u1_cos_u2, cos_lambda, sin_u1_sin_u2); - sigma = atan2(sin_sigma, cos_sigma); + sin_sigma = sqrt(fma(y, y, x * x)); + cos_sigma = fma(cos_u1_cos_u2, cos_lambda, sin_u1_sin_u2); + sigma = atan2(sin_sigma, cos_sigma); - sin_alpha = (cos_u1_cos_u2 * sin_lambda) / sin_sigma; - cos2_alpha = fma(-sin_alpha, sin_alpha, 1.0); - /* - * If sin σ = 0 the value of sin α is indeterminate. - * It represents an end point coincident with, or - * diametrically opposed to, the start point. - */ + sin_alpha = (cos_u1_cos_u2 * sin_lambda) / sin_sigma; + if (!isfinite(sin_alpha)) { + fprintf(stderr, "BAD\n"); + if (distance_out) { + /* Dot product of vectors is 0 if coincidental, but π if antipodal */ + /* TODO it may be fast to see if latitude1+latitude2≋0 and longitude2-longitude1≋π */ + t = sin(latitude1) * sin(latitude2); + t = fma(cos(latitude1) * cos(latitude2), cos(longitude2 - longitude1), t); + if (t > 0.0) { + fprintf(stderr, " COINCIDENTAL\n"); + *distance_out = 0; + } else { + fprintf(stderr, " ANTIPODAL\n"); + *distance_out = LIBTELLURIAN_MERIDIONAL_CIRCUMFERENCE; + } + } + if (azimuth1_out) + *azimuth1_out = nan(""); + if (azimuth2_out) + *azimuth2_out = nan(""); + return; + } + cos2_alpha = fma(-sin_alpha, sin_alpha, 1.0); - cos_2sigma_m = fma(-2.0 / cos2_alpha, sin_u1_sin_u2, cos_sigma); - C = 0.25 * f * cos2_alpha * fma(f, fma(-0.75, cos2_alpha, 1.0), 1.0); + cos_2sigma_m = fma(-2.0 / cos2_alpha, sin_u1_sin_u2, cos_sigma); + C = 0.25 * f * cos2_alpha * fma(f, fma(-0.75, cos2_alpha, 1.0), 1.0); - cos2_2sigma_m = cos_2sigma_m * cos_2sigma_m; - t = fma(2.0, cos2_2sigma_m, -1.0); - t = fma(C * cos_sigma, t, cos_2sigma_m); - t = fma(C * sin_sigma, t, sigma); - lambda = fma(fma(C, f, -f) * sin_alpha, t, L); + cos2_2sigma_m = cos_2sigma_m * cos_2sigma_m; + t = fma(2.0, cos2_2sigma_m, -1.0); + t = fma(C * cos_sigma, t, cos_2sigma_m); + t = fma(C * sin_sigma, t, sigma); + lambda = fma(fma(C, f, -f) * sin_alpha, t, L); - /* } repeat until lambda converges */ + } while (lambda != lambda_prev && --max_interations); /* repeat until lambda converges */ + fprintf(stderr, "%lg => %lg -> %lg : %lg (%i)\n", L, lambda_prev, lambda, lambda - lambda_prev, max_interations); if (distance_out) { - uu = cos2_alpha * fma(a / b, a / b, -1.0); + uu = cos2_alpha * fma(c, c, -1.0); A = fma(fma(fma(fma(-175.0, uu, 320.0), uu, -768.0), uu, 4096.0), uu / 16384.0, 1.0); B = fma(fma(fma(-47.0, uu, 74.0), uu, -128.0), uu, 256.0) * (uu / 1024.0); @@ -89,7 +108,7 @@ libtellurian_vincenty_inverse__(double latitude1, double longitude1, double lati if (azimuth2_out) { y = cos_u1 * sin_lambda; x = fma(cos_u1_sin_u2, cos_lambda, -sin_u1_cos_u2); - *azimuth1_out = atan2(y, x); + *azimuth2_out = atan2(y, x); } /* @@ -106,5 +125,77 @@ libtellurian_vincenty_inverse__(double latitude1, double longitude1, double lati #else -TODO_TEST + + +# define PE LIBTELLURIAN_EQUATORIAL_CIRCUMFERENCE +# define PM LIBTELLURIAN_MERIDIONAL_CIRCUMFERENCE + +static int +approx(double a, double e) +{ + return fabs((a / e) - 1.0) <= 1.0e-8; +} + +int +main(void) +{ + double s, a1, a2; + +#define RESET ((void)(s = 111, a1 = 222, a2 = 333)) + + libtellurian_vincenty_inverse__(1, 2, 3, 4, &s, &a1, &a2); + libtellurian_vincenty_inverse__(1, 2, 3, 4, &s, &a1, NULL); + libtellurian_vincenty_inverse__(1, 2, 3, 4, &s, NULL, &a2); + libtellurian_vincenty_inverse__(1, 2, 3, 4, &s, NULL, NULL); + libtellurian_vincenty_inverse__(1, 2, 3, 4, NULL, &a1, &a2); + libtellurian_vincenty_inverse__(1, 2, 3, 4, NULL, &a1, NULL); + libtellurian_vincenty_inverse__(1, 2, 3, 4, NULL, NULL, &a2); + libtellurian_vincenty_inverse__(1, 2, 3, 4, NULL, NULL, NULL); + + RESET; + libtellurian_vincenty_inverse__(0, 0, 0, 0, &s, &a1, &a2); + ASSERT(s == 0); + ASSERT(isnan(a1)); + ASSERT(isnan(a2)); + + RESET; + libtellurian_vincenty_inverse__(0, 0, 0, 0, NULL, &a1, &a2); + ASSERT(isnan(a1)); + ASSERT(isnan(a2)); + + RESET; + libtellurian_vincenty_inverse__(0, 0, 0, 0, &s, NULL, &a2); + ASSERT(s == 0); + ASSERT(isnan(a2)); + + RESET; + libtellurian_vincenty_inverse__(0, 0, 0, 0, NULL, NULL, &a2); + ASSERT(isnan(a2)); + + RESET; + libtellurian_vincenty_inverse__(0, 0, 0, 0, &s, &a1, NULL); + ASSERT(s == 0); + ASSERT(isnan(a1)); + + RESET; + libtellurian_vincenty_inverse__(0, 0, 0, 0, NULL, &a1, NULL); + ASSERT(isnan(a1)); + + RESET; + libtellurian_vincenty_inverse__(0, 0, 0, 0, &s, NULL, NULL); + ASSERT(s == 0); + + libtellurian_vincenty_inverse__(0, 0, 0, 0, NULL, NULL, NULL); + + fprintf(stderr, "--------\n"); + RESET; + libtellurian_vincenty_inverse__(0, 0, D45, 0, &s, &a1, &a2); + ASSERT(approx(s, PM / 4)); + ASSERT(isnan(a1)); + ASSERT(isnan(a2)); + + return 0; +} + + #endif |