 Timestamp:
 Nov 25, 2020, 10:03:12 PM (3 years ago)
 Location:
 branches/BNC_2.12/src
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

branches/BNC_2.12/src/PPP_SSR_I/pppFilter.cpp
r9086 r9281 409 409 xyz2ell(xyz, ell); 410 410 double height = ell[2]; 411 // Prevent pp from causing segmentation fault (Loukis) 412 if (height > 40000.0 ) { 413 return 0.000000001; 414 } 411 415 412 416 double pp = 1013.25 * pow(1.0  2.26e5 * height, 5.225); 
branches/BNC_2.12/src/bncwindow.cpp
r9251 r9281 188 188 _onTheFlyComboBox = new QComboBox(); 189 189 _onTheFlyComboBox>setEditable(false); 190 _onTheFlyComboBox>addItems(QString("no, 190 _onTheFlyComboBox>addItems(QString("no,1 day,1 hour,5 min,1 min").split(",")); 191 191 int ii = _onTheFlyComboBox>findText(settings.value("onTheFlyInterval").toString()); 192 192 if (ii != 1) { 
branches/BNC_2.12/src/pppModel.cpp
r7259 r9281 36 36 * Created: 29Jul2014 37 37 * 38 * Changes: 38 * Changes: 39 39 * 40 40 * */ … … 90 90 double Mjd_0 = floor(Mjd_UT1); 91 91 double UT1 = Secs*(Mjd_UT1Mjd_0); 92 double T_0 = (Mjd_0 MJD_J2000)/36525.0; 93 double T = (Mjd_UT1MJD_J2000)/36525.0; 92 double T_0 = (Mjd_0 MJD_J2000)/36525.0; 93 double T = (Mjd_UT1MJD_J2000)/36525.0; 94 94 95 95 double gmst = 24110.54841 + 8640184.812866*T_0 + 1.002737909350795*UT1 … … 126 126 const double T = (Mjd_1MJD_J2000)/36525.0; 127 127 const double dT = (Mjd_2Mjd_1)/36525.0; 128 128 129 129 double zeta = ( (2306.2181+(1.396560.000139*T)*T)+ 130 130 ((0.301880.000344*T)+0.017998*dT)*dT )*dT/RHO_SEC; … … 134 134 135 135 return rotZ(z) * rotY(theta) * rotZ(zeta); 136 } 136 } 137 137 138 138 // Sun's position … … 144 144 145 145 double M = 2.0*M_PI * Frac ( 0.9931267 + 99.9973583*T); 146 double L = 2.0*M_PI * Frac ( 0.7859444 + M/2.0/M_PI + 146 double L = 2.0*M_PI * Frac ( 0.7859444 + M/2.0/M_PI + 147 147 (6892.0*sin(M)+72.0*sin(2.0*M)) / 1296.0e3); 148 148 double r = 149.619e9  2.499e9*cos(M)  0.021e9*cos(2*M); 149 150 ColumnVector r_Sun(3); 149 150 ColumnVector r_Sun(3); 151 151 r_Sun << r*cos(L) << r*sin(L) << 0.0; r_Sun = rotX(eps) * r_Sun; 152 152 153 153 return rotZ(GMST(Mjd_TT)) 154 * NutMatrix(Mjd_TT) 154 * NutMatrix(Mjd_TT) 155 155 * PrecMatrix(MJD_J2000, Mjd_TT) 156 156 * r_Sun; … … 169 169 double D = 2.0*M_PI*Frac ( 0.827361 + 1236.853086*T ); 170 170 double F = 2.0*M_PI*Frac ( 0.259086 + 1342.227825*T ); 171 172 double dL = +22640*sin(l)  4586*sin(l2*D) + 2370*sin(2*D) + 769*sin(2*l) 171 172 double dL = +22640*sin(l)  4586*sin(l2*D) + 2370*sin(2*D) + 769*sin(2*l) 173 173 668*sin(lp)  412*sin(2*F)  212*sin(2*l2*D) 206*sin(l+lp2*D) 174 174 +192*sin(l+2*D)  165*sin(lp2*D)  125*sin(D)  110*sin(l+lp) … … 177 177 double L = 2.0*M_PI * Frac( L_0 + dL/1296.0e3 ); 178 178 179 double S = F + (dL+412*sin(2*F)+541*sin(lp)) / RHO_SEC; 179 double S = F + (dL+412*sin(2*F)+541*sin(lp)) / RHO_SEC; 180 180 double h = F2*D; 181 double N = 526*sin(h) + 44*sin(l+h)  31*sin(l+h)  23*sin(lp+h) 181 double N = 526*sin(h) + 44*sin(l+h)  31*sin(l+h)  23*sin(lp+h) 182 182 +11*sin(lp+h)  25*sin(2*l+F) + 21*sin(l+F); 183 183 184 184 double B = ( 18520.0*sin(S) + N ) / RHO_SEC; 185 185 186 186 double cosB = cos(B); 187 187 188 188 double R = 385000e3  20905e3*cos(l)  3699e3*cos(2*Dl)  2956e3*cos(2*D) 189 570e3*cos(2*l) + 246e3*cos(2*l2*D)  205e3*cos(lp2*D) 190 171e3*cos(l+2*D)  152e3*cos(l+lp2*D); 191 192 ColumnVector r_Moon(3); 189 570e3*cos(2*l) + 246e3*cos(2*l2*D)  205e3*cos(lp2*D) 190 171e3*cos(l+2*D)  152e3*cos(l+lp2*D); 191 192 ColumnVector r_Moon(3); 193 193 r_Moon << R*cos(L)*cosB << R*sin(L)*cosB << R*sin(B); 194 194 r_Moon = rotX(eps) * r_Moon; 195 196 return rotZ(GMST(Mjd_TT)) 197 * NutMatrix(Mjd_TT) 195 196 return rotZ(GMST(Mjd_TT)) 197 * NutMatrix(Mjd_TT) 198 198 * PrecMatrix(MJD_J2000, Mjd_TT) 199 199 * r_Moon; 200 200 } 201 201 202 // Tidal Correction 202 // Tidal Correction 203 203 //////////////////////////////////////////////////////////////////////////// 204 204 ColumnVector t_tides::displacement(const bncTime& time, const ColumnVector& xyz) { … … 239 239 double x2Sun = 3.0 * L2 * scSun; 240 240 double x2Moon = 3.0 * L2 * scMoon; 241 241 242 242 const double gmWGS = 398.6005e12; 243 243 const double gms = 1.3271250e20; 244 244 const double gmm = 4.9027890e12; 245 245 246 double facSun = gms / gmWGS * 246 double facSun = gms / gmWGS * 247 247 (rRec * rRec * rRec * rRec) / (_rSun * _rSun * _rSun); 248 248 249 double facMoon = gmm / gmWGS * 249 double facMoon = gmm / gmWGS * 250 250 (rRec * rRec * rRec * rRec) / (_rMoon * _rMoon * _rMoon); 251 251 252 ColumnVector dX = facSun * (x2Sun * _xSun + p2Sun * xyzUnit) + 252 ColumnVector dX = facSun * (x2Sun * _xSun + p2Sun * xyzUnit) + 253 253 facMoon * (x2Moon * _xMoon + p2Moon * xyzUnit); 254 254 … … 267 267 // Phase WindUp Correction 268 268 /////////////////////////////////////////////////////////////////////////// 269 double t_windUp::value(const bncTime& etime, const ColumnVector& rRec, 269 double t_windUp::value(const bncTime& etime, const ColumnVector& rRec, 270 270 t_prn prn, const ColumnVector& rSat) { 271 271 … … 276 276 ColumnVector rho = rRec  rSat; 277 277 rho /= rho.norm_Frobenius(); 278 278 279 279 // GPS Satellite unit Vectors sz, sy, sx 280 280 //  … … 289 289 // Effective Dipole of the GPS Satellite Antenna 290 290 //  291 ColumnVector dipSat = sx  rho * DotProduct(rho,sx) 291 ColumnVector dipSat = sx  rho * DotProduct(rho,sx) 292 292  crossproduct(rho, sy); 293 293 294 294 // Receiver unit Vectors rx, ry 295 295 //  … … 299 299 double recEll[3]; xyz2ell(rRec.data(), recEll) ; 300 300 double neu[3]; 301 301 302 302 neu[0] = 1.0; 303 303 neu[1] = 0.0; 304 304 neu[2] = 0.0; 305 305 neu2xyz(recEll, neu, rx.data()); 306 306 307 307 neu[0] = 0.0; 308 308 neu[1] = 1.0; 309 309 neu[2] = 0.0; 310 310 neu2xyz(recEll, neu, ry.data()); 311 311 312 312 // Effective Dipole of the Receiver Antenna 313 313 //  314 ColumnVector dipRec = rx  rho * DotProduct(rho,rx) 314 ColumnVector dipRec = rx  rho * DotProduct(rho,rx) 315 315 + crossproduct(rho, ry); 316 316 317 317 // Resulting Effect 318 318 //  319 double alpha = DotProduct(dipSat,dipRec) / 319 double alpha = DotProduct(dipSat,dipRec) / 320 320 (dipSat.norm_Frobenius() * dipRec.norm_Frobenius()); 321 321 322 322 if (alpha > 1.0) alpha = 1.0; 323 323 if (alpha < 1.0) alpha = 1.0; 324 324 325 325 double dphi = acos(alpha) / 2.0 / M_PI; // in cycles 326 326 327 327 if ( DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0 ) { 328 328 dphi = dphi; … … 339 339 } 340 340 341 return sumWind[prn.toInt()]; 341 return sumWind[prn.toInt()]; 342 342 } 343 343 … … 352 352 } 353 353 354 double ell[3]; 354 double ell[3]; 355 355 xyz2ell(xyz.data(), ell); 356 356 double height = ell[2]; 357 // Prevent pp from causing segmentation fault (Loukis) 358 if (height > 40000.0 ) { 359 return 0.000000001; 360 } 357 361 358 362 double pp = 1013.25 * pow(1.0  2.26e5 * height, 5.225); … … 362 366 363 367 double h_km = height / 1000.0; 364 368 365 369 if (h_km < 0.0) h_km = 0.0; 366 370 if (h_km > 5.0) h_km = 5.0; 367 371 int ii = int(h_km + 1); if (ii > 5) ii = 5; 368 372 double href = ii  1; 369 370 double bCor[6]; 373 374 double bCor[6]; 371 375 bCor[0] = 1.156; 372 376 bCor[1] = 1.006; … … 375 379 bCor[4] = 0.654; 376 380 bCor[5] = 0.563; 377 381 378 382 double BB = bCor[ii1] + (bCor[ii]bCor[ii1]) * (h_km  href); 379 383 380 384 double zen = M_PI/2.0  Ele; 381 385
Note:
See TracChangeset
for help on using the changeset viewer.