Changeset 7996 in ntrip for trunk/BNC/src/pppModel.cpp

Timestamp:

Aug 4, 2016, 9:51:01 AM (8 years ago)

Author:

stuerze

Message:

reading of ocean loading files is added

File:

• trunk/BNC/src/pppModel.cpp (modified) (16 diffs)

trunk/BNC/src/pppModel.cpp

@@ -36 +36 @@
  * Created:    29-Jul-2014
  *
- * Changes:    
+ * Changes:
  *
  * -----------------------------------------------------------------------*/
@@ -90 +90 @@
   double Mjd_0 = floor(Mjd_UT1);
   double UT1   = Secs*(Mjd_UT1-Mjd_0);
-  double T_0   = (Mjd_0  -MJD_J2000)/36525.0; 
-  double T     = (Mjd_UT1-MJD_J2000)/36525.0; 
+  double T_0   = (Mjd_0  -MJD_J2000)/36525.0;
+  double T     = (Mjd_UT1-MJD_J2000)/36525.0;
 
   double gmst  = 24110.54841 + 8640184.812866*T_0 + 1.002737909350795*UT1
@@ -126 +126 @@
   const double T  = (Mjd_1-MJD_J2000)/36525.0;
   const double dT = (Mjd_2-Mjd_1)/36525.0;
-  
+
   double zeta  =  ( (2306.2181+(1.39656-0.000139*T)*T)+
                         ((0.30188-0.000344*T)+0.017998*dT)*dT )*dT/RHO_SEC;
@@ -134 +134 @@
 
   return rotZ(-z) * rotY(theta) * rotZ(-zeta);
-}    
+}
 
 // Sun's position
@@ -144 +144 @@
 
   double M = 2.0*M_PI * Frac ( 0.9931267 + 99.9973583*T);
-  double L = 2.0*M_PI * Frac ( 0.7859444 + M/2.0/M_PI + 
+  double L = 2.0*M_PI * Frac ( 0.7859444 + M/2.0/M_PI +
                         (6892.0*sin(M)+72.0*sin(2.0*M)) / 1296.0e3);
   double r = 149.619e9 - 2.499e9*cos(M) - 0.021e9*cos(2*M);
-  
-  ColumnVector r_Sun(3); 
+
+  ColumnVector r_Sun(3);
   r_Sun << r*cos(L) << r*sin(L) << 0.0; r_Sun = rotX(-eps) * r_Sun;
 
   return    rotZ(GMST(Mjd_TT))
-          * NutMatrix(Mjd_TT) 
+          * NutMatrix(Mjd_TT)
           * PrecMatrix(MJD_J2000, Mjd_TT)
           * r_Sun;
@@ -169 +169 @@
   double D   = 2.0*M_PI*Frac ( 0.827361 + 1236.853086*T );
   double F   = 2.0*M_PI*Frac ( 0.259086 + 1342.227825*T );
-    
-  double dL = +22640*sin(l) - 4586*sin(l-2*D) + 2370*sin(2*D) +  769*sin(2*l) 
+
+  double dL = +22640*sin(l) - 4586*sin(l-2*D) + 2370*sin(2*D) +  769*sin(2*l)
               -668*sin(lp) - 412*sin(2*F) - 212*sin(2*l-2*D)- 206*sin(l+lp-2*D)
               +192*sin(l+2*D) - 165*sin(lp-2*D) - 125*sin(D) - 110*sin(l+lp)
@@ -177 +177 @@
   double L = 2.0*M_PI * Frac( L_0 + dL/1296.0e3 );
 
-  double S  = F + (dL+412*sin(2*F)+541*sin(lp)) / RHO_SEC; 
+  double S  = F + (dL+412*sin(2*F)+541*sin(lp)) / RHO_SEC;
   double h  = F-2*D;
-  double N  = -526*sin(h) + 44*sin(l+h) - 31*sin(-l+h) - 23*sin(lp+h) 
+  double N  = -526*sin(h) + 44*sin(l+h) - 31*sin(-l+h) - 23*sin(lp+h)
               +11*sin(-lp+h) - 25*sin(-2*l+F) + 21*sin(-l+F);
 
   double B = ( 18520.0*sin(S) + N ) / RHO_SEC;
-    
+
   double cosB = cos(B);
 
   double R = 385000e3 - 20905e3*cos(l) - 3699e3*cos(2*D-l) - 2956e3*cos(2*D)
-      -570e3*cos(2*l) + 246e3*cos(2*l-2*D) - 205e3*cos(lp-2*D) 
-      -171e3*cos(l+2*D) - 152e3*cos(l+lp-2*D);   
-
-  ColumnVector r_Moon(3); 
+      -570e3*cos(2*l) + 246e3*cos(2*l-2*D) - 205e3*cos(lp-2*D)
+      -171e3*cos(l+2*D) - 152e3*cos(l+lp-2*D);
+
+  ColumnVector r_Moon(3);
   r_Moon << R*cos(L)*cosB << R*sin(L)*cosB << R*sin(B);
   r_Moon = rotX(-eps) * r_Moon;
-    
-  return    rotZ(GMST(Mjd_TT)) 
-          * NutMatrix(Mjd_TT) 
+
+  return    rotZ(GMST(Mjd_TT))
+          * NutMatrix(Mjd_TT)
           * PrecMatrix(MJD_J2000, Mjd_TT)
           * r_Moon;
 }
 
