tmp

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)

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)

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;
 }

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

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,