Apply MarceloBanik fix for #8717 (#8720)

Author: joa-quim

src/gmt_map.c

@@ -9274,32 +9274,35 @@ void gmt_ECEF_forward (struct GMT_CTRL *GMT, double in[], double out[]) {
 }
 
 /*! . */
-void gmt_ECEF_inverse (struct GMT_CTRL *GMT, double in[], double out[]) {
+void gmt_ECEF_inverse(struct GMT_CTRL *GMT, double in[], double out[]) {
 	/* Convert ECEF coordinates to geodetic lon, lat, height given the datum parameters.
 	 * GMT->current.proj.datum.from is always the ellipsoid to use */
 
 	unsigned int i;
-	double in_p[3], sin_lat, cos_lat, N, p, theta, sin_theta, cos_theta;
+	double in_p[3], sin_lat, cos_lat, N, p, r, theta, sin_theta, cos_theta, slam, clam, sphi, cphi;
 
 	/* First remove the xyz shifts, us in_p to avoid changing in */
 	for (i = 0; i < 3; i++) in_p[i] = in[i] - GMT->current.proj.datum.from.xyz[i];
 
 	p = hypot (in_p[GMT_X], in_p[GMT_Y]);
-	if (p > 0.0) {	/* Not the S|N pole so we can invert */
-		theta = atan (in_p[GMT_Z] * GMT->current.proj.datum.from.a / (p * GMT->current.proj.datum.from.b));
-		sincos (theta, &sin_theta, &cos_theta);
-		out[GMT_X] = d_atan2d (in_p[GMT_Y], in_p[GMT_X]);
-		out[GMT_Y] = atand ((in_p[GMT_Z] + GMT->current.proj.datum.from.ep_squared * GMT->current.proj.datum.from.b * pow (sin_theta, 3.0)) / (p - GMT->current.proj.datum.from.e_squared * GMT->current.proj.datum.from.a * pow (cos_theta, 3.0)));
-		if (in_p[GMT_Z] > 0.0 && out[GMT_Y] < 0.0) out[GMT_Y] = -out[GMT_Y];
-		if (in_p[GMT_Z] < 0.0 && out[GMT_Y] > 0.0) out[GMT_Y] = -out[GMT_Y];
+	slam = p != 0 ? in_p[GMT_Y] / p : 0;
+	clam = p != 0 ? in_p[GMT_X] / p : 1;
+	theta = atan(in_p[GMT_Z] * GMT->current.proj.datum.from.a / (p * GMT->current.proj.datum.from.b));
+	sincos(theta, &sin_theta, &cos_theta);
+	sphi = in_p[GMT_Z] + GMT->current.proj.datum.from.ep_squared * GMT->current.proj.datum.from.b * pow(sin_theta, 3.0);
+	cphi = p - GMT->current.proj.datum.from.e_squared * GMT->current.proj.datum.from.a * pow(cos_theta, 3.0);
+	out[GMT_X] = d_atan2d(slam, clam);
+	out[GMT_Y] = d_atan2d(sphi, cphi);
+	if (in_p[GMT_Z] > 0.0 && out[GMT_Y] < 0.0) out[GMT_Y] = -out[GMT_Y];
+	if (in_p[GMT_Z] < 0.0 && out[GMT_Y] > 0.0) out[GMT_Y] = -out[GMT_Y];
+	if (fabs(out[GMT_Y]) < 89.999) {
 		sincosd (out[GMT_Y], &sin_lat, &cos_lat);
-		N = GMT->current.proj.datum.from.a / sqrt (1.0 - GMT->current.proj.datum.from.e_squared * sin_lat * sin_lat);
+		N = GMT->current.proj.datum.from.a / sqrt(1.0 - GMT->current.proj.datum.from.e_squared * sin_lat * sin_lat);
 		out[GMT_Z] = (p / cos_lat) - N;
 	}
-	else {	/* S or N pole, use sign of in_p[GMT_Z] to set latitude and height */
-		out[GMT_X] = 0.0;	/* Might as well pick0 since any longitude will work */
-		out[GMT_Y] = (in_p[GMT_Z] > 0.0) ? 90.0 : -90.0;	/* EIther at north or south pole, check via Z coordinate */
-		out[GMT_Z] = in_p[GMT_Z] - copysign (GMT->current.proj.datum.from.b, in_p[GMT_Z]);
+	else {		/* very close to the poles */
+		r = hypot(p, in_p[GMT_Z]);
+		out[GMT_Z] = r - hypot(GMT->current.proj.datum.from.b, p);
 	}
 }