1 //============================================================================== 2 // 3 // This file is part of GPSTk, the GPS Toolkit. 4 // 5 // The GPSTk is free software; you can redistribute it and/or modify 6 // it under the terms of the GNU Lesser General Public License as published 7 // by the Free Software Foundation; either version 3.0 of the License, or 8 // any later version. 9 // 10 // The GPSTk is distributed in the hope that it will be useful, 11 // but WITHOUT ANY WARRANTY; without even the implied warranty of 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 // GNU Lesser General Public License for more details. 14 // 15 // You should have received a copy of the GNU Lesser General Public 16 // License along with GPSTk; if not, write to the Free Software Foundation, 17 // Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA 18 // 19 // This software was developed by Applied Research Laboratories at the 20 // University of Texas at Austin. 21 // Copyright 2004-2020, The Board of Regents of The University of Texas System 22 // 23 //============================================================================== 24 25 //============================================================================== 26 // 27 // This software was developed by Applied Research Laboratories at the 28 // University of Texas at Austin, under contract to an agency or agencies 29 // within the U.S. Department of Defense. The U.S. Government retains all 30 // rights to use, duplicate, distribute, disclose, or release this software. 31 // 32 // Pursuant to DoD Directive 523024 33 // 34 // DISTRIBUTION STATEMENT A: This software has been approved for public 35 // release, distribution is unlimited. 36 // 37 //============================================================================== 38 39 /// @file SolidEarthTides.cpp 40 /// Implement the formula for the displacement of a point fixed to the solid Earth 41 /// due to the solid Earth tides resulting from the influence of the Sun and Moon. 42 /// Reference IERS Conventions (1996) found in IERS Technical Note 21 (IERS). 43 /// NB. Currently only the largest terms are implemented, yielding a result accurate 44 /// to the millimeter level. Specifically, IERS pg 61 eq 8 and IERS pg 65 eq 17, 45 /// (this includes removing the permanent component). 46 /// Class SolarSystem may be used to get Solar and Lunar ephemeris information, 47 /// including position and mass ratios. 48 49 //------------------------------------------------------------------------------------ 50 // GPSTk 51 //#include "geometry.hpp" // for DEG_TO_RAD 52 // geomatics 53 #include "logstream.hpp" 54 #include "SolidEarthTides.hpp" 55 56 using namespace std; 57 58 namespace gpstk 59 { 60 //--------------------------------------------------------------------------------- 61 // Compute the site displacement due to solid Earth tides for the given Position 62 // (assumed to be fixed to the solid Earth) at the given time, given the position 63 // of the site of interest, positions and mass ratios of the sun and moon. 64 // Return a Triple containing the site displacement in ECEF XYZ coordinates with 65 // units meters. 66 // Reference IERS Conventions (1996) found in IERS Technical Note 21 67 // and IERS Conventions (2003) found in IERS Technical Note 32 68 // and IERS Conventions (2010) found in IERS Technical Note 36. 69 // NB. Currently only the largest terms are implemented, yielding a result 70 // accurate to the millimeter level. Specifically, TN21 pg 61 eq 8 and 71 // TN21 pg 65 eq 17. 72 // param site Position Nominal position of the site of interest. 73 // param ttag EphTime Time of interest. 74 // param Sun Position Position of the Sun at time 75 // param Moon Position Position of the Moon at time 76 // param EMRAT double Earth-to-Moon mass ratio (default to DE405 value) 77 // param SERAT double Sun-to-Earth mass ratio (default to DE405 value) 78 // param IERSConvention IERS convention to use (default IERS2010) 79 // return Triple Displacement vector, ECEF XYZ in meters. computeSolidEarthTides(const Position site,const EphTime ttag,const Position Sun,const Position Moon,const double EMRAT,const double SERAT,const IERSConvention iers)80 Triple computeSolidEarthTides(const Position site, 81 const EphTime ttag, 82 const Position Sun, 83 const Position Moon, 84 const double EMRAT, 85 const double SERAT, 86 const IERSConvention iers) 87 { 88 try { 89 //if(ttag.getTimeSystem() == TimeSystem::Unknown) { 90 // Exception e("Time system is unknown"); 91 // GPSTK_THROW(e); 92 //} 93 94 // Use REarth from solid.f example program 95 static const double REarth=6378136.55; 96 static const bool debug = (LOGlevel >= DEBUG7); 97 // NB icount is a dummy used in test vs solid.f 98 int i,icount(-1); 99 double RSun, RMoon, Rx, Love, Shida, sunFactor, moonFactor; 100 double sunDOTrx, moonDOTrx, REoRS, REoRM; 101 double lat, lon, sinlat, coslat, sinlon, coslon; 102 double latSun, lonSun, latMoon, lonMoon; 103 Triple disp, sunUnit, moonUnit, rx, tSun, tMoon, north, east, up; 104 // quantities for debug printing only 105 Triple northGD, eastGD, upGD, tmp, tmp2, tmp3, tmp4; 106 107 LOG(DEBUG7) << "Sun position " << ttag.asGPSString() 108 << fixed << setprecision(3) 109 << setw(23) << Sun.X() << setw(23) << Sun.Y() << setw(23) << Sun.Z(); 110 LOG(DEBUG7) << "Moon position" << ttag.asGPSString() 111 << fixed << setprecision(3) 112 << setw(23) << Moon.X() << setw(23) << Moon.Y() << setw(23) << Moon.Z(); 113 114 // distances (m) 115 RSun = Sun.radius(); 116 RMoon = Moon.radius(); 117 Rx = site.radius(); 118 119 // unit vectors 120 sunUnit = Triple(Sun.X()/RSun, Sun.Y()/RSun, Sun.Z()/RSun); 121 moonUnit = Triple(Moon.X()/RMoon, Moon.Y()/RMoon, Moon.Z()/RMoon); 122 rx = Triple(site.X()/Rx, site.Y()/Rx, site.Z()/Rx); 123 124 // generate geodetic transformation first - for debug 125 if(debug) { 126 lat = site.getGeodeticLatitude()*DEG_TO_RAD; 127 lon = site.getLongitude()*DEG_TO_RAD; 128 sinlat = ::sin(lat); 129 coslat = ::cos(lat); 130 sinlon = ::sin(lon); 131 coslon = ::cos(lon); 132 133 // transform X=(x,y,z) into (R*X)(north,east,up) using geodetic longitude 134 northGD = Triple(-sinlat*coslon, -sinlat*sinlon, coslat); 135 eastGD = Triple( -sinlon, coslon, 0.0); 136 upGD = Triple( coslat*coslon, coslat*sinlon, sinlat); 137 } 138 139 // use geocentric latitude for formulas 140 latSun = Sun.getGeocentricLatitude()*DEG_TO_RAD; 141 lonSun = Sun.getLongitude()*DEG_TO_RAD; 142 latMoon = Moon.getGeocentricLatitude()*DEG_TO_RAD; 143 lonMoon = Moon.getLongitude()*DEG_TO_RAD; 144 lat = site.getGeocentricLatitude()*DEG_TO_RAD; 145 lon = site.getLongitude()*DEG_TO_RAD; 146 sinlat = ::sin(lat); 147 coslat = ::cos(lat); 148 sinlon = ::sin(lon); 149 coslon = ::cos(lon); 150 151 // transform X=(x,y,z) into (R*X)(north,east,up) using geocentric longitude 152 north = Triple(-sinlat*coslon, -sinlat*sinlon, coslat); 153 east = Triple( -sinlon, coslon, 0.0); 154 up = Triple( coslat*coslon, coslat*sinlon, sinlat); 155 156 // GM*R factors 157 REoRS = REarth/RSun; // ratio Earth/Sun radius 158 sunFactor = REarth*REoRS*REoRS*REoRS*SERAT; // = (GMS/GME)*RE^4/RS^3 159 REoRM = REarth/RMoon; // ratio Earth/Moon radius 160 moonFactor = REarth*REoRM*REoRM*REoRM/EMRAT; // = (GMM/GME)*RE^4/RM^3 161 // E/M mass ratio (403) 81.300584999999998 (405) 81.300560000000004 162 // S/E mass ratio (403) 332946.048630181234330 (405) 332946.050894783285912 163 //LOG(INFO) << " E/M mass ratio " << fixed << setprecision(15) << EMRAT; 164 //LOG(INFO) << " S/E mass ratio " << fixed << setprecision(15) << SERAT; 165 166 // dot products 167 sunDOTrx = sunUnit.dot(rx); 168 moonDOTrx = moonUnit.dot(rx); 169 170 // transverse to radial direction - not unit vectors 171 tSun = sunUnit - sunDOTrx * rx; 172 tMoon = moonUnit - moonDOTrx * rx; 173 174 // ------------------------------------------------------------------------- 175 // compute displacements 176 disp = Triple(0,0,0); 177 // Steps and equations refer to IERS(1996), esp table on page 60; 178 // formulas are generally repeated in other IERS technical notes. 179 180 // Step 1a IERS(1996) eq. (8) pg 61. 181 // nominal degree 2 Love and Shida numbers pg 60 182 double poly = sinlat; 183 poly = (3.0*poly*poly-1.0)/2.0; 184 185 // here is the only difference between 1996 and 2003/10 186 if(iers == IERSConvention::IERS1996) { 187 Love = 0.6026 - 0.0006*poly; 188 Shida = 0.0831 + 0.0002*poly; 189 } 190 else { // 2003 or 2010 191 Love = 0.6078 - 0.0006*poly; 192 Shida = 0.0847 + 0.0002*poly; 193 } 194 LOG(DEBUG6) << "H2L2 " << setw(4) << icount << fixed << setprecision(15) 195 << " " << setw(18) << Love 196 << " " << setw(18) << Shida 197 << " " << setw(18) << poly; 198 LOG(DEBUG6) << "P2 " << setw(4) << icount << fixed << setprecision(15) 199 << " " << setw(18) << 3*(Love/2-Shida)*sunDOTrx*sunDOTrx-0.5*Love 200 << " " << setw(18) << 3*(Love/2-Shida)*moonDOTrx*moonDOTrx-0.5*Love; 201 202 tmp = sunFactor * (Love * (1.5*sunDOTrx*sunDOTrx-0.5) * rx 203 + 3.0*Shida*sunDOTrx*tSun); 204 for(i=0; i<3; i++) disp[i] += tmp[i]; 205 206 tmp = moonFactor * (Love * (1.5*moonDOTrx*moonDOTrx-0.5) * rx 207 + 3.0*Shida*moonDOTrx*tMoon); 208 for(i=0; i<3; i++) disp[i] += tmp[i]; 209 210 // Step 1b IERS(1996) eq. (9) pg 61. 211 // nominal degree 3 Love and Shida numbers pg 60 212 double Shida2=Shida; 213 Love = 0.292; 214 Shida = 0.015; 215 tmp = sunFactor*REoRS * (Love * (2.5*sunDOTrx*sunDOTrx - 1.5) * sunDOTrx * rx 216 + Shida * (7.5*sunDOTrx*sunDOTrx - 1.5) * tSun); 217 for(i=0; i<3; i++) disp[i] += tmp[i]; 218 219 LOG(DEBUG6) << "P3 " << setw(4) << icount << fixed << setprecision(15) 220 << " " << setw(18) << 2.5*(Love-3*Shida)*sunDOTrx*sunDOTrx*sunDOTrx 221 +1.5*(Shida-Love)*sunDOTrx 222 << " " << setw(18) << 2.5*(Love-3*Shida)*moonDOTrx*moonDOTrx*moonDOTrx 223 +1.5*(Shida-Love)*moonDOTrx; 224 LOG(DEBUG6) << "X2 " << setw(4) << icount << fixed << setprecision(15) 225 << " " << setw(18) << 3*Shida2*sunDOTrx 226 << " " << setw(18) << 3*Shida2*moonDOTrx; 227 LOG(DEBUG6) << "X3 " << setw(4) << icount << fixed << setprecision(15) 228 << " " << setw(18) << 1.5*Shida*(5*sunDOTrx*sunDOTrx-1) 229 << " " << setw(18) << 1.5*Shida*(5*moonDOTrx*moonDOTrx-1); 230 LOG(DEBUG6) << "RAT " << setw(4) << icount << fixed << setprecision(6) 231 << " " << setw(18) << SERAT << setprecision(15) 232 << " " << setw(22) << EMRAT << setprecision(2) 233 << " " << setw(11) << REarth; 234 LOG(DEBUG6) << "FACT2 " << setw(4) << icount << fixed << setprecision(15) 235 << " " << setw(18) << sunFactor 236 << " " << setw(18) << moonFactor; 237 LOG(DEBUG6) << "FACT3 " << setw(4) << icount << fixed << setprecision(15) 238 << " " << setw(18) << sunFactor*REoRS 239 << " " << setw(18) << moonFactor*REoRM; 240 241 tmp = moonFactor*REoRM * (Love * (2.5*moonDOTrx*moonDOTrx-1.5) * moonDOTrx * rx 242 + Shida * (7.5*moonDOTrx*moonDOTrx-1.5) * tMoon); 243 for(i=0; i<3; i++) disp[i] += tmp[i]; 244 245 // all of 8 and 9 246 if(debug) { 247 tmp2[0] = northGD[0]*disp[0] + northGD[1]*disp[1] + northGD[2]*disp[2]; 248 tmp2[1] = eastGD[0]*disp[0] + eastGD[1]*disp[1] + eastGD[2]*disp[2]; 249 tmp2[2] = upGD[0]*disp[0] + upGD[1]*disp[1] + upGD[2]*disp[2]; 250 LOG(DEBUG7) << "7SET solar/lunar/2nd/3rd " << ttag.asGPSString() 251 << fixed << setprecision(9) 252 << " " << disp[0] << " " << disp[1] << " " << disp[2] 253 << " " << tmp2[0] << " " << tmp2[1] << " " << tmp2[2]; 254 } 255 LOG(DEBUG6) << "DX0 " << setw(4) << icount << fixed << setprecision(15) 256 << " " << setw(18) << disp[0] 257 << " " << setw(18) << disp[1] 258 << " " << setw(18) << disp[2]; 259 260 // Step 1c IERS(1996) eq. (13) pg 63. diurnal tides 261 Love = -0.0025; 262 Shida = -0.0007; 263 tmp = -0.75*Love*::sin(2*lat) * // radial 13a 264 ( sunFactor*::sin(2*latSun)*::sin(lon-lonSun) 265 + moonFactor*::sin(2*latMoon)*::sin(lon-lonMoon)) * rx 266 267 - 1.5*Shida*::cos(2*lat) * // north 13b 268 ( sunFactor*::sin(2*latSun)*::sin(lon-lonSun) 269 + moonFactor*::sin(2*latMoon)*::sin(lon-lonMoon)) * north 270 271 - 1.5*Shida*sinlat * // east 13b 272 ( sunFactor*::sin(2*latSun)*::cos(lon-lonSun) 273 + moonFactor*::sin(2*latMoon)*::cos(lon-lonMoon)) * east; 274 275 LOG(DEBUG6) << "DX1 " << setw(4) << icount << fixed << setprecision(15) 276 << " " << setw(18) << tmp[0] 277 << " " << setw(18) << tmp[1] 278 << " " << setw(18) << tmp[2]; 279 280 for(i=0; i<3; i++) disp[i] += tmp[i]; 281 282 if(debug) { 283 tmp2[0] = northGD[0]*tmp[0] + northGD[1]*tmp[1] + northGD[2]*tmp[2]; 284 tmp2[1] = eastGD[0]*tmp[0] + eastGD[1]*tmp[1] + eastGD[2]*tmp[2]; 285 tmp2[2] = upGD[0]*tmp[0] + upGD[1]*tmp[1] + upGD[2]*tmp[2]; 286 LOG(DEBUG7) << "7SET diurnal-band " << ttag.asGPSString() 287 << fixed << setprecision(9) 288 << " " << tmp[0] << " " << tmp[1] << " " << tmp[2] 289 << " " << tmp2[0] << " " << tmp2[1] << " " << tmp2[2]; 290 } 291 292 // Step 1d IERS(1996) eq. (14) pg 63. semidiurnal tides 293 Love = -0.0022; 294 Shida = -0.0007; 295 tmp = -0.75*Love*coslat*coslat * // radial 14a 296 ( sunFactor*::cos(latSun)*::cos(latSun)*::sin(2*(lon-lonSun)) 297 + moonFactor*::cos(latMoon)*::cos(latMoon)*::sin(2*(lon-lonMoon)))*rx 298 299 + 0.75*Shida*::sin(2*lat) * // north 14b 300 ( sunFactor*::cos(latSun)*::cos(latSun)*::sin(2*(lon-lonSun)) 301 + moonFactor*::cos(latMoon)*::cos(latMoon)*::sin(2*(lon-lonMoon)))*north 302 303 - 1.50*Shida*coslat * // east 14b 304 ( sunFactor*::cos(latSun)*::cos(latSun)*::cos(2*(lon-lonSun)) 305 + moonFactor*::cos(latMoon)*::cos(latMoon)*::cos(2*(lon-lonMoon)))*east; 306 307 LOG(DEBUG6) << "DX2 " << setw(4) << icount << fixed << setprecision(15) 308 << " " << setw(18) << tmp[0] 309 << " " << setw(18) << tmp[1] 310 << " " << setw(18) << tmp[2]; 311 312 for(i=0; i<3; i++) disp[i] += tmp[i]; 313 314 if(debug) { 315 tmp2[0] = northGD[0]*tmp[0] + northGD[1]*tmp[1] + northGD[2]*tmp[2]; 316 tmp2[1] = eastGD[0]*tmp[0] + eastGD[1]*tmp[1] + eastGD[2]*tmp[2]; 317 tmp2[2] = upGD[0]*tmp[0] + upGD[1]*tmp[1] + upGD[2]*tmp[2]; 318 LOG(DEBUG7) << "7SET semi-diurnal-band " << ttag.asGPSString() 319 << fixed << setprecision(9) 320 << " " << tmp[0] << " " << tmp[1] << " " << tmp[2] 321 << " " << tmp2[0] << " " << tmp2[1] << " " << tmp2[2]; 322 } 323 324 // Step 1e IERS(1996) eq. (11) pg 62. latitude dependence of diurnal band 325 Shida = 0.0012; 326 tmp = - 3.0*Shida*sinlat*sinlat * // north 327 ( sunFactor*::cos(latSun)*::sin(latSun)*::cos(lon-lonSun) 328 + moonFactor*::cos(latMoon)*::sin(latMoon)*::cos(lon-lonMoon)) * north 329 330 + 3.0*Shida*sinlat*::cos(2*lat) * // east 331 ( sunFactor*::cos(latSun)*::sin(latSun)*::sin(lon-lonSun) 332 + moonFactor*::cos(latMoon)*::sin(latMoon)*::sin(lon-lonMoon)) * east; 333 334 for(i=0; i<3; i++) { tmp4[i] = tmp[i]; disp[i] += tmp[i]; } 335 336 if(debug) { 337 tmp3 = tmp; // save for below 338 tmp2[0] = northGD[0]*tmp[0] + northGD[1]*tmp[1] + northGD[2]*tmp[2]; 339 tmp2[1] = eastGD[0]*tmp[0] + eastGD[1]*tmp[1] + eastGD[2]*tmp[2]; 340 tmp2[2] = upGD[0]*tmp[0] + upGD[1]*tmp[1] + upGD[2]*tmp[2]; 341 LOG(DEBUG7) << "7SET latitude-diurnal-band " << ttag.asGPSString() 342 << fixed << setprecision(9) 343 << " " << tmp[0] << " " << tmp[1] << " " << tmp[2] 344 << " " << tmp2[0] << " " << tmp2[1] << " " << tmp2[2]; 345 } 346 347 // Step 1f IERS(1996) eq. (12) pg 62. semidiurnal band 348 Shida = 0.0024; 349 tmp = - 1.5*Shida*sinlat*coslat * // north 350 ( sunFactor*::cos(latSun)*::cos(latSun)*::cos(2*(lon-lonSun)) 351 + moonFactor*::cos(latMoon)*::cos(latMoon)*::cos(2*(lon-lonMoon)))*north 352 353 - 1.5*Shida*sinlat*sinlat*coslat * // east 354 ( sunFactor*::cos(latSun)*::cos(latSun)*::sin(2*(lon-lonSun)) 355 + moonFactor*::cos(latMoon)*::cos(latMoon)*::sin(2*(lon-lonMoon)))*east; 356 357 LOG(DEBUG6) << "DX3 " << setw(4) << icount << fixed << setprecision(15) 358 << " " << setw(18) << tmp4[0]+tmp[0] 359 << " " << setw(18) << tmp4[1]+tmp[1] 360 << " " << setw(18) << tmp4[2]+tmp[2]; 361 362 for(i=0; i<3; i++) disp[i] += tmp[i]; 363 364 if(debug) { 365 tmp2[0] = northGD[0]*tmp[0] + northGD[1]*tmp[1] + northGD[2]*tmp[2]; 366 tmp2[1] = eastGD[0]*tmp[0] + eastGD[1]*tmp[1] + eastGD[2]*tmp[2]; 367 tmp2[2] = upGD[0]*tmp[0] + upGD[1]*tmp[1] + upGD[2]*tmp[2]; 368 LOG(DEBUG7) << "7SET latitude-semi-diurnal " << ttag.asGPSString() 369 << fixed << setprecision(9) 370 << " " << tmp[0] << " " << tmp[1] << " " << tmp[2] 371 << " " << tmp2[0] << " " << tmp2[1] << " " << tmp2[2]; 372 373 // combine latitude-dependent terms for comparison with solid.f 374 tmp = tmp + tmp3; 375 tmp2[0] = northGD[0]*tmp[0] + northGD[1]*tmp[1] + northGD[2]*tmp[2]; 376 tmp2[1] = eastGD[0]*tmp[0] + eastGD[1]*tmp[1] + eastGD[2]*tmp[2]; 377 tmp2[2] = upGD[0]*tmp[0] + upGD[1]*tmp[1] + upGD[2]*tmp[2]; 378 LOG(DEBUG7) << "7SET latitude-dependent " << ttag.asGPSString() 379 << fixed << setprecision(9) 380 << " " << tmp[0] << " " << tmp[1] << " " << tmp[2] 381 << " " << tmp2[0] << " " << tmp2[1] << " " << tmp2[2]; 382 } 383 384 // Step 2a IERS(1996) eq. (15) pg 63. 385 // frequency dependence of Love and Shida from diurnal band 386 static double step2diurnalData[9*31] = { 387 -3., 0., 2., 0., 0.,-0.01,-0.01, 0.0, 0.0, 388 -3., 2., 0., 0., 0.,-0.01,-0.01, 0.0, 0.0, 389 -2., 0., 1.,-1., 0.,-0.02,-0.01, 0.0, 0.0, 390 -2., 0., 1., 0., 0.,-0.08, 0.00, 0.01, 0.01, 391 -2., 2.,-1., 0., 0.,-0.02,-0.01, 0.0, 0.0, 392 -1., 0., 0.,-1., 0.,-0.10, 0.00, 0.00, 0.00, 393 -1., 0., 0., 0., 0.,-0.51, 0.00,-0.02, 0.03, 394 -1., 2., 0., 0., 0., 0.01, 0.0, 0.0, 0.0, 395 0.,-2., 1., 0., 0., 0.01, 0.0, 0.0, 0.0, 396 0., 0.,-1., 0., 0., 0.02, 0.01, 0.0, 0.0, 397 0., 0., 1., 0., 0., 0.06, 0.00, 0.00, 0.00, 398 0., 0., 1., 1., 0., 0.01, 0.0, 0.0, 0.0, 399 0., 2.,-1., 0., 0., 0.01, 0.0, 0.0, 0.0, 400 1.,-3., 0., 0., 1.,-0.06, 0.00, 0.00, 0.00, 401 1.,-2., 0., 1., 0., 0.01, 0.0, 0.0, 0.0, 402 1.,-2., 0., 0., 0.,-1.23,-0.07, 0.06, 0.01, 403 1.,-1., 0., 0.,-1., 0.02, 0.0, 0.0, 0.0, 404 1.,-1., 0., 0., 1., 0.04, 0.0, 0.0, 0.0, 405 1., 0., 0.,-1., 0.,-0.22, 0.01, 0.01, 0.00, 406 1., 0., 0., 0., 0.,12.00,-0.78,-0.67,-0.03, 407 1., 0., 0., 1., 0., 1.73,-0.12,-0.10, 0.00, 408 1., 0., 0., 2., 0.,-0.04, 0.0, 0.0, 0.0, 409 1., 1., 0., 0.,-1.,-0.50,-0.01, 0.03, 0.00, 410 1., 1., 0., 0., 1., 0.01, 0.0, 0.0, 0.0, 411 1., 1., 0., 1.,-1.,-0.01, 0.0, 0.0, 0.0, 412 1., 2.,-2., 0., 0.,-0.01, 0.0, 0.0, 0.0, 413 1., 2., 0., 0., 0.,-0.11, 0.01, 0.01, 0.00, 414 2.,-2., 1., 0., 0.,-0.01, 0.0, 0.0, 0.0, 415 2., 0.,-1., 0., 0.,-0.02, 0.02, 0.0, 0.01, 416 3., 0., 0., 0., 0., 0.0, 0.01, 0.0, 0.01, 417 3., 0., 0., 1., 0., 0.0, 0.01, 0.0, 0.0 }; 418 419 // times 420 EphTime TT(ttag); 421 TT.convertSystemTo(TimeSystem::TT); 422 double T,fhr,fmjd = TT.dMJD(); 423 T = (fmjd-51544.0)/36525.0; // MJD of J2000 is 51544.0 424 fhr = (fmjd-int(fmjd))*24.0; 425 426 // compute standard arguments 427 double s,tau,pr,h,p,zns,ps; 428 { 429 double T2 = T*T; 430 double T3 = T2*T; 431 double T4 = T3*T; 432 s = 218.31664563 + 481267.88194*T - 0.0014663889*T2 + 0.00000185139*T3; 433 tau = fhr*15. + 280.4606184 + 36000.7700536*T + 0.00038793*T2 434 - 0.0000000258*T3; 435 tau = tau - s; 436 pr = 1.396971278*T + 0.000308889*T2 + 0.000000021*T3 + 0.000000007*T4; 437 s = s + pr; 438 h = 280.46645 + 36000.7697489*T + 0.00030322222*T2 + 0.000000020*T3 439 - 0.00000000654*T4; 440 p = 83.35324312 + 4069.01363525*T - 0.01032172222*T2 - 0.0000124991*T3 441 + 0.00000005263*T4; 442 zns = 234.95544499 + 1934.13626197*T - 0.00207561111*T2 - 0.00000213944*T3 443 + 0.00000001650*T4; 444 ps = 282.93734098 + 1.71945766667*T + 0.00045688889*T2 - 0.00000001778*T3 445 - 0.00000000334*T4; 446 s = fmod(s, 360.0); 447 tau = fmod(tau,360.0); 448 h = fmod(h, 360.0); 449 p = fmod(p, 360.0); 450 zns = fmod(zns,360.0); 451 ps = fmod(ps, 360.0); 452 } 453 454 double thetaf,ctl,stl,dr,dn,de; 455 tmp = Triple(0,0,0); 456 for(i=0; i<31; i++) { 457 thetaf = (tau + step2diurnalData[0+9*i] * s 458 + step2diurnalData[1+9*i] * h 459 + step2diurnalData[2+9*i] * p 460 + step2diurnalData[3+9*i] * zns 461 + step2diurnalData[4+9*i] * ps) * DEG_TO_RAD; 462 ctl = ::cos(thetaf+lon); 463 stl = ::sin(thetaf+lon); 464 dr = (step2diurnalData[5+9*i] * stl 465 +step2diurnalData[6+9*i] * ctl) * 2 * sinlat * coslat; 466 dn = (step2diurnalData[7+9*i] * stl 467 +step2diurnalData[8+9*i] * ctl) * (coslat*coslat - sinlat*sinlat); 468 de = (step2diurnalData[7+9*i] * ctl 469 -step2diurnalData[8+9*i] * stl) * sinlat; 470 tmp[0] += dr*up[0] + de*east[0] + dn*north[0]; 471 tmp[1] += dr*up[1] + de*east[1] + dn*north[1]; 472 tmp[2] += dr*up[2] + dn*north[2]; 473 } 474 475 // convert total to meters 476 for(i=0; i<3; i++) tmp[i] /= 1000.0; // mm -> m 477 478 LOG(DEBUG6) << "DX4 " << setw(4) << icount << fixed << setprecision(15) 479 << " " << setw(18) << tmp[0] 480 << " " << setw(18) << tmp[1] 481 << " " << setw(18) << tmp[2]; 482 483 for(i=0; i<3; i++) disp[i] += tmp[i]; 484 485 if(debug) { 486 tmp2[0] = northGD[0]*tmp[0] + northGD[1]*tmp[1] + northGD[2]*tmp[2]; 487 tmp2[1] = eastGD[0]*tmp[0] + eastGD[1]*tmp[1] + eastGD[2]*tmp[2]; 488 tmp2[2] = upGD[0]*tmp[0] + upGD[1]*tmp[1] + upGD[2]*tmp[2]; 489 LOG(DEBUG7) << "7SET diurnal-band-corrections " << ttag.asGPSString() 490 << fixed << setprecision(9) 491 << " " << tmp[0] << " " << tmp[1] << " " << tmp[2] 492 << " " << tmp2[0] << " " << tmp2[1] << " " << tmp2[2]; 493 } 494 495 //// Ref. Kouba 2009 and IERS technical note 3, 1989 (out of print) 496 //// Kouba (pers comm 8/12/09) says it is not necessary unless using IERS(1989) 497 //double Tg; // Greenwhich sidereal time 498 //Tg = TT.JD()-2451545.; 499 //if(Tg <= 0.0) Tg -= 1.0; 500 //Tg /= 36525.; 501 //Tg = 0.279057273264 + 100.0021390378009*Tg 502 // + (0.093104 - 6.2e-6*Tg)*Tg*Tg/86400.0; 503 //Tg += TT.secOfDay()/86400.0; // days 504 //Tg = fmod(Tg,1.0); 505 //while(Tg < 0.0) Tg += 1.0; 506 //Tg *= TWO_PI; // radians 507 //dr = -25. * sinlat * coslat * ::sin(Tg + lon); 508 //for(i=0; i<3; i++) tmp[i] -= dr*up[i]; 509 510 // Step 2b IERS(1996) eq. (16) pg 64. 511 // frequency dependence of Love and Shida from the long period band 512 static double step2longData[9*5] = { 513 0, 0, 0, 1, 0, 0.47, 0.23, 0.16, 0.07, 514 0, 2, 0, 0, 0, -0.20,-0.12,-0.11,-0.05, 515 1, 0,-1, 0, 0, -0.11,-0.08,-0.09,-0.04, 516 2, 0, 0, 0, 0, -0.13,-0.11,-0.15,-0.07, 517 2, 0, 0, 1, 0, -0.05,-0.05,-0.06,-0.03 }; 518 519 tmp = Triple(0,0,0); 520 for(i=0; i<5; i++) { 521 thetaf = ( step2longData[0+9*i] * s 522 + step2longData[1+9*i] * h 523 + step2longData[2+9*i] * p 524 + step2longData[3+9*i] * zns 525 + step2longData[4+9*i] * ps) * DEG_TO_RAD; 526 ctl = ::cos(thetaf); 527 stl = ::sin(thetaf); 528 dr = (step2longData[5+9*i] * ctl 529 +step2longData[7+9*i] * stl) * (3*sinlat*sinlat-1)/2; 530 dn = (step2longData[6+9*i] * ctl 531 +step2longData[8+9*i] * stl) * 2*sinlat*coslat; 532 //de = 0.0; remove from next 3 lines 533 tmp[0] += dr*up[0] + dn*north[0]; 534 tmp[1] += dr*up[1] + dn*north[1]; 535 tmp[2] += dr*up[2] + dn*north[2]; 536 } 537 for(i=0; i<3; i++) tmp[i] /= 1000.0; // mm -> m 538 539 LOG(DEBUG6) << "DX5 " << setw(4) << icount << fixed << setprecision(15) 540 << " " << setw(18) << tmp[0] 541 << " " << setw(18) << tmp[1] 542 << " " << setw(18) << tmp[2]; 543 544 for(i=0; i<3; i++) disp[i] += tmp[i]; 545 546 if(debug) { 547 tmp2[0] = northGD[0]*tmp[0] + northGD[1]*tmp[1] + northGD[2]*tmp[2]; 548 tmp2[1] = eastGD[0]*tmp[0] + eastGD[1]*tmp[1] + eastGD[2]*tmp[2]; 549 tmp2[2] = upGD[0]*tmp[0] + upGD[1]*tmp[1] + upGD[2]*tmp[2]; 550 LOG(DEBUG7) << "7SET long-period-band-corr.s " << ttag.asGPSString() 551 << fixed << setprecision(9) 552 << " " << tmp[0] << " " << tmp[1] << " " << tmp[2] 553 << " " << tmp2[0] << " " << tmp2[1] << " " << tmp2[2]; 554 } 555 556 // remove permanent deformation IERS(1996) eq. 17 pg 65. 557 //tmp = -0.1196*(1.5*sinlat*sinlat-0.5)*rx - 0.0247*::sin(2*lat)*north; 558 //for(i=0; i<3; i++) disp[i] += tmp[i]; 559 560 if(debug) { 561 tmp = -0.1196*(1.5*sinlat*sinlat-0.5)*rx - 0.0247*::sin(2*lat)*north; 562 tmp2[0] = northGD[0]*tmp[0] + northGD[1]*tmp[1] + northGD[2]*tmp[2]; 563 tmp2[1] = eastGD[0]*tmp[0] + eastGD[1]*tmp[1] + eastGD[2]*tmp[2]; 564 tmp2[2] = upGD[0]*tmp[0] + upGD[1]*tmp[1] + upGD[2]*tmp[2]; 565 LOG(DEBUG7) << "7SET permanent-tide-not-incl. " << ttag.asGPSString() 566 << fixed << setprecision(9) 567 << " " << tmp[0] << " " << tmp[1] << " " << tmp[2] 568 << " " << tmp2[0] << " " << tmp2[1] << " " << tmp2[2]; 569 } 570 571 if(debug) { 572 for(i=0; i<3; i++) tmp[i] = disp[i]; 573 tmp2[0] = northGD[0]*tmp[0] + northGD[1]*tmp[1] + northGD[2]*tmp[2]; 574 tmp2[1] = eastGD[0]*tmp[0] + eastGD[1]*tmp[1] + eastGD[2]*tmp[2]; 575 tmp2[2] = upGD[0]*tmp[0] + upGD[1]*tmp[1] + upGD[2]*tmp[2]; 576 LOG(DEBUG7) << "7SET total " << ttag.asGPSString() 577 << fixed << setprecision(9) 578 << " " << tmp[0] << " " << tmp[1] << " " << tmp[2] 579 << " " << tmp2[0] << " " << tmp2[1] << " " << tmp2[2]; 580 } 581 582 return disp; 583 } 584 catch(Exception& e) { GPSTK_RETHROW(e); } 585 catch(exception& e) { 586 Exception E("std except: "+string(e.what())); 587 GPSTK_THROW(E); 588 } 589 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); } 590 591 } // end computeSolidEarthTides() 592 593 //--------------------------------------------------------------------------------- 594 /// Compute the site displacement due to rotational deformation due to polar motion 595 /// for the given Position (assumed to fixed to the solid Earth) at the given time, 596 /// given the polar motion angles at time (cf.EarthOrientation). 597 /// Return a Triple containing the site displacement in WGS84 ECEF XYZ coordinates 598 /// with units meters. 599 /// Reference (1996) IERS Technical Note 21 (IERS), ch. 7 page 67. 600 /// Reference (2003) IERS Technical Note 32 (IERS), ch. 7 page 83-84. 601 /// Reference (2010) IERS Technical Note 36 (IERS), ch. 7 page 114-116. 602 /// param site Nominal position of the site of interest. 603 /// param ttag Time of interest. 604 /// param iers IERSConvention IERS convention to use 605 /// param xp double Polar motion angle in arcsec (cf. EarthOrientation) 606 /// param yp double Polar motion angle in arcsec (cf. EarthOrientation) 607 /// return disp Triple disp Displacement vector, ECEF XYZ meters. computePolarTides(const Position site,const EphTime ttag,const double xp,const double yp,const IERSConvention iers)608 Triple computePolarTides(const Position site, const EphTime ttag, 609 const double xp, const double yp, 610 const IERSConvention iers) 611 { 612 try { 613 double m1, m2, upcoef; 614 615 if(iers == IERSConvention::IERS1996) { // 1996 616 m1 = xp; // arcsec 617 m2 = yp; // arcsec 618 upcoef = 0.032; 619 } 620 else { // 2003 and 2010 621 // compute time since J2000 in years and mean pole wander 622 double dt((ttag.dMJD()-51544.5)/365.25); 623 double xmean, ymean; 624 // mean sums in milliarcsec 625 if(iers == IERSConvention::IERS2003) { 626 xmean = (0.054 + 0.00083*dt)/1000.0; // convert to arcsec 627 ymean = (0.357 + 0.00395*dt)/1000.0; // convert to arcsec 628 upcoef = 0.032; 629 } 630 else { // 2003 and 2010 631 // mean sums are different until 2010 and after 2010 (in milliarcsec) 632 if(ttag.year() > 2010) { 633 xmean = 23.513 + 7.6141 * dt; 634 ymean = 358.891 - 0.6287 * dt; 635 } 636 else { 637 xmean = 55.974 + (1.8243 + (0.18413 + 0.007024 * dt) * dt) * dt; 638 ymean = 346.346 + (1.7896 - (0.10729 + 0.000908 * dt) * dt) * dt; 639 } 640 xmean /= 1000.0; // convert to arcsec 641 ymean /= 1000.0; // convert to arcsec 642 // is this 33 a typo in Tech Note 36? other years are 32 643 upcoef = 0.033; 644 } 645 m1 = (xp - xmean); // arcsec 646 m2 = -(yp - ymean); // arcsec 647 } 648 LOG(DEBUG7) << " poletide means " << iers 649 << fixed << setprecision(15) << " " << m1 << " " << m2; 650 651 // the rest is nearly identical in all conventions 652 double lat, lon, theta, sinlat, coslat, sinlon, coslon; 653 Triple disp, dispXYZ; 654 655 lat = site.getGeocentricLatitude(); 656 lon = site.getLongitude(); 657 sinlat = ::sin(lat*DEG_TO_RAD); 658 coslat = ::cos(lat*DEG_TO_RAD); 659 sinlon = ::sin(lon*DEG_TO_RAD); 660 coslon = ::cos(lon*DEG_TO_RAD); 661 theta = (90.0-lat)*DEG_TO_RAD; 662 663 // NEU components (r==Up, theta=S, lambda=E) 664 disp[0] = 0.009 * ::cos(2*theta) * (m1 * coslon + m2 * sinlon); // -S = N 665 disp[1] = 0.009 * ::cos(theta) * (m1 * sinlon - m2 * coslon); // E 666 disp[2] = -upcoef * ::sin(2*theta) * (m1 * coslon + m2 * sinlon); // U 667 668 LOG(DEBUG7) << " poletide " << iers << " (NEU) " 669 << ttag.asGPSString() 670 << fixed << setprecision(9) 671 << disp[0] << " " << disp[1] << " " << disp[2]; 672 673 //// transform X=(x,y,z) into (R*X)(north,east,up) 674 //R(0,0) = -sa*co; R(0,1) = -sa*so; R(0,2) = ca; 675 //R(1,0) = -so; R(1,1) = co; R(1,2) = 0.; 676 //R(2,0) = ca*co; R(2,1) = ca*so; R(2,2) = sa; 677 678 dispXYZ[0] = - sinlat*coslon*disp[0] - sinlon*disp[1] + coslat*coslon*disp[2]; 679 dispXYZ[1] = - sinlat*sinlon*disp[0] + coslon*disp[1] + coslat*sinlon*disp[2]; 680 dispXYZ[2] = coslat*disp[0] + sinlat*disp[2]; 681 682 return dispXYZ; 683 } 684 catch(Exception& e) { GPSTK_RETHROW(e); } 685 } 686 687 } // end namespace gpstk 688 689 //------------------------------------------------------------------------------------ 690 //------------------------------------------------------------------------------------ 691