From 45d38f5f50b95256874c633a0de8364da08a47f7 Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Thu, 6 Mar 2025 11:00:59 +0800 Subject: [PATCH] tmp --- example/CMakeLists.txt | 2 +- example/refellipsoid_ex.cpp | 21 ++++++++++++--------- lib/geometry/refellipsoid.cpp | 20 ++++++++------------ lib/geometry/refellipsoid.h | 8 ++++---- lib/maths/mathfunc.cpp | 4 ++-- 5 files changed, 27 insertions(+), 28 deletions(-) diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index a02bcba..7b50a73 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -21,7 +21,7 @@ add_example(fir_filter_ex OFF) add_example(fft_filter_ex OFF) add_example(windowfunc_ex OFF) add_example(legendre_ex OFF) -add_example(refellipsoid_ex OFF) +add_example(refellipsoid_ex ON) add_example(kde_ex OFF) add_example(meshio_ex OFF) add_example(autodiff_ex OFF) diff --git a/example/refellipsoid_ex.cpp b/example/refellipsoid_ex.cpp index 650dd64..0f297e4 100644 --- a/example/refellipsoid_ex.cpp +++ b/example/refellipsoid_ex.cpp @@ -53,16 +53,19 @@ int main(int argc, char const *argv[]) try point3ds ps; ps.set_io_precision(17); - // 注意这是一个地心纬度为40的椭球上的点 + // 大地坐标 (110,40,0) ps.lon = 110; - ps.lat = 40; - ps.rad = ellip.radius_at(ps.lat); - - point3dc pc; - pc.set_io_precision(17); - pc = ps.s2c(); - ellip.xyz2geodetic(pc, ps.lon, ps.lat, ps.rad); - std::cout << "Geodetic: " << ps << std::endl; + ellip.geodetic2spherical(40.0, 0.0, ps.lat, ps.rad); + // 球心坐标 (110,39.810615323061491,6369344.5441424493) + std::cout << "Geocentric: " << ps << std::endl; + // 转换为大地坐标 + point3ds pd; + ellip.xyz2geodetic(ps.s2c(), pd.lon, pd.lat, pd.rad); + std::cout << "Geodetic: " << pd << std::endl; + // 400km高 + pd.rad = 400000.0; + ellip.geodetic2spherical(pd.lat, pd.rad, pd.lat, pd.rad); + std::cout << "Geocentric: " << pd << std::endl; return 0; } catch(std::exception &e) diff --git a/lib/geometry/refellipsoid.cpp b/lib/geometry/refellipsoid.cpp index 97e6beb..5ca1d37 100644 --- a/lib/geometry/refellipsoid.cpp +++ b/lib/geometry/refellipsoid.cpp @@ -52,7 +52,7 @@ void gctl::refellipsoid::set(refellipsoid_type_e refellipsoid) else throw std::invalid_argument("Invalid reference system type for gctl::refellipsoid::set(...)"); f_ = (R_ - r_)/R_; - e_ = sqrt(R_*R_ - r_*r_)/R_; + e_ = sqrt(1.0 - (r_*r_)/(R_*R_)); return; } @@ -63,14 +63,13 @@ void gctl::refellipsoid::set(double R, double r_or_flat, bool is_flat) else r_ = r_or_flat; f_ = (R_ - r_)/R_; - e_ = sqrt(R_*R_ - r_*r_)/R_; + e_ = sqrt(1.0 - (r_*r_)/(R_*R_)); return; } -double gctl::refellipsoid::radius_at(double geocentric_lati) +double gctl::refellipsoid::geodetic_radius(double geodetic_lati) { - return ellipse_radius_2d(R_, r_, geocentric_lati*GCTL_Pi/180.0); - //return R_/sqrt(1.0 + e_*e_*sind(geocentric_lati)*sind(geocentric_lati)); + return R_/sqrt(1.0 - e_*e_*sind(geodetic_lati)*sind(geodetic_lati)); } void gctl::refellipsoid::geodetic2spherical(double geodetic_lati, @@ -85,18 +84,15 @@ void gctl::refellipsoid::geodetic2spherical(double geodetic_lati, CosLat = cosd(geodetic_lati); SinLat = sind(geodetic_lati); - double eps = sqrt(1.0 - (r_*r_)/(R_*R_)); - double epssq = eps*eps; - // compute the local rRdius of curvature on the WGS-84 reference ellipsoid - rc = R_ / sqrt(1.0 - epssq * SinLat * SinLat); + rc = R_/sqrt(1.0 - e_*e_*SinLat*SinLat); // compute ECEF Cartesian coordinates of specified point (for longitude=0) - xp = (rc + geodetic_hei) * CosLat; - zp = (rc*(1.0 - epssq) + geodetic_hei) * SinLat; + xp = (rc + geodetic_hei)*CosLat; + zp = (rc*(1.0 - e_*e_) + geodetic_hei)*SinLat; // compute spherical rRdius and angle lambda and phi of specified point - sph_rad = sqrt(xp * xp + zp * zp); + sph_rad = sqrt(xp*xp + zp*zp); sph_lati = deg(asin(zp/sph_rad)); return; } diff --git a/lib/geometry/refellipsoid.h b/lib/geometry/refellipsoid.h index 07f37b4..8e7ab0e 100644 --- a/lib/geometry/refellipsoid.h +++ b/lib/geometry/refellipsoid.h @@ -78,12 +78,12 @@ namespace gctl void set(double R, double r_or_flat, bool is_flat = false); /** - * @brief Get the ellipsoidal radius at the inquiring latitude + * @brief Compute the local rRdius of curvature on the reference ellipsoid * - * @param geocentric_lati the inquiring latitude - * @return ellipsoidal radius + * @param geodetic_lati the inquiring latitude + * @return radius */ - double radius_at(double geocentric_lati); + double geodetic_radius(double geodetic_lati); /** * @brief Converts Geodetic coordinates to Spherical coordinates diff --git a/lib/maths/mathfunc.cpp b/lib/maths/mathfunc.cpp index 1bf4399..6eb84a5 100644 --- a/lib/maths/mathfunc.cpp +++ b/lib/maths/mathfunc.cpp @@ -128,12 +128,12 @@ double gctl::geographic_distance(double lon1, double lon2, double lat1, double l double gctl::ellipse_radius_2d(double x_len, double y_len, double arc, double x_arc) { - if (fabs(x_len - y_len) < 1e-8) // 就是个圆 直接加 + if (fabs(x_len - y_len) < 1e-16) // 就是个圆 直接加 { return 0.5*(x_len + y_len); } - return sqrt(pow(x_len*cos(arc-x_arc),2) + pow(y_len*sin(arc-x_arc),2)); + return sqrt(power2(x_len*cos(arc - x_arc)) + power2(y_len*sin(arc - x_arc))); } void gctl::ellipse_plus_elevation_2d(double x_len, double y_len, double arc, double elev,