source: ntrip/trunk/BNC/src/PPP/pppModel.cpp@ 6004

Last change on this file since 6004 was 5828, checked in by mervart, 10 years ago
File size: 11.5 KB
Line 
1
2// Part of BNC, a utility for retrieving decoding and
3// converting GNSS data streams from NTRIP broadcasters.
4//
5// Copyright (C) 2007
6// German Federal Agency for Cartography and Geodesy (BKG)
7// http://www.bkg.bund.de
8// Czech Technical University Prague, Department of Geodesy
9// http://www.fsv.cvut.cz
10//
11// Email: euref-ip@bkg.bund.de
12//
13// This program is free software; you can redistribute it and/or
14// modify it under the terms of the GNU General Public License
15// as published by the Free Software Foundation, version 2.
16//
17// This program is distributed in the hope that it will be useful,
18// but WITHOUT ANY WARRANTY; without even the implied warranty of
19// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20// GNU General Public License for more details.
21//
22// You should have received a copy of the GNU General Public License
23// along with this program; if not, write to the Free Software
24// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
25
26/* -------------------------------------------------------------------------
27 * BKG NTRIP Client
28 * -------------------------------------------------------------------------
29 *
30 * Class: t_astro, t_tides, t_tropo
31 *
32 * Purpose: Observation model
33 *
34 * Author: L. Mervart
35 *
36 * Created: 29-Jul-2014
37 *
38 * Changes:
39 *
40 * -----------------------------------------------------------------------*/
41
42
43#include <cmath>
44
45#include "pppModel.h"
46#include "bncutils.h"
47
48using namespace BNC_PPP;
49using namespace std;
50
51Matrix t_astro::rotX(double Angle) {
52 const double C = cos(Angle);
53 const double S = sin(Angle);
54 Matrix UU(3,3);
55 UU[0][0] = 1.0; UU[0][1] = 0.0; UU[0][2] = 0.0;
56 UU[1][0] = 0.0; UU[1][1] = +C; UU[1][2] = +S;
57 UU[2][0] = 0.0; UU[2][1] = -S; UU[2][2] = +C;
58 return UU;
59}
60
61Matrix t_astro::rotY(double Angle) {
62 const double C = cos(Angle);
63 const double S = sin(Angle);
64 Matrix UU(3,3);
65 UU[0][0] = +C; UU[0][1] = 0.0; UU[0][2] = -S;
66 UU[1][0] = 0.0; UU[1][1] = 1.0; UU[1][2] = 0.0;
67 UU[2][0] = +S; UU[2][1] = 0.0; UU[2][2] = +C;
68 return UU;
69}
70
71Matrix t_astro::rotZ(double Angle) {
72 const double C = cos(Angle);
73 const double S = sin(Angle);
74 Matrix UU(3,3);
75 UU[0][0] = +C; UU[0][1] = +S; UU[0][2] = 0.0;
76 UU[1][0] = -S; UU[1][1] = +C; UU[1][2] = 0.0;
77 UU[2][0] = 0.0; UU[2][1] = 0.0; UU[2][2] = 1.0;
78 return UU;
79}
80
81// Greenwich Mean Sidereal Time
82///////////////////////////////////////////////////////////////////////////
83double t_astro::GMST(double Mjd_UT1) {
84
85 const double Secs = 86400.0;
86
87 double Mjd_0 = floor(Mjd_UT1);
88 double UT1 = Secs*(Mjd_UT1-Mjd_0);
89 double T_0 = (Mjd_0 -MJD_J2000)/36525.0;
90 double T = (Mjd_UT1-MJD_J2000)/36525.0;
91
92 double gmst = 24110.54841 + 8640184.812866*T_0 + 1.002737909350795*UT1
93 + (0.093104-6.2e-6*T)*T*T;
94
95 return 2.0*M_PI*Frac(gmst/Secs);
96}
97
98// Nutation Matrix
99///////////////////////////////////////////////////////////////////////////
100Matrix t_astro::NutMatrix(double Mjd_TT) {
101
102 const double T = (Mjd_TT-MJD_J2000)/36525.0;
103
104 double ls = 2.0*M_PI*Frac(0.993133+ 99.997306*T);
105 double D = 2.0*M_PI*Frac(0.827362+1236.853087*T);
106 double F = 2.0*M_PI*Frac(0.259089+1342.227826*T);
107 double N = 2.0*M_PI*Frac(0.347346- 5.372447*T);
108
109 double dpsi = ( -17.200*sin(N) - 1.319*sin(2*(F-D+N)) - 0.227*sin(2*(F+N))
110 + 0.206*sin(2*N) + 0.143*sin(ls) ) / RHO_SEC;
111 double deps = ( + 9.203*cos(N) + 0.574*cos(2*(F-D+N)) + 0.098*cos(2*(F+N))
112 - 0.090*cos(2*N) ) / RHO_SEC;
113
114 double eps = 0.4090928-2.2696E-4*T;
115
116 return rotX(-eps-deps)*rotZ(-dpsi)*rotX(+eps);
117}
118
119// Precession Matrix
120///////////////////////////////////////////////////////////////////////////
121Matrix t_astro::PrecMatrix(double Mjd_1, double Mjd_2) {
122
123 const double T = (Mjd_1-MJD_J2000)/36525.0;
124 const double dT = (Mjd_2-Mjd_1)/36525.0;
125
126 double zeta = ( (2306.2181+(1.39656-0.000139*T)*T)+
127 ((0.30188-0.000344*T)+0.017998*dT)*dT )*dT/RHO_SEC;
128 double z = zeta + ( (0.79280+0.000411*T)+0.000205*dT)*dT*dT/RHO_SEC;
129 double theta = ( (2004.3109-(0.85330+0.000217*T)*T)-
130 ((0.42665+0.000217*T)+0.041833*dT)*dT )*dT/RHO_SEC;
131
132 return rotZ(-z) * rotY(theta) * rotZ(-zeta);
133}
134
135// Sun's position
136///////////////////////////////////////////////////////////////////////////
137ColumnVector t_astro::Sun(double Mjd_TT) {
138
139 const double eps = 23.43929111/RHO_DEG;
140 const double T = (Mjd_TT-MJD_J2000)/36525.0;
141
142 double M = 2.0*M_PI * Frac ( 0.9931267 + 99.9973583*T);
143 double L = 2.0*M_PI * Frac ( 0.7859444 + M/2.0/M_PI +
144 (6892.0*sin(M)+72.0*sin(2.0*M)) / 1296.0e3);
145 double r = 149.619e9 - 2.499e9*cos(M) - 0.021e9*cos(2*M);
146
147 ColumnVector r_Sun(3);
148 r_Sun << r*cos(L) << r*sin(L) << 0.0; r_Sun = rotX(-eps) * r_Sun;
149
150 return rotZ(GMST(Mjd_TT))
151 * NutMatrix(Mjd_TT)
152 * PrecMatrix(MJD_J2000, Mjd_TT)
153 * r_Sun;
154}
155
156// Moon's position
157///////////////////////////////////////////////////////////////////////////
158ColumnVector t_astro::Moon(double Mjd_TT) {
159
160 const double eps = 23.43929111/RHO_DEG;
161 const double T = (Mjd_TT-MJD_J2000)/36525.0;
162
163 double L_0 = Frac ( 0.606433 + 1336.851344*T );
164 double l = 2.0*M_PI*Frac ( 0.374897 + 1325.552410*T );
165 double lp = 2.0*M_PI*Frac ( 0.993133 + 99.997361*T );
166 double D = 2.0*M_PI*Frac ( 0.827361 + 1236.853086*T );
167 double F = 2.0*M_PI*Frac ( 0.259086 + 1342.227825*T );
168
169 double dL = +22640*sin(l) - 4586*sin(l-2*D) + 2370*sin(2*D) + 769*sin(2*l)
170 -668*sin(lp) - 412*sin(2*F) - 212*sin(2*l-2*D)- 206*sin(l+lp-2*D)
171 +192*sin(l+2*D) - 165*sin(lp-2*D) - 125*sin(D) - 110*sin(l+lp)
172 +148*sin(l-lp) - 55*sin(2*F-2*D);
173
174 double L = 2.0*M_PI * Frac( L_0 + dL/1296.0e3 );
175
176 double S = F + (dL+412*sin(2*F)+541*sin(lp)) / RHO_SEC;
177 double h = F-2*D;
178 double N = -526*sin(h) + 44*sin(l+h) - 31*sin(-l+h) - 23*sin(lp+h)
179 +11*sin(-lp+h) - 25*sin(-2*l+F) + 21*sin(-l+F);
180
181 double B = ( 18520.0*sin(S) + N ) / RHO_SEC;
182
183 double cosB = cos(B);
184
185 double R = 385000e3 - 20905e3*cos(l) - 3699e3*cos(2*D-l) - 2956e3*cos(2*D)
186 -570e3*cos(2*l) + 246e3*cos(2*l-2*D) - 205e3*cos(lp-2*D)
187 -171e3*cos(l+2*D) - 152e3*cos(l+lp-2*D);
188
189 ColumnVector r_Moon(3);
190 r_Moon << R*cos(L)*cosB << R*sin(L)*cosB << R*sin(B);
191 r_Moon = rotX(-eps) * r_Moon;
192
193 return rotZ(GMST(Mjd_TT))
194 * NutMatrix(Mjd_TT)
195 * PrecMatrix(MJD_J2000, Mjd_TT)
196 * r_Moon;
197}
198
199// Tidal Correction
200////////////////////////////////////////////////////////////////////////////
201ColumnVector t_tides::displacement(const bncTime& time, const ColumnVector& xyz) {
202
203 double Mjd = time.mjd() + time.daysec() / 86400.0;
204
205 if (Mjd != _lastMjd) {
206 _lastMjd = Mjd;
207 _xSun = t_astro::Sun(Mjd);
208 _rSun = sqrt(DotProduct(_xSun,_xSun));
209 _xSun /= _rSun;
210 _xMoon = t_astro::Moon(Mjd);
211 _rMoon = sqrt(DotProduct(_xMoon,_xMoon));
212 _xMoon /= _rMoon;
213 }
214
215 double rRec = sqrt(DotProduct(xyz, xyz));
216 ColumnVector xyzUnit = xyz / rRec;
217
218 // Love's Numbers
219 // --------------
220 const double H2 = 0.6078;
221 const double L2 = 0.0847;
222
223 // Tidal Displacement
224 // ------------------
225 double scSun = DotProduct(xyzUnit, _xSun);
226 double scMoon = DotProduct(xyzUnit, _xMoon);
227
228 double p2Sun = 3.0 * (H2/2.0-L2) * scSun * scSun - H2/2.0;
229 double p2Moon = 3.0 * (H2/2.0-L2) * scMoon * scMoon - H2/2.0;
230
231 double x2Sun = 3.0 * L2 * scSun;
232 double x2Moon = 3.0 * L2 * scMoon;
233
234 const double gmWGS = 398.6005e12;
235 const double gms = 1.3271250e20;
236 const double gmm = 4.9027890e12;
237
238 double facSun = gms / gmWGS *
239 (rRec * rRec * rRec * rRec) / (_rSun * _rSun * _rSun);
240
241 double facMoon = gmm / gmWGS *
242 (rRec * rRec * rRec * rRec) / (_rMoon * _rMoon * _rMoon);
243
244 ColumnVector dX = facSun * (x2Sun * _xSun + p2Sun * xyzUnit) +
245 facMoon * (x2Moon * _xMoon + p2Moon * xyzUnit);
246
247 return dX;
248}
249
250// Constructor
251///////////////////////////////////////////////////////////////////////////
252t_windUp::t_windUp() {
253 for (unsigned ii = 0; ii <= t_prn::MAXPRN; ii++) {
254 sumWind[ii] = 0.0;
255 lastEtime[ii] = 0.0;
256 }
257}
258
259// Phase Wind-Up Correction
260///////////////////////////////////////////////////////////////////////////
261double t_windUp::value(const bncTime& etime, const ColumnVector& rRec,
262 t_prn prn, const ColumnVector& rSat) {
263
264 if (etime.mjddec() != lastEtime[prn.toInt()]) {
265
266 // Unit Vector GPS Satellite --> Receiver
267 // --------------------------------------
268 ColumnVector rho = rRec - rSat;
269 rho /= rho.norm_Frobenius();
270
271 // GPS Satellite unit Vectors sz, sy, sx
272 // -------------------------------------
273 ColumnVector sz = -rSat / rSat.norm_Frobenius();
274
275 ColumnVector xSun = t_astro::Sun(etime.mjddec());
276 xSun /= xSun.norm_Frobenius();
277
278 ColumnVector sy = crossproduct(sz, xSun);
279 ColumnVector sx = crossproduct(sy, sz);
280
281 // Effective Dipole of the GPS Satellite Antenna
282 // ---------------------------------------------
283 ColumnVector dipSat = sx - rho * DotProduct(rho,sx)
284 - crossproduct(rho, sy);
285
286 // Receiver unit Vectors rx, ry
287 // ----------------------------
288 ColumnVector rx(3);
289 ColumnVector ry(3);
290
291 double recEll[3]; xyz2ell(rRec.data(), recEll) ;
292 double neu[3];
293
294 neu[0] = 1.0;
295 neu[1] = 0.0;
296 neu[2] = 0.0;
297 neu2xyz(recEll, neu, rx.data());
298
299 neu[0] = 0.0;
300 neu[1] = -1.0;
301 neu[2] = 0.0;
302 neu2xyz(recEll, neu, ry.data());
303
304 // Effective Dipole of the Receiver Antenna
305 // ----------------------------------------
306 ColumnVector dipRec = rx - rho * DotProduct(rho,rx)
307 + crossproduct(rho, ry);
308
309 // Resulting Effect
310 // ----------------
311 double alpha = DotProduct(dipSat,dipRec) /
312 (dipSat.norm_Frobenius() * dipRec.norm_Frobenius());
313
314 if (alpha > 1.0) alpha = 1.0;
315 if (alpha < -1.0) alpha = -1.0;
316
317 double dphi = acos(alpha) / 2.0 / M_PI; // in cycles
318
319 if ( DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0 ) {
320 dphi = -dphi;
321 }
322
323 if (lastEtime[prn.toInt()] == 0.0) {
324 sumWind[prn.toInt()] = dphi;
325 }
326 else {
327 sumWind[prn.toInt()] = nint(sumWind[prn.toInt()] - dphi) + dphi;
328 }
329
330 lastEtime[prn.toInt()] = etime.mjddec();
331 }
332
333 return sumWind[prn.toInt()];
334}
335
336// Tropospheric Model (Saastamoinen)
337////////////////////////////////////////////////////////////////////////////
338double t_tropo::delay_saast(const ColumnVector& xyz, double Ele) {
339
340 Tracer tracer("bncModel::delay_saast");
341
342 if (xyz[0] == 0.0 && xyz[1] == 0.0 && xyz[2] == 0.0) {
343 return 0.0;
344 }
345
346 double ell[3];
347 xyz2ell(xyz.data(), ell);
348 double height = ell[2];
349
350 double pp = 1013.25 * pow(1.0 - 2.26e-5 * height, 5.225);
351 double TT = 18.0 - height * 0.0065 + 273.15;
352 double hh = 50.0 * exp(-6.396e-4 * height);
353 double ee = hh / 100.0 * exp(-37.2465 + 0.213166*TT - 0.000256908*TT*TT);
354
355 double h_km = height / 1000.0;
356
357 if (h_km < 0.0) h_km = 0.0;
358 if (h_km > 5.0) h_km = 5.0;
359 int ii = int(h_km + 1);
360 double href = ii - 1;
361
362 double bCor[6];
363 bCor[0] = 1.156;
364 bCor[1] = 1.006;
365 bCor[2] = 0.874;
366 bCor[3] = 0.757;
367 bCor[4] = 0.654;
368 bCor[5] = 0.563;
369
370 double BB = bCor[ii-1] + (bCor[ii]-bCor[ii-1]) * (h_km - href);
371
372 double zen = M_PI/2.0 - Ele;
373
374 return (0.002277/cos(zen)) * (pp + ((1255.0/TT)+0.05)*ee - BB*(tan(zen)*tan(zen)));
375}
376
Note: See TracBrowser for help on using the repository browser.