tmp

lib/core/macro.h

@@ -43,12 +43,13 @@
 /**
  * physical constants
  */
-#define GCTL_WGS84_PoleRadius 6356752.314 ///< WGS84椭球极半径
+#define GCTL_WGS84_PoleRadius 6356752.3142 ///< WGS84椭球极半径
-#define GCTL_WGS84_EquatorRadius 6378137 ///< WGS84椭球长半径
+#define GCTL_WGS84_EquatorRadius 6378136.460 ///< WGS84椭球长半径
-#define GCTL_WGS84_FLAT 3.35281066474748e-3 ///<  (1/298.25722356300)
+#define GCTL_WGS84_FLAT 3.35281066474748e-3 ///<  (1.0/298.257223563)
 #define GCTL_WGS84_NormPotential 6.263685171456948e+07 ///< m**2/s**2
 #define GCTL_Earth_GM 3.986004418e+14 ///< m**3/s**2
 #define GCTL_Earth_Radius 6371008.8 ///< 地球平均半径
+#define GCTL_Earth_RefRadius 6371200.0 ///< 地球参考半径（常用于地磁场建模中）
 #define GCTL_Moon_Radius 1738000 ///< 月球平均半径
 // Ardalan A. A., Karimi R. & Grafarend E. W. (2010). A new reference equipotential surface, and reference ellipsoid for the planet mars.
 // Earth Moon Planet, 106:1-13

lib/geometry/refellipsoid.cpp

@@ -44,21 +44,55 @@ gctl::refellipsoid::~refellipsoid(){}
 void gctl::refellipsoid::set(refellipsoid_type_e refellipsoid)
 {
     if (refellipsoid == Earth) {r_ = R_ = GCTL_Earth_Radius;}
+    else if (refellipsoid == MagEarth) {r_ = R_ = GCTL_Earth_RefRadius;}
     else if (refellipsoid == Moon) {r_ = R_ = GCTL_Moon_Radius;}
     else if (refellipsoid == Mars) {r_ = R_ = GCTL_Mars_Radius;}
     else if (refellipsoid == WGS84) {r_ = GCTL_WGS84_PoleRadius; R_ = GCTL_WGS84_EquatorRadius;}
     else if (refellipsoid == Ardalan2010) {r_ = GCTL_Mars_Ardalan2010_PoleRadius; R_ = GCTL_Mars_Ardalan2010_EquatorRadius;}
     else throw std::invalid_argument("Invalid reference system type for gctl::refellipsoid::set(...)");
+
+    eps_ = sqrt(1.0 - (r_*r_)/(R_*R_));
+    epssq_ = eps_*eps_;
+    return;
 }
 
 void gctl::refellipsoid::set(double R, double r_or_flat, bool is_flat)
 {
     R_ = R;
-    if (is_flat) r_ = R_*(1 - 1/r_or_flat);
+    if (is_flat) r_ = R_*(1.0 - 1.0/r_or_flat);
     else r_ = r_or_flat;
+
+    eps_ = sqrt(1.0 - (r_*r_)/(R_*R_));
+    epssq_ = eps_*eps_;
+    return;
 }
 
 double gctl::refellipsoid::radius_at(double lati)
 {
     return ellipse_radius_2d(R_, r_, lati*M_PI/180.0, 0.0);
 }
+
+void gctl::refellipsoid::geodetic2spherical(double geodetic_lati, 
+    double geodetic_hei, double& sph_lati, double& sph_rad)
+{
+    double CosLat, SinLat, rc, xp, zp;
+    /*
+    ** Convert geodetic coordinates, (defined by the WGS-84
+    ** reference ellipsoid), to Earth Centered Earth Fixed Cartesian
+    ** coordinates, and then to spherical coordinates.
+    */
+   CosLat = cosd(geodetic_lati);
+   SinLat = sind(geodetic_lati);
+
+   // compute the local rRdius of curvature on the WGS-84 reference ellipsoid
+   rc = R_ / sqrt(1.0 - epssq_ * 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;
+
+   // compute spherical rRdius and angle lambda and phi of specified point
+   sph_rad = sqrt(xp * xp + zp * zp);
+   sph_lati = deg(asin(zp/sph_rad));
+   return;
+}

lib/geometry/refellipsoid.h

@@ -38,6 +38,7 @@ namespace gctl
     {
         // reference system named as the planet name represents a mean radius reference sphere
         Earth,
+        MagEarth,
         Moon,
         Mars,
         // Reference ellipsoids 
@@ -82,8 +83,19 @@ namespace gctl
          */
         double radius_at(double lati);
 
+        /**
+         * @brief Converts Geodetic coordinates to Spherical coordinates
+         * 
+         * @param geodetic_lati Geodetic latitude
+         * @param geodetic_hei Height above the ellipsoid
+         * @param sph_lati Spherical latitude
+         * @param sph_rad Spherical radius
+         */
+        void geodetic2spherical(double geodetic_lati, double geodetic_hei, 
+            double& sph_lati, double& sph_rad);
+
     private:
-        double r_, R_;
+        double r_, R_, eps_, epssq_;
     };
     
 };

lib/maths/mathfunc_t.h

@@ -42,31 +42,31 @@ namespace gctl
     template <typename T>
     inline T arc(T deg)
     {
-        return deg*GCTL_Pi/180.0;
+        return deg*GCTL_Pi/180.0L;
     }
 
     template <typename T>
     inline T deg(T arc)
     {
-        return arc*180.0/GCTL_Pi;
+        return arc*180.0L/GCTL_Pi;
     }
 
     template <typename T>
     inline T sind(T deg)
     {
-        return sin(deg*GCTL_Pi/180.0);
+        return sin(deg*GCTL_Pi/180.0L);
     }
 
     template <typename T>
     inline T cosd(T deg)
     {
-        return cos(deg*GCTL_Pi/180.0);
+        return cos(deg*GCTL_Pi/180.0L);
     }
 
     template <typename T>
     inline T tand(T deg)
     {
-        return tan(deg*GCTL_Pi/180.0);
+        return tan(deg*GCTL_Pi/180.0L);
     }
 
 	template <typename T>