diff options
Diffstat (limited to '')
| -rw-r--r-- | Makefile | 4 | ||||
| -rw-r--r-- | TODO | 4 | ||||
| -rw-r--r-- | libtellurian.h | 68 | ||||
| -rw-r--r-- | libtellurian_coarse_distance.3 | 13 | ||||
| -rw-r--r-- | libtellurian_coarse_distance_radians.c | 4 | ||||
| l--------- | libtellurian_coarse_distance_radius.3 | 1 | ||||
| -rw-r--r-- | libtellurian_distance.3 | 15 | ||||
| -rw-r--r-- | libtellurian_end_point.3 | 25 | ||||
| -rw-r--r-- | libtellurian_end_point.c | 73 | ||||
| -rw-r--r-- | libtellurian_end_point_radians.c | 162 |
10 files changed, 318 insertions, 51 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,5 +1 @@ -Finish libtellurian_end_point_radians Finish libtellurian_vincenty_inverse__ - -Add a method for creating a locally optimised spherical model of the earth, -and add method for calculating distances on this model diff --git a/libtellurian.h b/libtellurian.h index 609f61d..b294e9c 100644 --- a/libtellurian.h +++ b/libtellurian.h @@ -132,6 +132,18 @@ /** + * 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 @@ -191,14 +203,14 @@ double libtellurian_coarse_distance_radians(double latitude1, double longitude1, * This function is 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 @@ -219,14 +231,14 @@ double libtellurian_distance(double latitude1, double longitude1, * This function is 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 @@ -247,14 +259,14 @@ double libtellurian_distance_radians(double latitude1, double longitude1, * This function is 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` * except it does not compute the distance between the @@ -270,14 +282,14 @@ void libtellurian_azimuth(double latitude1, double longitude1, * This function is 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` @@ -299,8 +311,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,8 +329,8 @@ 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); diff --git a/libtellurian_coarse_distance.3 b/libtellurian_coarse_distance.3 index 9a3973b..a490ba7 100644 --- a/libtellurian_coarse_distance.3 +++ b/libtellurian_coarse_distance.3 @@ -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..fc753af 100644 --- a/libtellurian_distance.3 +++ b/libtellurian_distance.3 @@ -104,6 +104,21 @@ functions do not return any value. .SH ERRORS None. +.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_end_point.3 b/libtellurian_end_point.3 index 15fe232..3bf78b2 100644 --- a/libtellurian_end_point.3 +++ b/libtellurian_end_point.3 @@ -8,12 +8,12 @@ libtellurian_end_point \- Calculate travel end-point 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 @@ -80,5 +83,21 @@ None. .SH ERRORS None. +.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..036ac99 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 = 0; + 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 |