Changeset 5802 in ntrip for trunk/BNC/src/PPP/pppModel.cpp
- Timestamp:
- Aug 6, 2014, 10:35:19 AM (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/BNC/src/PPP/pppModel.cpp
r5801 r5802 208 208 return dX; 209 209 } 210 211 // Constructor 212 /////////////////////////////////////////////////////////////////////////// 213 t_windUp::t_windUp() { 214 for (unsigned ii = 0; ii <= t_prn::MAXPRN; ii++) { 215 sumWind[ii] = 0.0; 216 lastEtime[ii] = 0.0; 217 } 218 } 219 220 // Phase Wind-Up Correction 221 /////////////////////////////////////////////////////////////////////////// 222 double t_windUp::value(const bncTime& etime, const ColumnVector& rRec, 223 t_prn prn, const ColumnVector& rSat) { 224 225 if (etime.mjddec() != lastEtime[prn.toInt()]) { 226 227 // Unit Vector GPS Satellite --> Receiver 228 // -------------------------------------- 229 ColumnVector rho = rRec - rSat; 230 rho /= rho.norm_Frobenius(); 231 232 // GPS Satellite unit Vectors sz, sy, sx 233 // ------------------------------------- 234 ColumnVector sz = -rSat / rSat.norm_Frobenius(); 235 236 ColumnVector xSun = Sun(etime.mjddec()); 237 xSun /= xSun.norm_Frobenius(); 238 239 ColumnVector sy = crossproduct(sz, xSun); 240 ColumnVector sx = crossproduct(sy, sz); 241 242 // Effective Dipole of the GPS Satellite Antenna 243 // --------------------------------------------- 244 ColumnVector dipSat = sx - rho * DotProduct(rho,sx) 245 - crossproduct(rho, sy); 246 247 // Receiver unit Vectors rx, ry 248 // ---------------------------- 249 ColumnVector rx(3); 250 ColumnVector ry(3); 251 252 double recEll[3]; xyz2ell(rRec.data(), recEll) ; 253 double neu[3]; 254 255 neu[0] = 1.0; 256 neu[1] = 0.0; 257 neu[2] = 0.0; 258 neu2xyz(recEll, neu, rx.data()); 259 260 neu[0] = 0.0; 261 neu[1] = -1.0; 262 neu[2] = 0.0; 263 neu2xyz(recEll, neu, ry.data()); 264 265 // Effective Dipole of the Receiver Antenna 266 // ---------------------------------------- 267 ColumnVector dipRec = rx - rho * DotProduct(rho,rx) 268 + crossproduct(rho, ry); 269 270 // Resulting Effect 271 // ---------------- 272 double alpha = DotProduct(dipSat,dipRec) / 273 (dipSat.norm_Frobenius() * dipRec.norm_Frobenius()); 274 275 if (alpha > 1.0) alpha = 1.0; 276 if (alpha < -1.0) alpha = -1.0; 277 278 double dphi = acos(alpha) / 2.0 / M_PI; // in cycles 279 280 if ( DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0 ) { 281 dphi = -dphi; 282 } 283 284 if (lastEtime[prn.toInt()] == 0.0) { 285 sumWind[prn.toInt()] = dphi; 286 } 287 else { 288 sumWind[prn.toInt()] = nint(sumWind[prn.toInt()] - dphi) + dphi; 289 } 290 291 lastEtime[prn.toInt()] = etime.mjddec(); 292 } 293 294 return sumWind[prn.toInt()]; 295 }
Note:
See TracChangeset
for help on using the changeset viewer.