Changeset 5801 in ntrip for trunk/BNC/src/PPP/pppModel.cpp

Timestamp:

Aug 6, 2014, 10:24:02 AM (10 years ago)

Author:

mervart

Message:

 

File:

• trunk/BNC/src/PPP/pppModel.cpp (moved) (moved from trunk/BNC/src/bnctides.cpp ) (8 diffs)

trunk/BNC/src/PPP/pppModel.cpp

@@ -2 +2 @@
 #include <cmath>
 
-#include "bnctides.h"
+#include "pppModel.h"
 #include "bncutils.h"
 
 using namespace std;
 
-// Auxiliary Functions
-///////////////////////////////////////////////////////////////////////////
-namespace { 
-
-  static const double RHO_DEG   = 180.0 / M_PI;
-  static const double RHO_SEC   = 3600.0 * RHO_DEG;
-  static const double MJD_J2000 = 51544.5;
-
-  double Frac (double x) { return x-floor(x); };
-  double Modulo (double x, double y) { return y*Frac(x/y); }
-
-  Matrix rotX(double Angle) {
-    const double C = cos(Angle);
-    const double S = sin(Angle);
-    Matrix UU(3,3);
-    UU[0][0] = 1.0;  UU[0][1] = 0.0;  UU[0][2] = 0.0;
-    UU[1][0] = 0.0;  UU[1][1] =  +C;  UU[1][2] =  +S;
-    UU[2][0] = 0.0;  UU[2][1] =  -S;  UU[2][2] =  +C;
-    return UU;
-  }
-  
-  Matrix rotY(double Angle) {
-    const double C = cos(Angle);
-    const double S = sin(Angle);
-    Matrix UU(3,3);
-    UU[0][0] =  +C;  UU[0][1] = 0.0;  UU[0][2] =  -S;
-    UU[1][0] = 0.0;  UU[1][1] = 1.0;  UU[1][2] = 0.0;
-    UU[2][0] =  +S;  UU[2][1] = 0.0;  UU[2][2] =  +C;
-    return UU;
-  }
-  
-  Matrix rotZ(double Angle) {
-    const double C = cos(Angle);
-    const double S = sin(Angle);
-    Matrix UU(3,3);
-    UU[0][0] =  +C;  UU[0][1] =  +S;  UU[0][2] = 0.0;
-    UU[1][0] =  -S;  UU[1][1] =  +C;  UU[1][2] = 0.0;
-    UU[2][0] = 0.0;  UU[2][1] = 0.0;  UU[2][2] = 1.0;
-    return UU;
-  }
+double Frac (double x) { return x-floor(x); };
+double Modulo (double x, double y) { return y*Frac(x/y); }
+
+Matrix t_astro::rotX(double Angle) {
+  const double C = cos(Angle);
+  const double S = sin(Angle);
+  Matrix UU(3,3);
+  UU[0][0] = 1.0;  UU[0][1] = 0.0;  UU[0][2] = 0.0;
+  UU[1][0] = 0.0;  UU[1][1] =  +C;  UU[1][2] =  +S;
+  UU[2][0] = 0.0;  UU[2][1] =  -S;  UU[2][2] =  +C;
+  return UU;
+}
+
+Matrix t_astro::rotY(double Angle) {
+  const double C = cos(Angle);
+  const double S = sin(Angle);
+  Matrix UU(3,3);
+  UU[0][0] =  +C;  UU[0][1] = 0.0;  UU[0][2] =  -S;
+  UU[1][0] = 0.0;  UU[1][1] = 1.0;  UU[1][2] = 0.0;
+  UU[2][0] =  +S;  UU[2][1] = 0.0;  UU[2][2] =  +C;
+  return UU;
+}
+
+Matrix t_astro::rotZ(double Angle) {
+  const double C = cos(Angle);
+  const double S = sin(Angle);
+  Matrix UU(3,3);
+  UU[0][0] =  +C;  UU[0][1] =  +S;  UU[0][2] = 0.0;
+  UU[1][0] =  -S;  UU[1][1] =  +C;  UU[1][2] = 0.0;
+  UU[2][0] = 0.0;  UU[2][1] = 0.0;  UU[2][2] = 1.0;
+  return UU;
 }
 
 // Greenwich Mean Sidereal Time
 ///////////////////////////////////////////////////////////////////////////
-double GMST(double Mjd_UT1) {
+double t_astro::GMST(double Mjd_UT1) {
 
   const double Secs = 86400.0;
@@ -68 +59 @@
 // Nutation Matrix
 ///////////////////////////////////////////////////////////////////////////
-Matrix NutMatrix(double Mjd_TT) {
+Matrix t_astro::NutMatrix(double Mjd_TT) {
 
   const double T  = (Mjd_TT-MJD_J2000)/36525.0;
@@ -89 +80 @@
 // Precession Matrix
 ///////////////////////////////////////////////////////////////////////////
-Matrix PrecMatrix (double Mjd_1, double Mjd_2) {
+Matrix t_astro::PrecMatrix(double Mjd_1, double Mjd_2) {
 
   const double T  = (Mjd_1-MJD_J2000)/36525.0;
@@ -105 +96 @@
 // Sun's position
 ///////////////////////////////////////////////////////////////////////////
-ColumnVector Sun(double Mjd_TT) {
+ColumnVector t_astro::Sun(double Mjd_TT) {
 
   const double eps = 23.43929111/RHO_DEG;
@@ -126 +117 @@
 // Moon's position
 ///////////////////////////////////////////////////////////////////////////
-ColumnVector Moon(double Mjd_TT) {
+ColumnVector t_astro::Moon(double Mjd_TT) {
 
   const double eps = 23.43929111/RHO_DEG;
@@ -169 +160 @@
 // Tidal Correction 
 ////////////////////////////////////////////////////////////////////////////
-void tides(const bncTime& time, ColumnVector& xyz) {
-
-  static double       lastMjd = 0.0;
-  static ColumnVector xSun;
-  static ColumnVector xMoon;
-  static double       rSun;
-  static double       rMoon;
+ColumnVector t_tides::displacement(const bncTime& time, const ColumnVector& xyz) {
 
   double Mjd = time.mjd() + time.daysec() / 86400.0;
 
-  if (Mjd != lastMjd) {
-    lastMjd = Mjd;
-    xSun = Sun(Mjd);
-    rSun = sqrt(DotProduct(xSun,xSun));
-    xSun /= rSun;
-    xMoon = Moon(Mjd);
-    rMoon = sqrt(DotProduct(xMoon,xMoon));
-    xMoon /= rMoon;
+  if (Mjd != _lastMjd) {
+    _lastMjd = Mjd;
+    _xSun = t_astro::Sun(Mjd);
+    _rSun = sqrt(DotProduct(_xSun,_xSun));
+    _xSun /= _rSun;
+    _xMoon = t_astro::Moon(Mjd);
+    _rMoon = sqrt(DotProduct(_xMoon,_xMoon));
+    _xMoon /= _rMoon;
   }
 
@@ -199 +184 @@
   // Tidal Displacement
   // ------------------
-  double scSun  = DotProduct(xyzUnit, xSun);
-  double scMoon = DotProduct(xyzUnit, xMoon);
+  double scSun  = DotProduct(xyzUnit, _xSun);
+  double scMoon = DotProduct(xyzUnit, _xMoon);
 
   double p2Sun  = 3.0 * (H2/2.0-L2) * scSun  * scSun  - H2/2.0;
@@ -213 +198 @@
 
   double facSun  = gms / gmWGS * 
-                   (rRec * rRec * rRec * rRec) / (rSun * rSun * rSun);
+                   (rRec * rRec * rRec * rRec) / (_rSun * _rSun * _rSun);
 
   double facMoon = gmm / gmWGS * 
-                   (rRec * rRec * rRec * rRec) / (rMoon * rMoon * rMoon);
-
-  ColumnVector dX = facSun  * (x2Sun  * xSun  + p2Sun  * xyzUnit) + 
-                    facMoon * (x2Moon * xMoon + p2Moon * xyzUnit);
-
-  xyz += dX;
-}
+                   (rRec * rRec * rRec * rRec) / (_rMoon * _rMoon * _rMoon);
+
+  ColumnVector dX = facSun  * (x2Sun  * _xSun  + p2Sun  * xyzUnit) + 
+                    facMoon * (x2Moon * _xMoon + p2Moon * xyzUnit);
+
+  return dX;
+}