-// Tidal Correction 
+// Tidal Correction
 ////////////////////////////////////////////////////////////////////////////
 ColumnVector t_tides::displacement(const bncTime& time, const ColumnVector& xyz) {
@@ -239 +239 @@
   double x2Sun  = 3.0 * L2 * scSun;
   double x2Moon = 3.0 * L2 * scMoon;
-  
+
   const double gmWGS = 398.6005e12;
   const double gms   = 1.3271250e20;
   const double gmm   = 4.9027890e12;
 
-  double facSun  = gms / gmWGS * 
+  double facSun  = gms / gmWGS *
                    (rRec * rRec * rRec * rRec) / (_rSun * _rSun * _rSun);
 
-  double facMoon = gmm / gmWGS * 
+  double facMoon = gmm / gmWGS *
                    (rRec * rRec * rRec * rRec) / (_rMoon * _rMoon * _rMoon);
 
-  ColumnVector dX = facSun  * (x2Sun  * _xSun  + p2Sun  * xyzUnit) + 
+  ColumnVector dX = facSun  * (x2Sun  * _xSun  + p2Sun  * xyzUnit) +
                     facMoon * (x2Moon * _xMoon + p2Moon * xyzUnit);
 
   return dX;
+}
+
+// Constructor
+///////////////////////////////////////////////////////////////////////////
+t_loading::t_loading(const QString& fileName) {
+
+  newBlqData = 0;
+  readFile(fileName);
+  printAll();
+}
+
+t_loading::~t_loading() {
+
+  if (newBlqData) {
+    delete newBlqData;
+  }
+
+  QMapIterator<QString, t_blqData*> it(blqMap);
+  while (it.hasNext()) {
+    it.next();
+    delete it.value();
+  }
+}
+
+t_irc t_loading::readFile(const QString& fileName) {
+  QFile inFile(fileName);
+  inFile.open(QIODevice::ReadOnly | QIODevice::Text);
+  QTextStream in(&inFile);
+  int row = 0;
+  QString site = QString();
+
+  while ( !in.atEnd() ) {
+
+    QString line = in.readLine();
+
+    // skip empty lines and comments
+    if (line.indexOf("$$") != -1 || line.size() == 0) {
+      continue;
+    }
+
+    line = line.trimmed();
+    QTextStream inLine(line.toAscii(), QIODevice::ReadOnly);
+
+    switch (row) {
+      case 0:
+        site = line;
+        site = site.toUpper();
+        if (!newBlqData) {
+          newBlqData = new t_blqData;
+          newBlqData->amplitudes.ReSize(3,11);
+          newBlqData->phases.ReSize(3,11);
+        }
+        break;
+      case 1: case 2: case 3:
+        for (int ii = 0; ii < 11; ii++) {
+          inLine >> newBlqData->amplitudes[row-1][ii];
+        }
+        break;
+      case 4: case 5: case 6:
+        for (int ii = 0; ii < 11; ii++) {
+          inLine >> newBlqData->phases[row-4][ii];
+        }
+        break;
+      case 7:
+        if (newBlqData && !site.isEmpty()) {
+          blqMap[site] = newBlqData;
+          site = QString();
+          row = -1;
+          newBlqData = 0;
+        }
+        break;
+    }
+    row++;
+  }
+  inFile.close();
+  return success;
+}
+
+// Print
+////////////////////////////////////////////////////////////////////////////
+void t_loading::printAll() const {
+
+  QMapIterator<QString, t_blqData*> it(blqMap);
+  while (it.hasNext()) {
+    it.next();
+    t_blqData* blq = it.value();
+    QString site = it.key();
+    cout << site.toStdString().c_str() << "\n";
+    for (int rr = 0; rr < 3; rr++) {
+      for (int cc = 0; cc < 11; cc++) {
+        cout << blq->amplitudes[rr][cc] << " ";
+      }
+      cout << endl;
+    }
+    for (int rr = 0; rr < 3; rr++) {
+      for (int cc = 0; cc < 11; cc++) {
+        cout << blq->phases[rr][cc] << " ";
+      }
+      cout << endl;
+    }
+  }
 }
 
