SolidEarthTides.cpp
Go to the documentation of this file.
1 //==============================================================================
2 //
3 // This file is part of GNSSTk, the ARL:UT GNSS Toolkit.
4 //
5 // The GNSSTk 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 GNSSTk 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 GNSSTk; 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-2022, 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 
50 //------------------------------------------------------------------------------------
51 // GNSSTk
52 // geomatics
53 #include "SolidEarthTides.hpp"
54 #include "logstream.hpp"
55 
56 using namespace std;
57 
58 namespace gnsstk
59 {
60  //---------------------------------------------------------------------------------
61  /* Compute the site displacement due to solid Earth tides for the given
62  Position (assumed to be fixed to the solid Earth) at the given time, given
63  the position of the site of interest, positions and mass ratios of the sun
64  and moon. Return a Triple containing the site displacement in ECEF XYZ
65  coordinates with units meters. Reference IERS Conventions (1996) found in
66  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. */
80  Triple computeSolidEarthTides(const Position& site, const EphTime& ttag,
81  const Position& Sun, const Position& Moon,
82  double EMRAT, double SERAT,
83  const IERSConvention& iers)
84  {
85  try
86  {
87  /* if(ttag.getTimeSystem() == TimeSystem::Unknown) {
88  Exception e("Time system is unknown");
89  GNSSTK_THROW(e);
90  } */
91 
92  // Use REarth from solid.f example program
93  static const double REarth = 6378136.55;
94  static bool debug = (LOGlevel >= DEBUG7);
95  // NB icount is a dummy used in test vs solid.f
96  int i, icount(-1);
97  double RSun, RMoon, Rx, Love, Shida, sunFactor, moonFactor;
98  double sunDOTrx, moonDOTrx, REoRS, REoRM;
99  double lat, lon, sinlat, coslat, sinlon, coslon;
100  double latSun, lonSun, latMoon, lonMoon;
101  Triple disp, sunUnit, moonUnit, rx, tSun, tMoon, north, east, up;
102  // quantities for debug printing only
103  Triple northGD, eastGD, upGD, tmp, tmp2, tmp3, tmp4;
104 
105  LOG(DEBUG7) << "Sun position " << ttag.asGPSString() << fixed
106  << setprecision(3) << setw(23) << Sun.X() << setw(23)
107  << Sun.Y() << setw(23) << Sun.Z();
108  LOG(DEBUG7) << "Moon position" << ttag.asGPSString() << fixed
109  << setprecision(3) << setw(23) << Moon.X() << setw(23)
110  << Moon.Y() << setw(23) << Moon.Z();
111 
112  // distances (m)
113  RSun = Sun.radius();
114  RMoon = Moon.radius();
115  Rx = site.radius();
116 
117  // unit vectors
118  sunUnit = Triple(Sun.X() / RSun, Sun.Y() / RSun, Sun.Z() / RSun);
119  moonUnit =
120  Triple(Moon.X() / RMoon, Moon.Y() / RMoon, Moon.Z() / RMoon);
121  rx = Triple(site.X() / Rx, site.Y() / Rx, site.Z() / Rx);
122 
123  // generate geodetic transformation first - for debug
124  if (debug)
125  {
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
134  // longitude
135  northGD = Triple(-sinlat * coslon, -sinlat * sinlon, coslat);
136  eastGD = Triple(-sinlon, coslon, 0.0);
137  upGD = Triple(coslat * coslon, coslat * sinlon, sinlat);
138  }
139 
140  // use geocentric latitude for formulas
141  latSun = Sun.getGeocentricLatitude() * DEG_TO_RAD;
142  lonSun = Sun.getLongitude() * DEG_TO_RAD;
143  latMoon = Moon.getGeocentricLatitude() * DEG_TO_RAD;
144  lonMoon = Moon.getLongitude() * DEG_TO_RAD;
145  lat = site.getGeocentricLatitude() * DEG_TO_RAD;
146  lon = site.getLongitude() * DEG_TO_RAD;
147  sinlat = ::sin(lat);
148  coslat = ::cos(lat);
149  sinlon = ::sin(lon);
150  coslon = ::cos(lon);
151 
152  /* transform X=(x,y,z) into (R*X)(north,east,up) using geocentric
153  longitude */
154  north = Triple(-sinlat * coslon, -sinlat * sinlon, coslat);
155  east = Triple(-sinlon, coslon, 0.0);
156  up = Triple(coslat * coslon, coslat * sinlon, sinlat);
157 
158  // GM*R factors
159  REoRS = REarth / RSun; // ratio Earth/Sun radius
160  sunFactor =
161  REarth * REoRS * REoRS * REoRS * SERAT; // = (GMS/GME)*RE^4/RS^3
162  REoRM = REarth / RMoon; // ratio Earth/Moon radius
163  moonFactor =
164  REarth * REoRM * REoRM * REoRM / EMRAT; // = (GMM/GME)*RE^4/RM^3
165  /* E/M mass ratio (403) 81.300584999999998 (405) 81.300560000000004
166  S/E mass ratio (403) 332946.048630181234330 (405) 332946.050894783285912
167  LOG(INFO) << " E/M mass ratio " << fixed << setprecision(15) <<
168  EMRAT; LOG(INFO) << " S/E mass ratio " << fixed << setprecision(15)
169  << SERAT; */
170 
171  // dot products
172  sunDOTrx = sunUnit.dot(rx);
173  moonDOTrx = moonUnit.dot(rx);
174 
175  // transverse to radial direction - not unit vectors
176  tSun = sunUnit - sunDOTrx * rx;
177  tMoon = moonUnit - moonDOTrx * rx;
178 
179  // ----------------------------------------------------------------------
180  // compute displacements
181  disp = Triple(0, 0, 0);
182  // Steps and equations refer to IERS(1996), esp table on page 60;
183  // formulas are generally repeated in other IERS technical notes.
184 
185  // Step 1a IERS(1996) eq. (8) pg 61.
186  // nominal degree 2 Love and Shida numbers pg 60
187  double poly = sinlat;
188  poly = (3.0 * poly * poly - 1.0) / 2.0;
189 
190  // here is the only difference between 1996 and 2003/10
191  if (iers == IERSConvention::IERS1996)
192  {
193  Love = 0.6026 - 0.0006 * poly;
194  Shida = 0.0831 + 0.0002 * poly;
195  }
196  else
197  { // 2003 or 2010
198  Love = 0.6078 - 0.0006 * poly;
199  Shida = 0.0847 + 0.0002 * poly;
200  }
201  LOG(DEBUG6) << "H2L2 " << setw(4) << icount << fixed
202  << setprecision(15) << " " << setw(18) << Love << " "
203  << setw(18) << Shida << " " << setw(18) << poly;
204  LOG(DEBUG6) << "P2 " << setw(4) << icount << fixed << setprecision(15)
205  << " " << setw(18)
206  << 3 * (Love / 2 - Shida) * sunDOTrx * sunDOTrx -
207  0.5 * Love
208  << " " << setw(18)
209  << 3 * (Love / 2 - Shida) * moonDOTrx * moonDOTrx -
210  0.5 * Love;
211 
212  tmp = sunFactor * (Love * (1.5 * sunDOTrx * sunDOTrx - 0.5) * rx +
213  3.0 * Shida * sunDOTrx * tSun);
214  for (i = 0; i < 3; i++)
215  disp[i] += tmp[i];
216 
217  tmp = moonFactor * (Love * (1.5 * moonDOTrx * moonDOTrx - 0.5) * rx +
218  3.0 * Shida * moonDOTrx * tMoon);
219  for (i = 0; i < 3; i++)
220  disp[i] += tmp[i];
221 
222  // Step 1b IERS(1996) eq. (9) pg 61.
223  // nominal degree 3 Love and Shida numbers pg 60
224  double Shida2 = Shida;
225  Love = 0.292;
226  Shida = 0.015;
227  tmp = sunFactor * REoRS *
228  (Love * (2.5 * sunDOTrx * sunDOTrx - 1.5) * sunDOTrx * rx +
229  Shida * (7.5 * sunDOTrx * sunDOTrx - 1.5) * tSun);
230  for (i = 0; i < 3; i++)
231  disp[i] += tmp[i];
232 
233  LOG(DEBUG6) << "P3 " << setw(4) << icount << fixed << setprecision(15)
234  << " " << setw(18)
235  << 2.5 * (Love - 3 * Shida) * sunDOTrx * sunDOTrx *
236  sunDOTrx +
237  1.5 * (Shida - Love) * sunDOTrx
238  << " " << setw(18)
239  << 2.5 * (Love - 3 * Shida) * moonDOTrx * moonDOTrx *
240  moonDOTrx +
241  1.5 * (Shida - Love) * moonDOTrx;
242  LOG(DEBUG6) << "X2 " << setw(4) << icount << fixed << setprecision(15)
243  << " " << setw(18) << 3 * Shida2 * sunDOTrx << " "
244  << setw(18) << 3 * Shida2 * moonDOTrx;
245  LOG(DEBUG6) << "X3 " << setw(4) << icount << fixed << setprecision(15)
246  << " " << setw(18)
247  << 1.5 * Shida * (5 * sunDOTrx * sunDOTrx - 1) << " "
248  << setw(18)
249  << 1.5 * Shida * (5 * moonDOTrx * moonDOTrx - 1);
250  LOG(DEBUG6) << "RAT " << setw(4) << icount << fixed << setprecision(6)
251  << " " << setw(18) << SERAT << setprecision(15) << " "
252  << setw(22) << EMRAT << setprecision(2) << " " << setw(11)
253  << REarth;
254  LOG(DEBUG6) << "FACT2 " << setw(4) << icount << fixed
255  << setprecision(15) << " " << setw(18) << sunFactor << " "
256  << setw(18) << moonFactor;
257  LOG(DEBUG6) << "FACT3 " << setw(4) << icount << fixed
258  << setprecision(15) << " " << setw(18) << sunFactor * REoRS
259  << " " << setw(18) << moonFactor * REoRM;
260 
261  tmp = moonFactor * REoRM *
262  (Love * (2.5 * moonDOTrx * moonDOTrx - 1.5) * moonDOTrx * rx +
263  Shida * (7.5 * moonDOTrx * moonDOTrx - 1.5) * tMoon);
264  for (i = 0; i < 3; i++)
265  disp[i] += tmp[i];
266 
267  // all of 8 and 9
268  if (debug)
269  {
270  tmp2[0] = northGD[0] * disp[0] + northGD[1] * disp[1] +
271  northGD[2] * disp[2];
272  tmp2[1] =
273  eastGD[0] * disp[0] + eastGD[1] * disp[1] + eastGD[2] * disp[2];
274  tmp2[2] = upGD[0] * disp[0] + upGD[1] * disp[1] + upGD[2] * disp[2];
275  LOG(DEBUG7) << "7SET solar/lunar/2nd/3rd " << ttag.asGPSString()
276  << fixed << setprecision(9) << " " << disp[0] << " "
277  << disp[1] << " " << disp[2] << " " << tmp2[0] << " "
278  << tmp2[1] << " " << tmp2[2];
279  }
280  LOG(DEBUG6) << "DX0 " << setw(4) << icount << fixed << setprecision(15)
281  << " " << setw(18) << disp[0] << " " << setw(18) << disp[1]
282  << " " << setw(18) << disp[2];
283 
284  // Step 1c IERS(1996) eq. (13) pg 63. diurnal tides
285  Love = -0.0025;
286  Shida = -0.0007;
287  tmp = -0.75 * Love * ::sin(2 * lat) * // radial 13a
288  (sunFactor * ::sin(2 * latSun) * ::sin(lon - lonSun) +
289  moonFactor * ::sin(2 * latMoon) * ::sin(lon - lonMoon)) *
290  rx
291 
292  - 1.5 * Shida * ::cos(2 * lat) * // north 13b
293  (sunFactor * ::sin(2 * latSun) * ::sin(lon - lonSun) +
294  moonFactor * ::sin(2 * latMoon) * ::sin(lon - lonMoon)) *
295  north
296 
297  - 1.5 * Shida * sinlat * // east 13b
298  (sunFactor * ::sin(2 * latSun) * ::cos(lon - lonSun) +
299  moonFactor * ::sin(2 * latMoon) * ::cos(lon - lonMoon)) *
300  east;
301 
302  LOG(DEBUG6) << "DX1 " << setw(4) << icount << fixed << setprecision(15)
303  << " " << setw(18) << tmp[0] << " " << setw(18) << tmp[1]
304  << " " << setw(18) << tmp[2];
305 
306  for (i = 0; i < 3; i++)
307  disp[i] += tmp[i];
308 
309  if (debug)
310  {
311  tmp2[0] =
312  northGD[0] * tmp[0] + northGD[1] * tmp[1] + northGD[2] * tmp[2];
313  tmp2[1] =
314  eastGD[0] * tmp[0] + eastGD[1] * tmp[1] + eastGD[2] * tmp[2];
315  tmp2[2] = upGD[0] * tmp[0] + upGD[1] * tmp[1] + upGD[2] * tmp[2];
316  LOG(DEBUG7) << "7SET diurnal-band " << ttag.asGPSString() << fixed
317  << setprecision(9) << " " << tmp[0] << " " << tmp[1]
318  << " " << tmp[2] << " " << tmp2[0] << " " << tmp2[1]
319  << " " << tmp2[2];
320  }
321 
322  // Step 1d IERS(1996) eq. (14) pg 63. semidiurnal tides
323  Love = -0.0022;
324  Shida = -0.0007;
325  tmp = -0.75 * Love * coslat * coslat * // radial 14a
326  (sunFactor * ::cos(latSun) * ::cos(latSun) *
327  ::sin(2 * (lon - lonSun)) +
328  moonFactor * ::cos(latMoon) * ::cos(latMoon) *
329  ::sin(2 * (lon - lonMoon))) *
330  rx
331 
332  + 0.75 * Shida * ::sin(2 * lat) * // north 14b
333  (sunFactor * ::cos(latSun) * ::cos(latSun) *
334  ::sin(2 * (lon - lonSun)) +
335  moonFactor * ::cos(latMoon) * ::cos(latMoon) *
336  ::sin(2 * (lon - lonMoon))) *
337  north
338 
339  - 1.50 * Shida * coslat * // east 14b
340  (sunFactor * ::cos(latSun) * ::cos(latSun) *
341  ::cos(2 * (lon - lonSun)) +
342  moonFactor * ::cos(latMoon) * ::cos(latMoon) *
343  ::cos(2 * (lon - lonMoon))) *
344  east;
345 
346  LOG(DEBUG6) << "DX2 " << setw(4) << icount << fixed << setprecision(15)
347  << " " << setw(18) << tmp[0] << " " << setw(18) << tmp[1]
348  << " " << setw(18) << tmp[2];
349 
350  for (i = 0; i < 3; i++)
351  disp[i] += tmp[i];
352 
353  if (debug)
354  {
355  tmp2[0] =
356  northGD[0] * tmp[0] + northGD[1] * tmp[1] + northGD[2] * tmp[2];
357  tmp2[1] =
358  eastGD[0] * tmp[0] + eastGD[1] * tmp[1] + eastGD[2] * tmp[2];
359  tmp2[2] = upGD[0] * tmp[0] + upGD[1] * tmp[1] + upGD[2] * tmp[2];
360  LOG(DEBUG7) << "7SET semi-diurnal-band " << ttag.asGPSString()
361  << fixed << setprecision(9) << " " << tmp[0] << " "
362  << tmp[1] << " " << tmp[2] << " " << tmp2[0] << " "
363  << tmp2[1] << " " << tmp2[2];
364  }
365 
366  // Step 1e IERS(1996) eq. (11) pg 62. latitude dependence of diurnal
367  // band
368  Shida = 0.0012;
369  tmp = -3.0 * Shida * sinlat * sinlat * // north
370  (sunFactor * ::cos(latSun) * ::sin(latSun) *
371  ::cos(lon - lonSun) +
372  moonFactor * ::cos(latMoon) * ::sin(latMoon) *
373  ::cos(lon - lonMoon)) *
374  north
375 
376  + 3.0 * Shida * sinlat * ::cos(2 * lat) * // east
377  (sunFactor * ::cos(latSun) * ::sin(latSun) *
378  ::sin(lon - lonSun) +
379  moonFactor * ::cos(latMoon) * ::sin(latMoon) *
380  ::sin(lon - lonMoon)) *
381  east;
382 
383  for (i = 0; i < 3; i++)
384  {
385  tmp4[i] = tmp[i];
386  disp[i] += tmp[i];
387  }
388 
389  if (debug)
390  {
391  tmp3 = tmp; // save for below
392  tmp2[0] =
393  northGD[0] * tmp[0] + northGD[1] * tmp[1] + northGD[2] * tmp[2];
394  tmp2[1] =
395  eastGD[0] * tmp[0] + eastGD[1] * tmp[1] + eastGD[2] * tmp[2];
396  tmp2[2] = upGD[0] * tmp[0] + upGD[1] * tmp[1] + upGD[2] * tmp[2];
397  LOG(DEBUG7) << "7SET latitude-diurnal-band " << ttag.asGPSString()
398  << fixed << setprecision(9) << " " << tmp[0] << " "
399  << tmp[1] << " " << tmp[2] << " " << tmp2[0] << " "
400  << tmp2[1] << " " << tmp2[2];
401  }
402 
403  // Step 1f IERS(1996) eq. (12) pg 62. semidiurnal band
404  Shida = 0.0024;
405  tmp = -1.5 * Shida * sinlat * coslat * // north
406  (sunFactor * ::cos(latSun) * ::cos(latSun) *
407  ::cos(2 * (lon - lonSun)) +
408  moonFactor * ::cos(latMoon) * ::cos(latMoon) *
409  ::cos(2 * (lon - lonMoon))) *
410  north
411 
412  - 1.5 * Shida * sinlat * sinlat * coslat * // east
413  (sunFactor * ::cos(latSun) * ::cos(latSun) *
414  ::sin(2 * (lon - lonSun)) +
415  moonFactor * ::cos(latMoon) * ::cos(latMoon) *
416  ::sin(2 * (lon - lonMoon))) *
417  east;
418 
419  LOG(DEBUG6) << "DX3 " << setw(4) << icount << fixed << setprecision(15)
420  << " " << setw(18) << tmp4[0] + tmp[0] << " " << setw(18)
421  << tmp4[1] + tmp[1] << " " << setw(18) << tmp4[2] + tmp[2];
422 
423  for (i = 0; i < 3; i++)
424  disp[i] += tmp[i];
425 
426  if (debug)
427  {
428  tmp2[0] =
429  northGD[0] * tmp[0] + northGD[1] * tmp[1] + northGD[2] * tmp[2];
430  tmp2[1] =
431  eastGD[0] * tmp[0] + eastGD[1] * tmp[1] + eastGD[2] * tmp[2];
432  tmp2[2] = upGD[0] * tmp[0] + upGD[1] * tmp[1] + upGD[2] * tmp[2];
433  LOG(DEBUG7) << "7SET latitude-semi-diurnal " << ttag.asGPSString()
434  << fixed << setprecision(9) << " " << tmp[0] << " "
435  << tmp[1] << " " << tmp[2] << " " << tmp2[0] << " "
436  << tmp2[1] << " " << tmp2[2];
437 
438  // combine latitude-dependent terms for comparison with solid.f
439  tmp = tmp + tmp3;
440  tmp2[0] =
441  northGD[0] * tmp[0] + northGD[1] * tmp[1] + northGD[2] * tmp[2];
442  tmp2[1] =
443  eastGD[0] * tmp[0] + eastGD[1] * tmp[1] + eastGD[2] * tmp[2];
444  tmp2[2] = upGD[0] * tmp[0] + upGD[1] * tmp[1] + upGD[2] * tmp[2];
445  LOG(DEBUG7) << "7SET latitude-dependent " << ttag.asGPSString()
446  << fixed << setprecision(9) << " " << tmp[0] << " "
447  << tmp[1] << " " << tmp[2] << " " << tmp2[0] << " "
448  << tmp2[1] << " " << tmp2[2];
449  }
450 
451  // Step 2a IERS(1996) eq. (15) pg 63.
452  // frequency dependence of Love and Shida from diurnal band
453  static double step2diurnalData[9 * 31] = {
454  -3., 0., 2., 0., 0., -0.01, -0.01, 0.0, 0.0,
455  -3., 2., 0., 0., 0., -0.01, -0.01, 0.0, 0.0,
456  -2., 0., 1., -1., 0., -0.02, -0.01, 0.0, 0.0,
457  -2., 0., 1., 0., 0., -0.08, 0.00, 0.01, 0.01,
458  -2., 2., -1., 0., 0., -0.02, -0.01, 0.0, 0.0,
459  -1., 0., 0., -1., 0., -0.10, 0.00, 0.00, 0.00,
460  -1., 0., 0., 0., 0., -0.51, 0.00, -0.02, 0.03,
461  -1., 2., 0., 0., 0., 0.01, 0.0, 0.0, 0.0,
462  0., -2., 1., 0., 0., 0.01, 0.0, 0.0, 0.0,
463  0., 0., -1., 0., 0., 0.02, 0.01, 0.0, 0.0,
464  0., 0., 1., 0., 0., 0.06, 0.00, 0.00, 0.00,
465  0., 0., 1., 1., 0., 0.01, 0.0, 0.0, 0.0,
466  0., 2., -1., 0., 0., 0.01, 0.0, 0.0, 0.0,
467  1., -3., 0., 0., 1., -0.06, 0.00, 0.00, 0.00,
468  1., -2., 0., 1., 0., 0.01, 0.0, 0.0, 0.0,
469  1., -2., 0., 0., 0., -1.23, -0.07, 0.06, 0.01,
470  1., -1., 0., 0., -1., 0.02, 0.0, 0.0, 0.0,
471  1., -1., 0., 0., 1., 0.04, 0.0, 0.0, 0.0,
472  1., 0., 0., -1., 0., -0.22, 0.01, 0.01, 0.00,
473  1., 0., 0., 0., 0., 12.00, -0.78, -0.67, -0.03,
474  1., 0., 0., 1., 0., 1.73, -0.12, -0.10, 0.00,
475  1., 0., 0., 2., 0., -0.04, 0.0, 0.0, 0.0,
476  1., 1., 0., 0., -1., -0.50, -0.01, 0.03, 0.00,
477  1., 1., 0., 0., 1., 0.01, 0.0, 0.0, 0.0,
478  1., 1., 0., 1., -1., -0.01, 0.0, 0.0, 0.0,
479  1., 2., -2., 0., 0., -0.01, 0.0, 0.0, 0.0,
480  1., 2., 0., 0., 0., -0.11, 0.01, 0.01, 0.00,
481  2., -2., 1., 0., 0., -0.01, 0.0, 0.0, 0.0,
482  2., 0., -1., 0., 0., -0.02, 0.02, 0.0, 0.01,
483  3., 0., 0., 0., 0., 0.0, 0.01, 0.0, 0.01,
484  3., 0., 0., 1., 0., 0.0, 0.01, 0.0, 0.0};
485 
486  // times
487  EphTime TT(ttag);
488  TT.convertSystemTo(TimeSystem::TT);
489  double T, fhr, fmjd = TT.dMJD();
490  T = (fmjd - 51544.0) / 36525.0; // MJD of J2000 is 51544.0
491  fhr = (fmjd - int(fmjd)) * 24.0;
492 
493  // compute standard arguments
494  double s, tau, pr, h, p, zns, ps;
495  {
496  double T2 = T * T;
497  double T3 = T2 * T;
498  double T4 = T3 * T;
499  s = 218.31664563 + 481267.88194 * T - 0.0014663889 * T2 +
500  0.00000185139 * T3;
501  tau = fhr * 15. + 280.4606184 + 36000.7700536 * T +
502  0.00038793 * T2 - 0.0000000258 * T3;
503  tau = tau - s;
504  pr = 1.396971278 * T + 0.000308889 * T2 + 0.000000021 * T3 +
505  0.000000007 * T4;
506  s = s + pr;
507  h = 280.46645 + 36000.7697489 * T + 0.00030322222 * T2 +
508  0.000000020 * T3 - 0.00000000654 * T4;
509  p = 83.35324312 + 4069.01363525 * T - 0.01032172222 * T2 -
510  0.0000124991 * T3 + 0.00000005263 * T4;
511  zns = 234.95544499 + 1934.13626197 * T - 0.00207561111 * T2 -
512  0.00000213944 * T3 + 0.00000001650 * T4;
513  ps = 282.93734098 + 1.71945766667 * T + 0.00045688889 * T2 -
514  0.00000001778 * T3 - 0.00000000334 * T4;
515  s = fmod(s, 360.0);
516  tau = fmod(tau, 360.0);
517  h = fmod(h, 360.0);
518  p = fmod(p, 360.0);
519  zns = fmod(zns, 360.0);
520  ps = fmod(ps, 360.0);
521  }
522 
523  double thetaf, ctl, stl, dr, dn, de;
524  tmp = Triple(0, 0, 0);
525  for (i = 0; i < 31; i++)
526  {
527  thetaf = (tau + step2diurnalData[0 + 9 * i] * s +
528  step2diurnalData[1 + 9 * i] * h +
529  step2diurnalData[2 + 9 * i] * p +
530  step2diurnalData[3 + 9 * i] * zns +
531  step2diurnalData[4 + 9 * i] * ps) *
532  DEG_TO_RAD;
533  ctl = ::cos(thetaf + lon);
534  stl = ::sin(thetaf + lon);
535  dr = (step2diurnalData[5 + 9 * i] * stl +
536  step2diurnalData[6 + 9 * i] * ctl) *
537  2 * sinlat * coslat;
538  dn = (step2diurnalData[7 + 9 * i] * stl +
539  step2diurnalData[8 + 9 * i] * ctl) *
540  (coslat * coslat - sinlat * sinlat);
541  de = (step2diurnalData[7 + 9 * i] * ctl -
542  step2diurnalData[8 + 9 * i] * stl) *
543  sinlat;
544  tmp[0] += dr * up[0] + de * east[0] + dn * north[0];
545  tmp[1] += dr * up[1] + de * east[1] + dn * north[1];
546  tmp[2] += dr * up[2] + dn * north[2];
547  }
548 
549  // convert total to meters
550  for (i = 0; i < 3; i++) {
551  tmp[i] /= 1000.0; // mm -> m
552  }
553 
554  LOG(DEBUG6) << "DX4 " << setw(4) << icount << fixed << setprecision(15)
555  << " " << setw(18) << tmp[0] << " " << setw(18) << tmp[1]
556  << " " << setw(18) << tmp[2];
557 
558  for (i = 0; i < 3; i++)
559  disp[i] += tmp[i];
560 
561  if (debug)
562  {
563  tmp2[0] =
564  northGD[0] * tmp[0] + northGD[1] * tmp[1] + northGD[2] * tmp[2];
565  tmp2[1] =
566  eastGD[0] * tmp[0] + eastGD[1] * tmp[1] + eastGD[2] * tmp[2];
567  tmp2[2] = upGD[0] * tmp[0] + upGD[1] * tmp[1] + upGD[2] * tmp[2];
568  LOG(DEBUG7) << "7SET diurnal-band-corrections "
569  << ttag.asGPSString() << fixed << setprecision(9) << " "
570  << tmp[0] << " " << tmp[1] << " " << tmp[2] << " "
571  << tmp2[0] << " " << tmp2[1] << " " << tmp2[2];
572  }
573 
574  /* Ref. Kouba 2009 and IERS technical note 3, 1989 (out of print)
575  Kouba (pers comm 8/12/09) says it is not necessary unless using
576  IERS(1989)
577  double Tg; Greenwhich sidereal time
578  Tg = TT.JD()-2451545.;
579  if(Tg <= 0.0) Tg -= 1.0;
580  Tg /= 36525.;
581  Tg = 0.279057273264 + 100.0021390378009*Tg
582  + (0.093104 - 6.2e-6*Tg)*Tg*Tg/86400.0;
583  Tg += TT.secOfDay()/86400.0; days
584  Tg = fmod(Tg,1.0);
585  while(Tg < 0.0) Tg += 1.0;
586  Tg *= TWO_PI; radians
587  dr = -25. * sinlat * coslat * ::sin(Tg + lon);
588  for(i=0; i<3; i++) tmp[i] -= dr*up[i];
589 
590  Step 2b IERS(1996) eq. (16) pg 64.
591  frequency dependence of Love and Shida from the long period band */
592  static double step2longData[9 * 5] = {
593  0, 0, 0, 1, 0, 0.47, 0.23, 0.16, 0.07,
594  0, 2, 0, 0, 0, -0.20, -0.12, -0.11, -0.05,
595  1, 0, -1, 0, 0, -0.11, -0.08, -0.09, -0.04,
596  2, 0, 0, 0, 0, -0.13, -0.11, -0.15, -0.07,
597  2, 0, 0, 1, 0, -0.05, -0.05, -0.06, -0.03};
598 
599  tmp = Triple(0, 0, 0);
600  for (i = 0; i < 5; i++)
601  {
602  thetaf =
603  (step2longData[0 + 9 * i] * s + step2longData[1 + 9 * i] * h +
604  step2longData[2 + 9 * i] * p + step2longData[3 + 9 * i] * zns +
605  step2longData[4 + 9 * i] * ps) *
606  DEG_TO_RAD;
607  ctl = ::cos(thetaf);
608  stl = ::sin(thetaf);
609  dr = (step2longData[5 + 9 * i] * ctl +
610  step2longData[7 + 9 * i] * stl) *
611  (3 * sinlat * sinlat - 1) / 2;
612  dn = (step2longData[6 + 9 * i] * ctl +
613  step2longData[8 + 9 * i] * stl) *
614  2 * sinlat * coslat;
615  // de = 0.0; remove from next 3 lines
616  tmp[0] += dr * up[0] + dn * north[0];
617  tmp[1] += dr * up[1] + dn * north[1];
618  tmp[2] += dr * up[2] + dn * north[2];
619  }
620  for (i = 0; i < 3; i++)
621  tmp[i] /= 1000.0; // mm -> m
622 
623  LOG(DEBUG6) << "DX5 " << setw(4) << icount << fixed << setprecision(15)
624  << " " << setw(18) << tmp[0] << " " << setw(18) << tmp[1]
625  << " " << setw(18) << tmp[2];
626 
627  for (i = 0; i < 3; i++)
628  disp[i] += tmp[i];
629 
630  if (debug)
631  {
632  tmp2[0] =
633  northGD[0] * tmp[0] + northGD[1] * tmp[1] + northGD[2] * tmp[2];
634  tmp2[1] =
635  eastGD[0] * tmp[0] + eastGD[1] * tmp[1] + eastGD[2] * tmp[2];
636  tmp2[2] = upGD[0] * tmp[0] + upGD[1] * tmp[1] + upGD[2] * tmp[2];
637  LOG(DEBUG7) << "7SET long-period-band-corr.s " << ttag.asGPSString()
638  << fixed << setprecision(9) << " " << tmp[0] << " "
639  << tmp[1] << " " << tmp[2] << " " << tmp2[0] << " "
640  << tmp2[1] << " " << tmp2[2];
641  }
642 
643  /* remove permanent deformation IERS(1996) eq. 17 pg 65.
644  tmp = -0.1196*(1.5*sinlat*sinlat-0.5)*rx -
645  0.0247*::sin(2*lat)*north; for(i=0; i<3; i++) disp[i] += tmp[i]; */
646 
647  if (debug)
648  {
649  tmp = -0.1196 * (1.5 * sinlat * sinlat - 0.5) * rx -
650  0.0247 * ::sin(2 * lat) * north;
651  tmp2[0] =
652  northGD[0] * tmp[0] + northGD[1] * tmp[1] + northGD[2] * tmp[2];
653  tmp2[1] =
654  eastGD[0] * tmp[0] + eastGD[1] * tmp[1] + eastGD[2] * tmp[2];
655  tmp2[2] = upGD[0] * tmp[0] + upGD[1] * tmp[1] + upGD[2] * tmp[2];
656  LOG(DEBUG7) << "7SET permanent-tide-not-incl. "
657  << ttag.asGPSString() << fixed << setprecision(9) << " "
658  << tmp[0] << " " << tmp[1] << " " << tmp[2] << " "
659  << tmp2[0] << " " << tmp2[1] << " " << tmp2[2];
660  }
661 
662  if (debug)
663  {
664  for (i = 0; i < 3; i++)
665  tmp[i] = disp[i];
666  tmp2[0] =
667  northGD[0] * tmp[0] + northGD[1] * tmp[1] + northGD[2] * tmp[2];
668  tmp2[1] =
669  eastGD[0] * tmp[0] + eastGD[1] * tmp[1] + eastGD[2] * tmp[2];
670  tmp2[2] = upGD[0] * tmp[0] + upGD[1] * tmp[1] + upGD[2] * tmp[2];
671  LOG(DEBUG7) << "7SET total " << ttag.asGPSString() << fixed
672  << setprecision(9) << " " << tmp[0] << " " << tmp[1]
673  << " " << tmp[2] << " " << tmp2[0] << " " << tmp2[1]
674  << " " << tmp2[2];
675  }
676 
677  return disp;
678  }
679  catch (Exception& e)
680  {
681  GNSSTK_RETHROW(e);
682  }
683  catch (exception& e)
684  {
685  Exception E("std except: " + string(e.what()));
686  GNSSTK_THROW(E);
687  }
688  catch (...)
689  {
690  Exception e("Unknown exception");
691  GNSSTK_THROW(e);
692  }
693 
694  } // end computeSolidEarthTides()
695 
696  //---------------------------------------------------------------------------------
697  /* Compute the site displacement due to rotational deformation due to polar
698  motion for the given Position (assumed to fixed to the solid Earth) at
699  the given time, given the polar motion angles at time
700  (cf.EarthOrientation). Return a Triple containing the site displacement
701  in WGS84 ECEF XYZ coordinates with units meters.
702  Reference (1996) IERS Technical Note 21 (IERS), ch. 7 page 67.
703  Reference (2003) IERS Technical Note 32 (IERS), ch. 7 page 83-84.
704  Reference (2010) IERS Technical Note 36 (IERS), ch. 7 page 114-116.
705  param site Nominal position of the site of interest.
706  param ttag Time of interest.
707  param iers IERSConvention IERS convention to use
708  param xp double Polar motion angle in arcsec (cf. EarthOrientation)
709  param yp double Polar motion angle in arcsec (cf. EarthOrientation)
710  return disp Triple disp Displacement vector, ECEF XYZ meters. */
711  Triple computePolarTides(const Position& site, const EphTime& ttag,
712  double xp, double yp,
713  const IERSConvention& iers)
714  {
715  try
716  {
717  double m1, m2, upcoef;
718 
719  if (iers == IERSConvention::IERS1996)
720  { // 1996
721  m1 = xp; // arcsec
722  m2 = yp; // arcsec
723  upcoef = 0.032;
724  }
725  else
726  { // 2003 and 2010
727  // compute time since J2000 in years and mean pole wander
728  double dt((ttag.dMJD() - 51544.5) / 365.25);
729  double xmean, ymean;
730  // mean sums in milliarcsec
731  if (iers == IERSConvention::IERS2003)
732  {
733  xmean = (0.054 + 0.00083 * dt) / 1000.0; // convert to arcsec
734  ymean = (0.357 + 0.00395 * dt) / 1000.0; // convert to arcsec
735  upcoef = 0.032;
736  }
737  else
738  { // 2003 and 2010
739  // mean sums are different until 2010 and after 2010 (in
740  // milliarcsec)
741  if (ttag.year() > 2010)
742  {
743  xmean = 23.513 + 7.6141 * dt;
744  ymean = 358.891 - 0.6287 * dt;
745  }
746  else
747  {
748  xmean =
749  55.974 + (1.8243 + (0.18413 + 0.007024 * dt) * dt) * dt;
750  ymean =
751  346.346 + (1.7896 - (0.10729 + 0.000908 * dt) * dt) * dt;
752  }
753  xmean /= 1000.0; // convert to arcsec
754  ymean /= 1000.0; // convert to arcsec
755  // is this 33 a typo in Tech Note 36? other years are 32
756  upcoef = 0.033;
757  }
758  m1 = (xp - xmean); // arcsec
759  m2 = -(yp - ymean); // arcsec
760  }
761  LOG(DEBUG7) << " poletide means " << iers << fixed << setprecision(15)
762  << " " << m1 << " " << m2;
763 
764  // the rest is nearly identical in all conventions
765  double lat, lon, theta, sinlat, coslat, sinlon, coslon;
766  Triple disp, dispXYZ;
767 
768  lat = site.getGeocentricLatitude();
769  lon = site.getLongitude();
770  sinlat = ::sin(lat * DEG_TO_RAD);
771  coslat = ::cos(lat * DEG_TO_RAD);
772  sinlon = ::sin(lon * DEG_TO_RAD);
773  coslon = ::cos(lon * DEG_TO_RAD);
774  theta = (90.0 - lat) * DEG_TO_RAD;
775 
776  // NEU components (r==Up, theta=S, lambda=E)
777  disp[0] =
778  0.009 * ::cos(2 * theta) * (m1 * coslon + m2 * sinlon); // -S = N
779  disp[1] = 0.009 * ::cos(theta) * (m1 * sinlon - m2 * coslon); // E
780  disp[2] =
781  -upcoef * ::sin(2 * theta) * (m1 * coslon + m2 * sinlon); // U
782 
783  LOG(DEBUG7) << " poletide " << iers << " (NEU) " << ttag.asGPSString()
784  << fixed << setprecision(9) << disp[0] << " " << disp[1]
785  << " " << disp[2];
786 
787  /* transform X=(x,y,z) into (R*X)(north,east,up)
788  R(0,0) = -sa*co; R(0,1) = -sa*so; R(0,2) = ca;
789  R(1,0) = -so; R(1,1) = co; R(1,2) = 0.;
790  R(2,0) = ca*co; R(2,1) = ca*so; R(2,2) = sa; */
791 
792  dispXYZ[0] = -sinlat * coslon * disp[0] - sinlon * disp[1] +
793  coslat * coslon * disp[2];
794  dispXYZ[1] = -sinlat * sinlon * disp[0] + coslon * disp[1] +
795  coslat * sinlon * disp[2];
796  dispXYZ[2] = coslat * disp[0] + sinlat * disp[2];
797 
798  return dispXYZ;
799  }
800  catch (Exception& e)
801  {
802  GNSSTK_RETHROW(e);
803  }
804  }
805 
806 } // end namespace gnsstk
807 
808 //------------------------------------------------------------------------------------
809 //------------------------------------------------------------------------------------
gnsstk::TimeSystem::TT
@ TT
Terrestrial time (used in IERS conventions)
SolidEarthTides.hpp
gnsstk::EphTime::asGPSString
std::string asGPSString(const int prec=2) const
Definition: EphTime.hpp:252
gnsstk::computePolarTides
Triple computePolarTides(const Position &site, const EphTime &ttag, double xp, double yp, const IERSConvention &iers)
Definition: SolidEarthTides.cpp:711
logstream.hpp
gnsstk::Position::Y
double Y() const noexcept
return Y coordinate (meters)
Definition: Position.cpp:364
gnsstk::Position::Z
double Z() const noexcept
return Z coordinate (meters)
Definition: Position.cpp:375
gnsstk::Triple
Definition: Triple.hpp:68
gnsstk
For Sinex::InputHistory.
Definition: BasicFramework.cpp:50
LOGlevel
#define LOGlevel
Definition: logstream.hpp:323
std::sin
double sin(gnsstk::Angle x)
Definition: Angle.hpp:144
gnsstk::Exception
Definition: Exception.hpp:151
gnsstk::Triple::dot
double dot(const Triple &right) const noexcept
Definition: Triple.cpp:107
gnsstk::Position::radius
double radius() const noexcept
Definition: Position.cpp:444
gnsstk::Position::getGeocentricLatitude
double getGeocentricLatitude() const noexcept
return geocentric latitude (deg N)
Definition: Position.hpp:459
debug
#define debug
Definition: Rinex3ClockHeader.cpp:51
gnsstk::Position::getGeodeticLatitude
double getGeodeticLatitude() const noexcept
return geodetic latitude (deg N)
Definition: Position.hpp:454
gnsstk::EphTime
Definition: EphTime.hpp:67
std::cos
double cos(gnsstk::Angle x)
Definition: Angle.hpp:146
gnsstk::IERSConvention
IERSConvention
Definition: IERSConvention.hpp:69
GNSSTK_RETHROW
#define GNSSTK_RETHROW(exc)
Definition: Exception.hpp:369
gnsstk::EphTime::year
int year() const
Definition: EphTime.hpp:180
gnsstk::Position::X
double X() const noexcept
return X coordinate (meters)
Definition: Position.cpp:353
LOG
#define LOG(level)
define the macro that is used to write to the log stream
Definition: logstream.hpp:315
std
Definition: Angle.hpp:142
gnsstk::Position::getLongitude
double getLongitude() const noexcept
return longitude (deg E) (either geocentric or geodetic)
Definition: Position.hpp:464
gnsstk::Position
Definition: Position.hpp:136
GNSSTK_THROW
#define GNSSTK_THROW(exc)
Definition: Exception.hpp:366
gnsstk::computeSolidEarthTides
Triple computeSolidEarthTides(const Position &site, const EphTime &ttag, const Position &Sun, const Position &Moon, double EMRAT, double SERAT, const IERSConvention &iers)
Definition: SolidEarthTides.cpp:80
gnsstk::EphTime::dMJD
double dMJD() const
Definition: EphTime.hpp:171
gnsstk::DEBUG7
@ DEBUG7
Definition: logstream.hpp:58
gnsstk::DEBUG6
@ DEBUG6
Definition: logstream.hpp:58
gnsstk::DEG_TO_RAD
static const double DEG_TO_RAD
Conversion Factor from degrees to radians (units: degrees^-1)
Definition: GNSSconstants.hpp:76

gnsstk
Author(s):
autogenerated on Wed Oct 25 2023 02:40:41