...

Signed-off-by: Mattias Andrée <m@maandree.se>

10 files changed, 318 insertions, 51 deletions

diff --git a/Makefile b/Makefile
index 99e7f4d..ef10117 100644
--- a/Makefile
+++ b/Makefile
@@ -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\
diff --git a/TODO b/TODO
index 0b30e62..f022c99 100644
--- a/TODO
+++ b/TODO
@@ -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