Changeset 2582 in ntrip for trunk/BNC/bncmodel.cpp
- Timestamp:
- Aug 29, 2010, 3:45:01 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/BNC/bncmodel.cpp
r2580 r2582 960 960 QQ << (SS.t() * SS); 961 961 } 962 963 // Phase Wind-Up Correction 964 /////////////////////////////////////////////////////////////////////////// 965 double bncModel::windUp(const QString& prn, const ColumnVector& rSat, 966 const ColumnVector& rRec) { 967 968 double Mjd = _time.mjd() + _time.daysec() / 86400.0; 969 970 // First time - initialize to zero 971 // ------------------------------- 972 if (!_windUpTime.contains(prn)) { 973 _windUpTime[prn] = Mjd; 974 _windUpSum[prn] = 0.0; 975 } 976 977 // Compute the correction for new time 978 // ----------------------------------- 979 else if (_windUpTime[prn] != Mjd) { 980 _windUpTime[prn] = Mjd; 981 982 // Unit Vector GPS Satellite --> Receiver 983 // -------------------------------------- 984 ColumnVector rho = rRec - rSat; 985 rho /= rho.norm_Frobenius(); 986 987 // GPS Satellite unit Vectors sz, sy, sx 988 // ------------------------------------- 989 ColumnVector sz = -rSat / rSat.norm_Frobenius(); 990 991 ColumnVector xSun = Sun(Mjd); 992 xSun /= xSun.norm_Frobenius(); 993 994 ColumnVector sy = crossproduct(sz, xSun); 995 ColumnVector sx = crossproduct(sy, sz); 996 997 // Effective Dipole of the GPS Satellite Antenna 998 // --------------------------------------------- 999 ColumnVector dipSat = sx - rho * DotProduct(rho,sx) 1000 - crossproduct(rho, sy); 1001 1002 // Receiver unit Vectors rx, ry 1003 // ---------------------------- 1004 ColumnVector rx(3); 1005 ColumnVector ry(3); 1006 1007 double recEll[3]; xyz2ell(rRec.data(), recEll) ; 1008 double neu[3]; 1009 1010 neu[0] = 1.0; 1011 neu[1] = 0.0; 1012 neu[2] = 0.0; 1013 neu2xyz(recEll, neu, rx.data()); 1014 1015 neu[0] = 0.0; 1016 neu[1] = -1.0; 1017 neu[2] = 0.0; 1018 neu2xyz(recEll, neu, ry.data()); 1019 1020 // Effective Dipole of the Receiver Antenna 1021 // ---------------------------------------- 1022 ColumnVector dipRec = rx - rho * DotProduct(rho,rx) 1023 + crossproduct(rho, ry); 1024 1025 // Resulting Effect 1026 // ---------------- 1027 double alpha = DotProduct(dipSat,dipRec) / 1028 (dipSat.norm_Frobenius() * dipRec.norm_Frobenius()); 1029 1030 if (alpha > 1.0) alpha = 1.0; 1031 if (alpha < -1.0) alpha = -1.0; 1032 1033 double dphi = acos(alpha) / 2.0 / M_PI; // in cycles 1034 1035 if ( DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0 ) { 1036 dphi = -dphi; 1037 } 1038 1039 _windUpSum[prn] = floor(_windUpSum[prn] - dphi + 0.5) + dphi; 1040 } 1041 1042 return _windUpSum[prn]; 1043 }
Note:
See TracChangeset
for help on using the changeset viewer.