@@ -267 +368 @@
 // Phase Wind-Up Correction
 ///////////////////////////////////////////////////////////////////////////
-double t_windUp::value(const bncTime& etime, const ColumnVector& rRec, 
+double t_windUp::value(const bncTime& etime, const ColumnVector& rRec,
                        t_prn prn, const ColumnVector& rSat) {
 
@@ -276 +377 @@
     ColumnVector rho = rRec - rSat;
     rho /= rho.norm_Frobenius();
-    
+
     // GPS Satellite unit Vectors sz, sy, sx
     // -------------------------------------
@@ -289 +390 @@
     // Effective Dipole of the GPS Satellite Antenna
     // ---------------------------------------------
-    ColumnVector dipSat = sx - rho * DotProduct(rho,sx) 
+    ColumnVector dipSat = sx - rho * DotProduct(rho,sx)
                                                 - crossproduct(rho, sy);
-    
+
     // Receiver unit Vectors rx, ry
     // ----------------------------
@@ -299 +400 @@
     double recEll[3]; xyz2ell(rRec.data(), recEll) ;
     double neu[3];
-    
+
     neu[0] = 1.0;
     neu[1] = 0.0;
     neu[2] = 0.0;
     neu2xyz(recEll, neu, rx.data());
-    
+
     neu[0] =  0.0;
     neu[1] = -1.0;
     neu[2] =  0.0;
     neu2xyz(recEll, neu, ry.data());
-    
+
     // Effective Dipole of the Receiver Antenna
     // ----------------------------------------
-    ColumnVector dipRec = rx - rho * DotProduct(rho,rx) 
+    ColumnVector dipRec = rx - rho * DotProduct(rho,rx)
                                                    + crossproduct(rho, ry);
-    
+
     // Resulting Effect
     // ----------------
-    double alpha = DotProduct(dipSat,dipRec) / 
+    double alpha = DotProduct(dipSat,dipRec) /
                       (dipSat.norm_Frobenius() * dipRec.norm_Frobenius());
-    
+
     if (alpha >  1.0) alpha =  1.0;
     if (alpha < -1.0) alpha = -1.0;
-    
+
     double dphi = acos(alpha) / 2.0 / M_PI;  // in cycles
-    
+
     if ( DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0 ) {
       dphi = -dphi;
@@ -339 +440 @@
   }
 
-  return sumWind[prn.toInt()];  
+  return sumWind[prn.toInt()];
 }
 
@@ -352 +453 @@
   }
 
-  double ell[3]; 
+  double ell[3];
   xyz2ell(xyz.data(), ell);
   double height = ell[2];
@@ -362 +463 @@
 
   double h_km = height / 1000.0;
-  
+
   if (h_km < 0.0) h_km = 0.0;
   if (h_km > 5.0) h_km = 5.0;
   int    ii   = int(h_km + 1); if (ii > 5) ii = 5;
   double href = ii - 1;
-  
-  double bCor[6]; 
+
+  double bCor[6];
   bCor[0] = 1.156;
   bCor[1] = 1.006;
@@ -375 +476 @@
   bCor[4] = 0.654;
   bCor[5] = 0.563;
-  
+
   double BB = bCor[ii-1] + (bCor[ii]-bCor[ii-1]) * (h_km - href);
-  
+
   double zen  = M_PI/2.0 - Ele;