Skip to content

Apply MarceloBanik fix for #8717 #8720

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
joa-quim merged 1 commit into from
Apr 26, 2025
Merged

Apply MarceloBanik fix for #8717 #8720

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 17 additions & 14 deletions src/gmt_map.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9271,32 +9271,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);
}
}

Expand Down
Loading