EarthOrientation.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 
46 //------------------------------------------------------------------------------------
47 // system includes
48 #include <fstream>
49 // GNSSTk
50 #include "MiscMath.hpp"
51 #include "logstream.hpp"
52 
53 #include "EarthOrientation.hpp"
54 
55 //------------------------------------------------------------------------------------
56 using namespace std;
57 
58 namespace gnsstk
59 {
60  //---------------------------------------------------------------------------------
61  // constants
62  //---------------------------------------------------------------------------------
63  /* JulianEpoch is the epoch for coordTransTime
64  const double EarthOrientation::JulianEpoch=2451545.0; // JD */
65  const double EarthOrientation::JulianEpoch = 51544.5; // MJD
66  const int EarthOrientation::intJulianEpoch = 51544;
67 
68  // pi
69  const double EarthOrientation::TWOPI = 6.283185307179586476925287;
70  const double EarthOrientation::PI = (TWOPI / 2.);
71  const double EarthOrientation::HALFPI = (TWOPI / 4.);
72 
73  // convert degrees to radians
74  const double EarthOrientation::DEG_TO_RAD = 0.0174532925199432957692369;
75 
76  // convert radians to degrees
77  const double EarthOrientation::RAD_TO_DEG = 57.29577951308232087679815;
78 
79  // convert arc seconds to radians
80  const double EarthOrientation::ARCSEC_TO_RAD = 4.848136811095359935899141e-6;
81 
82  // arcseconds in 360 degrees
83  const double EarthOrientation::ARCSEC_PER_CIRCLE = 1296000.0;
84 
85  //---------------------------------------------------------------------------------
86  ostream& operator<<(ostream& os, const EarthOrientation& eo)
87  {
88  os << " " << fixed << setw(10) << setprecision(6) << eo.xp << " "
89  << setw(10) << setprecision(6) << eo.yp << " " << setw(11)
90  << setprecision(7) << eo.UT1mUTC << " " << eo.convention;
91  return os;
92  }
93 
94  //---------------------------------------------------------------------------------
95  // EOP interpolation and correction routines - internal use only
96 
97  /* Compute fundamental arguments at coordTransTime T, which are
98  GMST+pi, L, Lp, F, D, Omega. These are valid for interpolating EOPs in
99  both IERS 2003 and 2010 conventions. param T coordTransTime for epoch of
100  interest param args array of 6 fundamental arguments at T, units radians */
101  static void computeFundamentalArgs(double T, double args[6]);
102 
103  /* Corrections to EOPs xp,yp,UT1mUTC for diurnal and semi-diurnal variations
104  due to ocean tides. Based on IERS routine ortho_eop.f from the USNO web
105  site for the 2010 conventions (NB the 2003 version of ortho_eop is
106  nominally different, however in fact the algorithm and constants are the
107  same; the only difference is that the 2003 algorithm (FTN) is single
108  precision). param mjd time (UTC) of interest param dxp correction to xp
109  in arcseconds param dyp correction to yp in arcseconds param dUT
110  correction to UT1mUTC in seconds */
111  static void correctEOPOceanTides(double mjd, double& dxp, double& dyp,
112  double& dUT);
113 
114  /* Corrections to Earth rotation due to zonal tides using USNO IERS2010
115  algorithm. param args array of 6 fundamental arguments at T, units
116  radians
117  NB result of computeFundamentalArgs(); args[0] (=GMST+pi) is not used
118  param dUT correction to UT1mUTC in seconds
119  param dld correction to length of day in seconds/day
120  param dom correction to yp in radians per second */
121  static void correctEarthRotationZonalTides(const double args[6], double& dUT,
122  double& dld, double& dom);
123 
124  /* Corrections to UT1 and length of day (LOD) due to zonal tides using USNO
125  IERS2003 algorithm.
126  param args array of 6 fundamental arguments at T, units radians
127  NB result of computeFundamentalArgs(); args[0] (=GMST+pi) is not used
128  param dUT correction to UT1mUTC in seconds
129  param dld correction to length of day in seconds/day
130  param dom correction to yp in radians per second */
131  static void correctEarthRotationZonalTides2003(const double args[6],
132  double& dUT, double& dld,
133  double& dom);
134 
135  /* Corrections to UT1 and length of day (LOD) due to subdiurnal librations
136  using USNO IERS2010 algorithm.
137  param args array of 6 fundamental arguments at T (coordTransTime), in
138  radians
139  NB result of computeFundamentalArgs()
140  param dUT correction to UT1mUTC in seconds
141  param dld correction to length of day in seconds/day */
142  static void correctEarthRotationLibrations(const double args[6], double& dUT,
143  double& dld);
144 
145  /* Given parallel arrays of length four containing the values from EOPStore
146  for time (int MJD) and EOPs xp, yp, and UT1mUTC, where the time of
147  interest t lies within the values of the time array, interpolate and apply
148  corrections to determine the EOPs at t, using the algorithm prescribed
149  by the given IERS convention.
150  param t EphTime at which to compute EOPs
151  param time vector of length 4 of consecutive MJDs from EOPStore;
152  t must lie within this timespan.
153  param X vector of length 4 of consecutive xp from EOPStore.
154  param Y vector of length 4 of consecutive yp from EOPStore.
155  param dT vector of length 4 of consecutive UT1mUTC from EOPStore.
156  param conv the IERSConvention to be used.
157  return EarthOrientation EOPs at time t. */
158  void EarthOrientation::interpolateEOP(const EphTime& t,
159  const vector<double>& time,
160  const vector<double>& X,
161  const vector<double>& Y,
162  vector<double>& dT,
163  const IERSConvention& in_conv)
164  {
165  int i, j;
166  double dxp, dyp, dUT, dlod, domega;
167  double args[6];
168 
169  // set the convention for this object
170  convention = in_conv;
171 
172  // first get MJD(UTC), for the Lagrange interpolation
173  EphTime ttag(t);
174  ttag.convertSystemTo(TimeSystem::UTC);
175  double mjdUTC(ttag.dMJD());
176 
177  // now convert to TT, for the corrections algorithms
178  ttag.convertSystemTo(TimeSystem::TT);
179  double mjd(ttag.dMJD());
180  double T = (mjd - 51544.5) / 36525.0;
181 
182  // ----------------------------------------------------------------
183  // step 1 : Lagrange interpolation of xp and yp
184  double err;
185  xp = LagrangeInterpolation(time, X, mjdUTC, err); // arcsec
186  yp = LagrangeInterpolation(time, Y, mjdUTC, err); // arcsec
187  // LOG(INFO) << " -> " << fixed << setprecision(10) << mjdUTC
188  // << " " << setprecision(15) << xp << " " << yp;
189 
190  /* 1a. remove long period tides from UT1-UTC data -------------------
191  NB don't need ttag anymore but do need mjd and T. */
192  for (i = 0; i < time.size(); i++)
193  {
194  ttag.setMJD(time[i]);
195  ttag.setTimeSystem(TimeSystem::UTC);
196  ttag.convertSystemTo(TimeSystem::TT);
197  double Ttemp = (ttag.dMJD() - 51544.5) / 36525.0;
198  computeFundamentalArgs(Ttemp, args);
199  if (convention == IERSConvention::IERS2010)
200  {
201  correctEarthRotationZonalTides(args, dUT, dlod, domega);
202  }
203  else
204  {
205  correctEarthRotationZonalTides2003(args, dUT, dlod, domega);
206  }
207  dT[i] -= dUT;
208  // LOG(INFO) << " UT " <<fixed<< setprecision(10) << time[i] << " " <<
209  // dT[i];
210  }
211 
212  // 1b. interpolate UT1-UTC -------------------------------------------
213  UT1mUTC = LagrangeInterpolation(time, dT, mjdUTC, err); // seconds
214  // LOG(INFO) << " -> " << fixed << setprecision(10) << mjdUTC
215  // << " " << setprecision(15) << UT1mUTC;
216 
217  /* ----------------------------------------------------------------
218  step 2 : Compute fundamental arguments for use in corrections */
219  computeFundamentalArgs(T, args);
220  /* LOG(INFO) << fixed
221  << "T(" << setprecision(2) << mjd << ")=" << setprecision(15) << T
222  << "\nGMST(" << setprecision(2) << mjd << ")="<<setprecision(15)<<
223  args[0]
224  << "\nL(" << setprecision(2) << mjd << ")=" << setprecision(15) <<
225  args[1]
226  << "\nLp(" << setprecision(2) << mjd << ")=" << setprecision(15) <<
227  args[2]
228  << "\nF(" << setprecision(2) << mjd << ")=" << setprecision(15) <<
229  args[3]
230  << "\nD(" << setprecision(2) << mjd << ")=" << setprecision(15) <<
231  args[4]
232  << "\nOm(" << setprecision(2) << mjd << ")=" << setprecision(15)<<
233  args[5]; */
234 
235  // ----------------------------------------------------------------
236  // step 3 : Compute corrections and apply to eop
237 
238  /* 3a. restore long period tides to UT1-UTC --------------------------
239  Corrections to Earth rotation due to zonal tides using USNO algorithm
240  differences between 2003 and 2010 are very small, and only in the Zonal
241  tides */
242  if (convention == IERSConvention::IERS2010)
243  {
244  correctEarthRotationZonalTides(args, dUT, dlod, domega);
245  }
246  else
247  {
248  correctEarthRotationZonalTides2003(args, dUT, dlod, domega);
249  }
250 
251  UT1mUTC += dUT;
252  // LOG(INFO) << " Zonal tides UT1 correction " <<fixed<<setprecision(15)
253  // << dUT;
254 
255  /* step 3b. --------------------------------------
256  Corrections to EOP due to diurnal and semidiurnal effects of ocean
257  tides using USNO algorithm */
258  correctEOPOceanTides(mjd, dxp, dyp, dUT);
259  xp += dxp;
260  yp += dyp;
261  UT1mUTC += dUT;
262  /* LOG(INFO) << " Tides corrections " << fixed << setprecision(15)
263  << dxp << " " << dyp << " " << dUT;
264 
265  step 3c. --------------------------------------
266  correct for librations - a few microarcseconds but high frequency
267  correctEarthRotationLibrations(args, dUT, dlod);
268  UT1mUTC += dUT;
269  LOG(INFO) << " Librations UT1 correction " << fixed<<setprecision(15)
270  << dUT; */
271  }
272 
273  //---------------------------------------------------------------------------------
274  void computeFundamentalArgs(double T, double args[6])
275  {
276 
277  /* compute arguments and their time derivatives GMST+pi, L, Lp, F, D,
278  Omega units are radians and radians/day NB to compute corrections to
279  length of day (LOD), remove '//LOD' everywhere
280  LOD double dargs[6];
281  LOD static const double AS_TO_RAD_PER_DAY=(ARCSEC_TO_RAD/36525.0);
282 
283  GMST+pi
284  GMST in seconds of day. NB 3155760000.0 = 876600.0*3600.0 */
285  double GMST =
286  ::fmod(67310.54841 + T * (3155760000.0 + 8640184.812866 +
287  T * T * (0.093104 + T * (-6.2e-6))),
288  86400.0); // seconds of day
289  /* convert to arcsec, then radians : 360*3600 arcsec / 24*3600 sec = 15
290  as/sec */
291  GMST *= 15.0 * EarthOrientation::ARCSEC_TO_RAD;
292  args[0] = ::fmod(GMST + EarthOrientation::PI, EarthOrientation::TWOPI);
293  // LOD dargs[0] = 15.0*(876600*3600 + 8640184.812866
294  // LOD + T*(2.*0.093104 - T*(3.*6.2e-6))) * AS_TO_RAD_PER_DAY;
295 
296  // fundamental arguments cf. IERS 2010 conventions TN36 5.7.2 eqn 5.43
297 
298  // L = mean anomaly of moon
299  args[1] = EarthOrientation::L(T);
300  // LOD dargs[1] = (1717915923.2178 + T*(2.*31.8792
301  // LOD + T*(3.*0.051635 + T*(-4.*0.00024470)))) *
302  // AS_TO_RAD_PER_DAY;
303 
304  // Lp = mean anomaly of sun
305  args[2] = EarthOrientation::Lp(T);
306  // LOD dargs[2] = (129596581.0481 + T*(-2.*0.5532
307  // LOD + T*(-3.*0.000136 + T*(-4.*0.00001149)))) * AS_TO_RAD_PER_DAY;
308 
309  // F = mean anomaly of moon minus Omega
310  args[3] = EarthOrientation::F(T);
311  // LOD dargs[3] = (1739527262.8478 + T*(-2.*12.7512
312  // LOD + T*(-3.*0.001037 + T*(4.*0.00000417)))) * AS_TO_RAD_PER_DAY;
313 
314  // D = mean elongation moon from sun
315  args[4] = EarthOrientation::D(T);
316  // LOD dargs[4] = (1602961601.2090 + T*(-2.*6.3706
317  // LOD + T*(+3.*0.006593 + T*(-4.*0.00003169)))) * AS_TO_RAD_PER_DAY;
318 
319  // Omega = mean longitude of ascending node of moon
320  args[5] = EarthOrientation::Omega2003(T);
321  // LOD dargs[5] = (-6962890.2665 + T*(2.*7.4722
322  // LOD + T*(3.*0.007702 + T*(-4.*0.00005939)))) * AS_TO_RAD_PER_DAY;
323 
324  } // end computeFundamentalArgs()
325 
326  //---------------------------------------------------------------------------------
327  void correctEOPOceanTides(double mjd, double& dxp, double& dyp,
328  double& dUT)
329  {
330  /* compute time dependent part of second degree diurnal and semidiurnal
331  tidal potential from the dominant spectral lines of the
332  Cartwright-Tayler-Edden harmonic decomposition. */
333 
334  // define the orthotide weight factors
335  static const double fact[][2] = {{0.0298, 0.0200}, {0.1408, 0.0905},
336  {0.0805, 0.0638}, {0.6002, 0.3476},
337  {0.3025, 0.1645}, {0.1517, 0.0923}};
338  /* for(int kk=0; kk<6; kk++)
339  LOG(INFO) << "SP(" << kk << ")=" << fixed << setprecision(5)
340  << setw(9) << fact[kk][0] << setw(9) << fact[kk][1]; */
341 
342  // tidal potential model for 71 diurnal and semidiurnal lines
343  typedef struct
344  {
345  int nj, mj; // nj is always 2!
346  double hs, phase, freq;
347  } Coeff;
348  static const Coeff C[] = {{2, 1, -1.94, 9.0899831, 5.18688050},
349  {2, 1, -1.25, 8.8234208, 5.38346657},
350  {2, 1, -6.64, 12.1189598, 5.38439079},
351  {2, 1, -1.51, 1.4425700, 5.41398343},
352  {2, 1, -8.02, 4.7381090, 5.41490765},
353  {2, 1, -9.47, 4.4715466, 5.61149372},
354  {2, 1, -50.20, 7.7670857, 5.61241794},
355  {2, 1, -1.80, -2.9093042, 5.64201057},
356  {2, 1, -9.54, 0.3862349, 5.64293479},
357  {2, 1, 1.52, -3.1758666, 5.83859664},
358  {2, 1, -49.45, 0.1196725, 5.83952086},
359  {2, 1, -262.21, 3.4152116, 5.84044508},
360  {2, 1, 1.70, 12.8946194, 5.84433381},
361  {2, 1, 3.43, 5.5137686, 5.87485066},
362  {2, 1, 1.94, 6.4441883, 6.03795537},
363  {2, 1, 1.37, -4.2322016, 6.06754801},
364  {2, 1, 7.41, -0.9366625, 6.06847223},
365  {2, 1, 20.62, 8.5427453, 6.07236095},
366  {2, 1, 4.14, 11.8382843, 6.07328517},
367  {2, 1, 3.94, 1.1618945, 6.10287781},
368  {2, 1, -7.14, 5.9693878, 6.24878055},
369  {2, 1, 1.37, -1.2032249, 6.26505830},
370  {2, 1, -122.03, 2.0923141, 6.26598252},
371  {2, 1, 1.02, -1.7847596, 6.28318449},
372  {2, 1, 2.89, 8.0679449, 6.28318613},
373  {2, 1, -7.30, 0.8953321, 6.29946388},
374  {2, 1, 368.78, 4.1908712, 6.30038810},
375  {2, 1, 50.01, 7.4864102, 6.30131232},
376  {2, 1, -1.08, 10.7819493, 6.30223654},
377  {2, 1, 2.93, 0.3137975, 6.31759007},
378  {2, 1, 5.25, 6.2894282, 6.33479368},
379  {2, 1, 3.95, 7.2198478, 6.49789839},
380  {2, 1, 20.62, -0.1610030, 6.52841524},
381  {2, 1, 4.09, 3.1345361, 6.52933946},
382  {2, 1, 3.42, 2.8679737, 6.72592553},
383  {2, 1, 1.69, -4.5128771, 6.75644239},
384  {2, 1, 11.29, 4.9665307, 6.76033111},
385  {2, 1, 7.23, 8.2620698, 6.76125533},
386  {2, 1, 1.51, 11.5576089, 6.76217955},
387  {2, 1, 2.16, 0.6146566, 6.98835826},
388  {2, 1, 1.38, 3.9101957, 6.98928248},
389  {2, 2, 1.80, 20.6617051, 11.45675174},
390  {2, 2, 4.67, 13.2808543, 11.48726860},
391  {2, 2, 16.01, 16.3098310, 11.68477889},
392  {2, 2, 19.32, 8.9289802, 11.71529575},
393  {2, 2, 1.30, 5.0519065, 11.73249771},
394  {2, 2, -1.02, 15.8350306, 11.89560406},
395  {2, 2, -4.51, 8.6624178, 11.91188181},
396  {2, 2, 120.99, 11.9579569, 11.91280603},
397  {2, 2, 1.13, 8.0808832, 11.93000800},
398  {2, 2, 22.98, 4.5771061, 11.94332289},
399  {2, 2, 1.06, 0.7000324, 11.96052486},
400  {2, 2, -1.90, 14.9869335, 12.11031632},
401  {2, 2, -2.18, 11.4831564, 12.12363121},
402  {2, 2, -23.58, 4.3105437, 12.13990896},
403  {2, 2, 631.92, 7.6060827, 12.14083318},
404  {2, 2, 1.92, 3.7290090, 12.15803515},
405  {2, 2, -4.66, 10.6350594, 12.33834347},
406  {2, 2, -17.86, 3.2542086, 12.36886033},
407  {2, 2, 4.47, 12.7336164, 12.37274905},
408  {2, 2, 1.97, 16.0291555, 12.37367327},
409  {2, 2, 17.20, 10.1602590, 12.54916865},
410  {2, 2, 294.00, 6.2831853, 12.56637061},
411  {2, 2, -2.46, 2.4061116, 12.58357258},
412  {2, 2, -1.02, 5.0862033, 12.59985198},
413  {2, 2, 79.96, 8.3817423, 12.60077620},
414  {2, 2, 23.83, 11.6772814, 12.60170041},
415  {2, 2, 2.59, 14.9728205, 12.60262463},
416  {2, 2, 4.47, 4.0298682, 12.82880334},
417  {2, 2, 1.95, 7.3254073, 12.82972756},
418  {2, 2, 1.17, 9.1574019, 13.06071921}};
419 
420  static const double dt(2.0);
421  static const double TWOPI(6.283185307179586476925287);
422  int j, k, n, m;
423  double a[2][3], b[2][3], dt60, pinm, alpha;
424 
425  // compute time dependent potential matrix
426  for (k = 0; k < 3; ++k)
427  {
428  dt60 = (mjd - (k - 1) * dt) - 37076.5;
429  a[0][k] = a[1][k] = 0.0;
430  b[0][k] = b[1][k] = 0.0;
431  for (j = 0; j < 71; ++j)
432  {
433  n = C[j].nj; // 2!
434  m = C[j].mj; // 1 or 2
435  pinm = double((n + m) % 2) * TWOPI / 4.0;
436  alpha = ::fmod(C[j].phase - pinm, TWOPI) +
437  ::fmod(C[j].freq * dt60, TWOPI);
438  a[m - 1][k] += C[j].hs * ::cos(alpha);
439  b[m - 1][k] -= C[j].hs * ::sin(alpha);
440  }
441  }
442 
443  // orthogonalize response terms
444  double ap, am, bp, bm, p[3][2], q[3][2];
445  for (m = 0; m < 2; ++m)
446  {
447  ap = a[m][2] + a[m][0];
448  am = a[m][2] - a[m][0];
449  bp = b[m][2] + b[m][0];
450  bm = b[m][2] - b[m][0];
451  p[0][m] = fact[0][m] * a[m][1];
452  p[1][m] = fact[1][m] * a[m][1] - fact[2][m] * ap;
453  p[2][m] = fact[3][m] * a[m][1] - fact[4][m] * ap + fact[5][m] * bm;
454  q[0][m] = fact[0][m] * b[m][1];
455  q[1][m] = fact[1][m] * b[m][1] - fact[2][m] * bp;
456  q[2][m] = fact[3][m] * b[m][1] - fact[4][m] * bp - fact[5][m] * am;
457  a[m][0] = p[0][m];
458  a[m][1] = p[1][m];
459  a[m][2] = p[2][m];
460  b[m][0] = q[0][m];
461  b[m][1] = q[1][m];
462  b[m][2] = q[2][m];
463  }
464 
465  // fill partials vector
466  double h[12];
467  for (j = 0, m = 0; m < 2; ++m)
468  {
469  for (k = 0; k < 3; ++k)
470  {
471  h[j] = a[m][k];
472  h[j + 1] = b[m][k];
473  j += 2;
474  }
475  }
476  // for(k=0; k<12; k++)
477  // LOG(INFO) << "H(" << k << ")=" << fixed << setprecision(15) << h[k];
478 
479  // orthoweights
480  static const double orthowts[][3] = {
481  {-6.77832, 14.86283, -1.76335}, {-14.86323, -6.77846, 1.03364},
482  {0.47884, 1.45234, -0.27553}, {-1.45303, 0.47888, 0.34569},
483  {0.16406, -0.42056, -0.12343}, {0.42030, 0.16469, -0.10146},
484  {0.09398, 15.30276, -0.47119}, {25.73054, -4.30615, 1.28997},
485  {-4.77974, 0.07564, -0.19336}, {0.28080, 2.28321, 0.02724},
486  {1.94539, -0.45717, 0.08955}, {-0.73089, -1.62010, 0.04726}};
487  // for(k=0; k<12; k++)
488  // LOG(INFO) << "ORTHOW(" << k << ")=" << fixed << setprecision(5)
489  // << setw(9) << orthowts[k][0]
490  // << setw(9) << orthowts[k][1]
491  // << setw(9) << orthowts[k][2];
492 
493  double eop[3];
494  for (k = 0; k < 3; ++k)
495  {
496  eop[k] = 0.0;
497  for (j = 0; j < 12; ++j)
498  {
499  eop[k] += h[j] * orthowts[j][k];
500  }
501  }
502 
503  // convert to arcsec and seconds
504  dxp = eop[0] * 1.e-6;
505  dyp = eop[1] * 1.e-6;
506  dUT = eop[2] * 1.e-6;
507  }
508 
509  //---------------------------------------------------------------------------------
510  /* Corrections to Earth rotation due to zonal tides using USNO algorithm
511  param args array of 6 fundamental arguments at T, units radians
512  NB args[0] (=GMST+pi) is not used
513  param dUT correction to UT1mUTC in seconds
514  param dld correction to length of day in seconds/day
515  param dom correction to yp in radians per second */
516  void correctEarthRotationZonalTides(const double args[6], double& dUT,
517  double& dld, double& dom)
518  {
519  // constants
520  static const double TWOPI(6.283185307179586476925287);
521 
522  /* luni-solar coefficients in argument multiplying L Lp F D Omega, 62
523  terms and coefficients of DUT sin, cos; DLOD cos, sin; DOMEGA cos, sin */
524  typedef struct
525  {
526  int nargs[5];
527  double dutsin, dutcos, dldcos, dldsin, domcos,
528  domsin; // yes its sc cs cs
529  } Coeff;
530 
531  static const Coeff C[] = {
532  {{1, 0, 2, 2, 2}, -0.0235, 0.0000, 0.2617, 0.0000, -0.2209, 0.0000},
533  {{2, 0, 2, 0, 1}, -0.0404, 0.0000, 0.3706, 0.0000, -0.3128, 0.0000},
534  {{2, 0, 2, 0, 2}, -0.0987, 0.0000, 0.9041, 0.0000, -0.7630, 0.0000},
535  {{0, 0, 2, 2, 1}, -0.0508, 0.0000, 0.4499, 0.0000, -0.3797, 0.0000},
536  {{0, 0, 2, 2, 2}, -0.1231, 0.0000, 1.0904, 0.0000, -0.9203, 0.0000},
537  {{1, 0, 2, 0, 0}, -0.0385, 0.0000, 0.2659, 0.0000, -0.2244, 0.0000},
538  {{1, 0, 2, 0, 1}, -0.4108, 0.0000, 2.8298, 0.0000, -2.3884, 0.0000},
539  {{1, 0, 2, 0, 2}, -0.9926, 0.0000, 6.8291, 0.0000, -5.7637, 0.0000},
540  {{3, 0, 0, 0, 0}, -0.0179, 0.0000, 0.1222, 0.0000, -0.1031, 0.0000},
541  {{-1, 0, 2, 2, 1}, -0.0818, 0.0000, 0.5384, 0.0000, -0.4544, 0.0000},
542  {{-1, 0, 2, 2, 2}, -0.1974, 0.0000, 1.2978, 0.0000, -1.0953, 0.0000},
543  {{1, 0, 0, 2, 0}, -0.0761, 0.0000, 0.4976, 0.0000, -0.4200, 0.0000},
544  {{2, 0, 2, -2, 2}, 0.0216, 0.0000, -0.1060, 0.0000, 0.0895, 0.0000},
545  {{0, 1, 2, 0, 2}, 0.0254, 0.0000, -0.1211, 0.0000, 0.1022, 0.0000},
546  {{0, 0, 2, 0, 0}, -0.2989, 0.0000, 1.3804, 0.0000, -1.1650, 0.0000},
547  {{0, 0, 2, 0, 1}, -3.1873, 0.2010, 14.6890, 0.9266, -12.3974, -0.7820},
548  {{0, 0, 2, 0, 2}, -7.8468, 0.5320, 36.0910, 2.4469, -30.4606, -2.0652},
549  {{2, 0, 0, 0, -1}, 0.0216, 0.0000, -0.0988, 0.0000, 0.0834, 0.0000},
550  {{2, 0, 0, 0, 0}, -0.3384, 0.0000, 1.5433, 0.0000, -1.3025, 0.0000},
551  {{2, 0, 0, 0, 1}, 0.0179, 0.0000, -0.0813, 0.0000, 0.0686, 0.0000},
552  {{0, -1, 2, 0, 2}, -0.0244, 0.0000, 0.1082, 0.0000, -0.0913, 0.0000},
553  {{0, 0, 0, 2, -1}, 0.0470, 0.0000, -0.2004, 0.0000, 0.1692, 0.0000},
554  {{0, 0, 0, 2, 0}, -0.7341, 0.0000, 3.1240, 0.0000, -2.6367, 0.0000},
555  {{0, 0, 0, 2, 1}, -0.0526, 0.0000, 0.2235, 0.0000, -0.1886, 0.0000},
556  {{0, -1, 0, 2, 0}, -0.0508, 0.0000, 0.2073, 0.0000, -0.1749, 0.0000},
557  {{1, 0, 2, -2, 1}, 0.0498, 0.0000, -0.1312, 0.0000, 0.1107, 0.0000},
558  {{1, 0, 2, -2, 2}, 0.1006, 0.0000, -0.2640, 0.0000, 0.2228, 0.0000},
559  {{1, 1, 0, 0, 0}, 0.0395, 0.0000, -0.0968, 0.0000, 0.0817, 0.0000},
560  {{-1, 0, 2, 0, 0}, 0.0470, 0.0000, -0.1099, 0.0000, 0.0927, 0.0000},
561  {{-1, 0, 2, 0, 1}, 0.1767, 0.0000, -0.4115, 0.0000, 0.3473, 0.0000},
562  {{-1, 0, 2, 0, 2}, 0.4352, 0.0000, -1.0093, 0.0000, 0.8519, 0.0000},
563  {{1, 0, 0, 0, -1}, 0.5339, 0.0000, -1.2224, 0.0000, 1.0317, 0.0000},
564  {{1, 0, 0, 0, 0}, -8.4046, 0.2500, 19.1647, 0.5701, -16.1749, -0.4811},
565  {{1, 0, 0, 0, 1}, 0.5443, 0.0000, -1.2360, 0.0000, 1.0432, 0.0000},
566  {{0, 0, 0, 1, 0}, 0.0470, 0.0000, -0.1000, 0.0000, 0.0844, 0.0000},
567  {{1, -1, 0, 0, 0}, -0.0555, 0.0000, 0.1169, 0.0000, -0.0987, 0.0000},
568  {{-1, 0, 0, 2, -1}, 0.1175, 0.0000, -0.2332, 0.0000, 0.1968, 0.0000},
569  {{-1, 0, 0, 2, 0}, -1.8236, 0.0000, 3.6018, 0.0000, -3.0399, 0.0000},
570  {{-1, 0, 0, 2, 1}, 0.1316, 0.0000, -0.2587, 0.0000, 0.2183, 0.0000},
571  {{1, 0, -2, 2, -1}, 0.0179, 0.0000, -0.0344, 0.0000, 0.0290, 0.0000},
572  {{-1, -1, 0, 2, 0}, -0.0855, 0.0000, 0.1542, 0.0000, -0.1302, 0.0000},
573  {{0, 2, 2, -2, 2}, -0.0573, 0.0000, 0.0395, 0.0000, -0.0333, 0.0000},
574  {{0, 1, 2, -2, 1}, 0.0329, 0.0000, -0.0173, 0.0000, 0.0146, 0.0000},
575  {{0, 1, 2, -2, 2}, -1.8847, 0.0000, 0.9726, 0.0000, -0.8209, 0.0000},
576  {{0, 0, 2, -2, 0}, 0.2510, 0.0000, -0.0910, 0.0000, 0.0768, 0.0000},
577  {{0, 0, 2, -2, 1}, 1.1703, 0.0000, -0.4135, 0.0000, 0.3490, 0.0000},
578  {{0, 0, 2, -2, 2}, -49.7174, 0.4330, 17.1056, 0.1490, -14.4370, -0.1257},
579  {{0, 2, 0, 0, 0}, -0.1936, 0.0000, 0.0666, 0.0000, -0.0562, 0.0000},
580  {{2, 0, 0, -2, -1}, 0.0489, 0.0000, -0.0154, 0.0000, 0.0130, 0.0000},
581  {{2, 0, 0, -2, 0}, -0.5471, 0.0000, 0.1670, 0.0000, -0.1409, 0.0000},
582  {{2, 0, 0, -2, 1}, 0.0367, 0.0000, -0.0108, 0.0000, 0.0092, 0.0000},
583  {{0, -1, 2, -2, 1}, -0.0451, 0.0000, 0.0082, 0.0000, -0.0069, 0.0000},
584  {{0, 1, 0, 0, -1}, 0.0921, 0.0000, -0.0167, 0.0000, 0.0141, 0.0000},
585  {{0, -1, 2, -2, 2}, 0.8281, 0.0000, -0.1425, 0.0000, 0.1202, 0.0000},
586  {{0, 1, 0, 0, 0}, -15.8887, 0.1530, 2.7332, 0.0267, -2.3068, -0.0222},
587  {{0, 1, 0, 0, 1}, -0.1382, 0.0000, 0.0225, 0.0000, -0.0190, 0.0000},
588  {{1, 0, 0, -1, 0}, 0.0348, 0.0000, -0.0053, 0.0000, 0.0045, 0.0000},
589  {{2, 0, -2, 0, 0}, -0.1372, 0.0000, -0.0079, 0.0000, 0.0066, 0.0000},
590  {{-2, 0, 2, 0, 1}, 0.4211, 0.0000, -0.0203, 0.0000, 0.0171, 0.0000},
591  {{-1, 1, 0, 1, 0}, -0.0404, 0.0000, 0.0008, 0.0000, -0.0007, 0.0000},
592  {{0, 0, 0, 0, 2}, 7.8998, 0.0000, 0.1460, 0.0000, -0.1232, 0.0000},
593  {{0, 0, 0, 0, 1}, -1617.2681, 0.0000, -14.9471, 0.0000, 12.6153, 0.0000}};
594 
595  // number of coefficients
596  static const int N = int(sizeof(C) / sizeof(Coeff));
597 
598  dUT = dld = dom = 0.0;
599  for (int i = 0; i < N; i++)
600  {
601  double arg(0.0);
602  for (int j = 0; j < 5; j++)
603  {
604  arg += double(C[i].nargs[j]) * args[j + 1]; // NB args[0] is not used
605  }
606  arg = ::fmod(arg, TWOPI);
607  double sarg = ::sin(arg);
608  double carg = ::cos(arg);
609 
610  dUT += C[i].dutsin * sarg + C[i].dutcos * carg;
611  dld += C[i].dldsin * sarg + C[i].dldcos * carg;
612  dom += C[i].domsin * sarg + C[i].domcos * carg;
613  }
614 
615  // change units to seconds, seconds/day, and radians/sec
616  dUT *= 1.e-4;
617  dld *= 1.e-5;
618  dom *= 1.e-14;
619 
620  } // end correctEarthRotationZonalTides()
621 
622  //---------------------------------------------------------------------------------
623  /* Corrections to Earth rotation due to zonal tides using USNO IERS2003
624  algorithm param args array of 6 fundamental arguments at T, units radians
625  NB args[0] (=GMST+pi) is not used
626  param dUT correction to UT1mUTC in seconds
627  param dld correction to length of day in seconds/day
628  param dom correction to yp in radians per second */
629  void correctEarthRotationZonalTides2003(const double args[6], double& dUT,
630  double& dld, double& dom)
631  {
632  // constants
633  static const double TWOPI(6.283185307179586476925287);
634 
635  /* luni-solar multipliers in argument L Lp F D Omega, 62 terms
636  and coefficients of DUT sin, cos; DLOD cos, sin; DOMEGA cos, sin */
637  typedef struct
638  {
639  int nargs[5];
640  double dutsin, dutcos, dldcos, dldsin, domcos,
641  domsin; // yes its sc cs cs
642  } Coeff;
643 
644  static const Coeff C[] = {
645  {{1, 0, 2, 2, 2}, -0.02, 0.00, 0.26, 0.00, -0.22, 0.00},
646  {{2, 0, 2, 0, 1}, -0.04, 0.00, 0.38, 0.00, -0.32, 0.00},
647  {{2, 0, 2, 0, 2}, -0.10, 0.00, 0.91, 0.00, -0.76, 0.00},
648  {{0, 0, 2, 2, 1}, -0.05, 0.00, 0.45, 0.00, -0.38, 0.00},
649  {{0, 0, 2, 2, 2}, -0.12, 0.00, 1.09, 0.01, -0.92, -0.01},
650  {{1, 0, 2, 0, 0}, -0.04, 0.00, 0.27, 0.00, -0.22, 0.00},
651  {{1, 0, 2, 0, 1}, -0.41, 0.00, 2.84, 0.02, -2.40, -0.01},
652  {{1, 0, 2, 0, 2}, -1.00, 0.01, 6.85, 0.04, -5.78, -0.03},
653  {{3, 0, 0, 0, 0}, -0.02, 0.00, 0.12, 0.00, -0.11, 0.00},
654  {{-1, 0, 2, 2, 1}, -0.08, 0.00, 0.54, 0.00, -0.46, 0.00},
655  {{-1, 0, 2, 2, 2}, -0.20, 0.00, 1.30, 0.01, -1.10, -0.01},
656  {{1, 0, 0, 2, 0}, -0.08, 0.00, 0.50, 0.00, -0.42, 0.00},
657  {{2, 0, 2, -2, 2}, 0.02, 0.00, -0.11, 0.00, 0.09, 0.00},
658  {{0, 1, 2, 0, 2}, 0.03, 0.00, -0.12, 0.00, 0.10, 0.00},
659  {{0, 0, 2, 0, 0}, -0.30, 0.00, 1.39, 0.01, -1.17, -0.01},
660  {{0, 0, 2, 0, 1}, -3.22, 0.02, 14.86, 0.09, -12.54, -0.08},
661  {{0, 0, 2, 0, 2}, -7.79, 0.05, 35.84, 0.22, -30.25, -0.18},
662  {{2, 0, 0, 0, -1}, 0.02, 0.00, -0.10, 0.00, 0.08, 0.00},
663  {{2, 0, 0, 0, 0}, -0.34, 0.00, 1.55, 0.01, -1.31, -0.01},
664  {{2, 0, 0, 0, 1}, 0.02, 0.00, -0.08, 0.00, 0.07, 0.00},
665  {{0, -1, 2, 0, 2}, -0.02, 0.00, 0.11, 0.00, -0.09, 0.00},
666  {{0, 0, 0, 2, -1}, 0.05, 0.00, -0.20, 0.00, 0.17, 0.00},
667  {{0, 0, 0, 2, 0}, -0.74, 0.00, 3.14, 0.02, -2.65, -0.02},
668  {{0, 0, 0, 2, 1}, -0.05, 0.00, 0.22, 0.00, -0.19, 0.00},
669  {{0, -1, 0, 2, 0}, -0.05, 0.00, 0.21, 0.00, -0.17, 0.00},
670  {{1, 0, 2, -2, 1}, 0.05, 0.00, -0.13, 0.00, 0.11, 0.00},
671  {{1, 0, 2, -2, 2}, 0.10, 0.00, -0.26, 0.00, 0.22, 0.00},
672  {{1, 1, 0, 0, 0}, 0.04, 0.00, -0.10, 0.00, 0.08, 0.00},
673  {{-1, 0, 2, 0, 0}, 0.05, 0.00, -0.11, 0.00, 0.09, 0.00},
674  {{-1, 0, 2, 0, 1}, 0.18, 0.00, -0.41, 0.00, 0.35, 0.00},
675  {{-1, 0, 2, 0, 2}, 0.44, 0.00, -1.02, -0.01, 0.86, 0.01},
676  {{1, 0, 0, 0, -1}, 0.54, 0.00, -1.23, -0.01, 1.04, 0.01},
677  {{1, 0, 0, 0, 0}, -8.33, 0.06, 18.99, 0.13, -16.03, -0.11},
678  {{1, 0, 0, 0, 1}, 0.55, 0.00, -1.25, -0.01, 1.05, 0.01},
679  {{0, 0, 0, 1, 0}, 0.05, 0.00, -0.11, 0.00, 0.09, 0.00},
680  {{1, -1, 0, 0, 0}, -0.06, 0.00, 0.12, 0.00, -0.10, 0.00},
681  {{-1, 0, 0, 2, -1}, 0.12, 0.00, -0.24, 0.00, 0.20, 0.00},
682  {{-1, 0, 0, 2, 0}, -1.84, 0.01, 3.63, 0.02, -3.07, -0.02},
683  {{-1, 0, 0, 2, 1}, 0.13, 0.00, -0.26, 0.00, 0.22, 0.00},
684  {{1, 0, -2, 2, -1}, 0.02, 0.00, -0.04, 0.00, 0.03, 0.00},
685  {{-1, -1, 0, 2, 0}, -0.09, 0.00, 0.16, 0.00, -0.13, 0.00},
686  {{0, 2, 2, -2, 2}, -0.06, 0.00, 0.04, 0.00, -0.03, 0.00},
687  {{0, 1, 2, -2, 1}, 0.03, 0.00, -0.02, 0.00, 0.01, 0.00},
688  {{0, 1, 2, -2, 2}, -1.91, 0.02, 0.98, 0.01, -0.83, -0.01},
689  {{0, 0, 2, -2, 0}, 0.26, 0.00, -0.09, 0.00, 0.08, 0.00},
690  {{0, 0, 2, -2, 1}, 1.18, -0.01, -0.42, 0.00, 0.35, 0.00},
691  {{0, 0, 2, -2, 2}, -49.06, 0.43, 16.88, 0.15, -14.25, -0.13},
692  {{0, 2, 0, 0, 0}, -0.20, 0.00, 0.07, 0.00, -0.06, 0.00},
693  {{2, 0, 0, -2, -1}, 0.05, 0.00, -0.02, 0.00, 0.01, 0.00},
694  {{2, 0, 0, -2, 0}, -0.56, 0.01, 0.17, 0.00, -0.14, 0.00},
695  {{2, 0, 0, -2, 1}, 0.04, 0.00, -0.01, 0.00, 0.01, 0.00},
696  {{0, -1, 2, -2, 1}, -0.05, 0.00, 0.01, 0.00, -0.01, 0.00},
697  {{0, 1, 0, 0, -1}, 0.09, 0.00, -0.02, 0.00, 0.01, 0.00},
698  {{0, -1, 2, -2, 2}, 0.82, -0.01, -0.14, 0.00, 0.12, 0.00},
699  {{0, 1, 0, 0, 0}, -15.65, 0.15, 2.69, 0.03, -2.27, -0.02},
700  {{0, 1, 0, 0, 1}, -0.14, 0.00, 0.02, 0.00, -0.02, 0.00},
701  {{1, 0, 0, -1, 0}, 0.03, 0.00, 0.00, 0.00, 0.00, 0.00},
702  {{2, 0, -2, 0, 0}, -0.14, 0.00, -0.02, 0.00, 0.02, 0.00},
703  {{-2, 0, 2, 0, 1}, 0.43, -0.01, -0.02, 0.00, 0.02, 0.00},
704  {{-1, 1, 0, 1, 0}, -0.04, 0.00, 0.00, 0.00, 0.00, 0.00},
705  {{0, 0, 0, 0, 2}, 8.20, 0.11, 0.15, 0.00, -0.13, 0.00},
706  {{0, 0, 0, 0, 1}, -1689.54, -25.04, -15.62, 0.23, 13.18, -0.20}};
707 
708  // number of coefficients (62)
709  static const int N = int(sizeof(C) / sizeof(Coeff));
710 
711  dUT = dld = dom = 0.0;
712  for (int i = 0; i < N; i++)
713  {
714  double arg(0.0);
715  for (int j = 0; j < 5; j++)
716  {
717  arg += double(C[i].nargs[j]) * args[j + 1]; // NB args[0] is not used
718  }
719 
720  arg = ::fmod(arg, TWOPI);
721  double sarg = ::sin(arg);
722  double carg = ::cos(arg);
723 
724  dUT += C[i].dutsin * sarg + C[i].dutcos * carg;
725  dld += C[i].dldsin * sarg + C[i].dldcos * carg;
726  dom += C[i].domsin * sarg + C[i].domcos * carg;
727  }
728 
729  // change units to seconds, seconds/day, and radians/sec
730  dUT *= 1.e-4;
731  dld *= 1.e-5;
732  dom *= 1.e-14;
733 
734  } // end correctEarthRotationZonalTides2003()
735 
736  //---------------------------------------------------------------------------------
737  /* Corrections to UT1 and length of day (LOD) due to subdiurnal librations
738  using USNO IERS2010 algorithm.
739  param args array of 6 fundamental arguments at T (coordTransTime), in
740  radians param dUT correction to UT1mUTC in seconds param dld correction to
741  length of day in seconds/day */
742  void correctEarthRotationLibrations(const double args[6], double& dUT,
743  double& dld)
744  {
745  // static const double TWOPI(6.283185307179586476925287);
746  // static const double ARCSEC_TO_RAD=4.848136811095359935899141e-6;
747  // static const double ARCSEC_PER_CIRCLE(1296000.0);
748 
749  /* Coefficients of the quasi semidiurnal terms in dUT1, dLOD
750  IERS Conventions (2010), Table 5.1b */
751  typedef struct
752  {
753  int nargs[6];
754  double period, dutsin, dutcos, dldsin, dldcos; // NB period is not used
755  } Coeff;
756 
757  static const Coeff C[] = {
758  {{2, -2, 0, -2, 0, -2}, 0.5377239, 0.05, -0.03, -0.3, -0.6},
759  {{2, 0, 0, -2, -2, -2}, 0.5363232, 0.06, -0.03, -0.4, -0.7},
760  {{2, -1, 0, -2, 0, -2}, 0.5274312, 0.35, -0.20, -2.4, -4.1},
761  {{2, 1, 0, -2, -2, -2}, 0.5260835, 0.07, -0.04, -0.5, -0.8},
762  {{2, 0, 0, -2, 0, -1}, 0.5175645, -0.07, 0.04, 0.5, 0.8},
763  {{2, 0, 0, -2, 0, -2}, 0.5175251, 1.75, -1.01, -12.2, -21.3},
764  {{2, 1, 0, -2, 0, -2}, 0.5079842, -0.05, 0.03, 0.3, 0.6},
765  {{2, 0, -1, -2, 2, -2}, 0.5006854, 0.04, -0.03, -0.3, -0.6},
766  {{2, 0, 0, -2, 2, -2}, 0.5000000, 0.76, -0.44, -5.5, -9.6},
767  {{2, 0, 0, 0, 0, 0}, 0.4986348, 0.21, -0.12, -1.5, -2.6},
768  {{2, 0, 0, 0, 0, -1}, 0.4985982, 0.06, -0.04, -0.4, -0.8}};
769 
770  // number of coefficients (11)
771  static const int N = int(sizeof(C) / sizeof(Coeff));
772 
773  dUT = dld = 0.0;
774  for (int i = 0; i < N; ++i)
775  {
776  double arg = 0.0;
777  for (int j = 0; j < 6; ++j)
778  {
779  arg += C[i].nargs[j] * args[j];
780  }
781  arg = ::fmod(arg, EarthOrientation::TWOPI);
782  double sarg = ::sin(arg);
783  double carg = ::cos(arg);
784 
785  dUT += C[i].dutsin * sarg + C[i].dutcos * carg;
786  dld += C[i].dldsin * sarg + C[i].dldcos * carg;
787  }
788 
789  // convert to seconds
790  dUT *= 1.e-6; // seconds
791  dld *= 1.e-6; // seconds per day
792 
793  } // end correctEarthRotationLibrations()
794 
795  //---------------------------------------------------------------------------------
796  // private member functions, used internally
797 
798  //---------------------------------------------------------------------------------
799  /* Compute the 'coordinate transformation time' which is the time since
800  epoch J2000 = January 1 2000 12h UT = 2451545.0JD, divided by 36525 days.
801  This quantity is used throughout the terrestrial / inertial coordinate
802  transformations.
803  Throws if the TimeSystem conversion fails (if
804  t.system==TimeSystem::Unknown) */
805  double EarthOrientation::coordTransTime(EphTime ttag)
806  {
807  try
808  {
809  EphTime t(ttag);
810  t.convertSystemTo(TimeSystem::TT);
811  // return (t.dMJD()-JulianEpoch)/36525.0;
812 
813  // do in a way to maximize precision
814  int in = int(t.dMJD() - 0.5) - intJulianEpoch;
815  double frac = 0.5 + t.secOfDay() / 86400.0;
816  if (frac > 1.0)
817  {
818  frac -= 1.0;
819  }
820  return (double(in) / 36525.0 + frac / 36525.0);
821  }
822  catch (Exception& e)
823  {
824  GNSSTK_RETHROW(e);
825  }
826  }
827 
828  //---------------------------------------------------------------------------------
829  /* Return mean longitude of lunar ascending node, in radians,
830  given T, the coordTransTime at the Epoch of interest. (Ref: F5 pg 23)
831  valid for IERS1996 */
832  double EarthOrientation::Omega(double T)
833  {
834  return ::fmod(
835  450160.398036 +
836  T * (-6962890.2665 // diff Omega2003 only in .2665 vs .5431
837  + T * (7.4722 + T * (0.007702 + T * (-0.00005939)))),
838  ARCSEC_PER_CIRCLE) *
839  ARCSEC_TO_RAD;
840  }
841 
842  //---------------------------------------------------------------------------------
843  /* Return mean longitude of lunar ascending node, in radians,
844  given T, the coordTransTime at the Epoch of interest.
845  valid for IERS 2003, 2010 */
846  double EarthOrientation::Omega2003(double T)
847  {
848  return ::fmod(
849  450160.398036 // 125.04455501 * 3600
850  + T * (-6962890.5431 // .2665
851  + T * (7.4722 + T * (0.007702 + T * (-0.00005939)))),
852  ARCSEC_PER_CIRCLE) *
853  ARCSEC_TO_RAD;
854  }
855 
856  //---------------------------------------------------------------------------------
857  /* Return mean longitude of the moon - Omega, in radians,
858  given T, the coordTransTime at the Epoch of interest. (Ref: F3 pg 23)
859  valid for IERS1996, IERS2003, IERS2010 */
860  double EarthOrientation::F(double T)
861  {
862  return ::fmod(
863  335779.526232 // 93.27209062 * 3600
864  + T * (1739527262.8478 +
865  T * (-12.7512 + T * (-0.001037 + T * (0.00000417)))),
866  ARCSEC_PER_CIRCLE) *
867  ARCSEC_TO_RAD;
868  }
869 
870  //---------------------------------------------------------------------------------
871  /* Return mean elongation of the moon from the sun, in radians,
872  given T, the coordTransTime at the Epoch of interest. (Ref: F4 pg 23)
873  valid for IERS1996, IERS2003 */
874  double EarthOrientation::D(double T)
875  {
876  return ::fmod(1072260.703692 // 297.85019547 * 3600
877  +
878  T * (1602961601.2090 +
879  T * (-6.3706 + T * (0.006593 + T * (-0.00003169)))),
880  ARCSEC_PER_CIRCLE) *
881  ARCSEC_TO_RAD;
882  }
883 
884  //---------------------------------------------------------------------------------
885  /* Return mean anomaly of the moon, in radians,
886  given T, the coordTransTime at the Epoch of interest. (Ref: F1 pg 23)
887  valid for IERS1996, IERS2003 */
888  double EarthOrientation::L(double T)
889  {
890  return ::fmod(485868.249036 // 134.96340251 * 3600
891  +
892  T * (1717915923.2178 +
893  T * (31.8792 + T * (0.051635 + T * (-0.00024470)))),
894  ARCSEC_PER_CIRCLE) *
895  ARCSEC_TO_RAD;
896  }
897 
898  //---------------------------------------------------------------------------------
899  /* Return mean anomaly of the sun, in radians, given T, the coordTransTime at
900  the time of interest. (Ref: F2 pg 23) valid for IERS1996, IERS2003 */
901  double EarthOrientation::Lp(double T)
902  {
903  return ::fmod(1287104.793048 // 357.52910918 * 3600
904  +
905  T * (129596581.0481 +
906  T * (-0.5532 + T * (0.000136 // NB this has a minus
907  // sign in interp.f
908  + T * (-0.00001149)))),
909  ARCSEC_PER_CIRCLE) *
910  ARCSEC_TO_RAD;
911  }
912 
913  //------------------------------------------------------------------------------
914  /* Compute the mean longitude of Mercury, in radians, given T, the
915  coordTransTime at the time of interest. Valid for IERS2003, IERS2010 param
916  T coordinate transformation time. return LMe in radians. */
917  double EarthOrientation::LMe(double T)
918  {
919  double lme(::fmod(4.402608842 + 2608.7903141574 * T, TWOPI));
920  return lme;
921  }
922 
923  //------------------------------------------------------------------------------
924  /* Compute the mean longitude of Venus, in radians, given T, the
925  coordTransTime at the time of interest. Valid for IERS2003, IERS2010 param
926  T coordinate transformation time. return LV in radians. */
927  double EarthOrientation::LV(double T)
928  {
929  double lv(::fmod(3.176146697 + 1021.3285546211 * T, TWOPI));
930  return lv;
931  }
932 
933  //------------------------------------------------------------------------------
934  /* Compute the mean longitude of Earth, in radians, given T, the
935  coordTransTime at the time of interest. Valid for IERS2003, IERS2010 param
936  T coordinate transformation time. return LE in radians. */
937  double EarthOrientation::LE(double T)
938  {
939  double le(::fmod(1.753470314 + 628.3075849991 * T, TWOPI));
940  return le;
941  }
942 
943  //------------------------------------------------------------------------------
944  /* Compute the mean longitude of Mars, in radians, given T, the
945  coordTransTime at the time of interest. Valid for IERS2003, IERS2010 param
946  T coordinate transformation time. return LMa in radians. */
947  double EarthOrientation::LMa(double T)
948  {
949  double lma(::fmod(6.203480913 + 334.0612426700 * T, TWOPI));
950  return lma;
951  }
952 
953  //------------------------------------------------------------------------------
954  /* Compute the mean longitude of Jupiter, in radians, given T, the
955  coordTransTime at the time of interest. Valid for IERS2003, IERS2010
956  param T coordinate transformation time.
957  return LJ in radians. */
958  double EarthOrientation::LJ(double T)
959  {
960  double lj(::fmod(0.599546497 + 52.9690962641 * T, TWOPI));
961  return lj;
962  }
963 
964  //------------------------------------------------------------------------------
965  /* Compute the mean longitude of Saturn, in radians, given T, the
966  coordTransTime at the time of interest. Valid for IERS2003, IERS2010
967  param T coordinate transformation time.
968  return LS in radians. */
969  double EarthOrientation::LS(double T)
970  {
971  double ls(::fmod(0.874016757 + 21.3299104960 * T, TWOPI));
972  return ls;
973  }
974 
975  //------------------------------------------------------------------------------
976  /* Compute the mean longitude of Uranus, in radians, given T, the
977  coordTransTime at the time of interest. Valid for IERS2003, IERS2010
978  param T coordinate transformation time.
979  return LU in radians. */
980  double EarthOrientation::LU(double T)
981  {
982  double lu(::fmod(5.481293872 + 7.4781598567 * T, TWOPI));
983  return lu;
984  }
985 
986  //------------------------------------------------------------------------------
987  /* Compute the mean longitude of Neptune, in radians, given T, the
988  coordTransTime at the time of interest. Valid for IERS2003, IERS2010
989  param T coordinate transformation time.
990  return LN in radians. */
991  double EarthOrientation::LN(double T)
992  {
993  double ln(::fmod(5.311886287 + 3.8133035638 * T, TWOPI));
994  return ln;
995  }
996 
997  //------------------------------------------------------------------------------
998  /* Compute the general precession in longitude, in radians,
999  given T, the coordTransTime at the time of interest.
1000  Valid for IERS2003
1001  param T coordinate transformation time.
1002  return Pa in radians. */
1003  double EarthOrientation::Pa(double T)
1004  {
1005  double pa((0.024381750 + 0.00000538691 * T) * T);
1006  return pa;
1007  }
1008 
1009  //---------------------------------------------------------------------------------
1010  /* Compute eps, the obliquity of the ecliptic, in radians,
1011  given T, the coordTransTime at the time of interest.
1012  throw if convention is not defined */
1013  double EarthOrientation::obliquity(double T)
1014  {
1015  if (convention == IERSConvention::IERS1996)
1016  {
1017  return obliquity1996(T);
1018  }
1019  else if (convention == IERSConvention::IERS2003)
1020  {
1021  return obliquity1996(T); // same as 96
1022  }
1023  else if (convention == IERSConvention::IERS2010)
1024  {
1025  return obliquity2010(T);
1026  }
1027  else
1028  {
1029  Exception e("IERS convention is not defined");
1030  GNSSTK_THROW(e);
1031  }
1032  }
1033 
1034  //---------------------------------------------------------------------------------
1035  /* Compute Greenwich Mean Sidereal Time, or the Greenwich hour angle of
1036  the mean vernal equinox (radians), given the coordinate time of interest,
1037  and this object's UT1-UTC (sec), which comes from the IERS bulletin.
1038  param t EphTime epoch of the rotation (UTC).
1039  param reduced, bool true when UT1mUTC is 'reduced', meaning assumes
1040  'no tides', as is the case with the NGA EOPs (default=F).
1041  throw if convention is not defined */
1042  double EarthOrientation::GMST(const EphTime& t, bool reduced)
1043  {
1044  try
1045  {
1046  if (convention == IERSConvention::IERS1996)
1047  {
1048  return GMST1996(t, UT1mUTC, reduced);
1049  }
1050  else if (convention == IERSConvention::IERS2003)
1051  {
1052  return GMST2003(t, UT1mUTC);
1053  }
1054  else if (convention == IERSConvention::IERS2010)
1055  {
1056  return GMST2010(t, UT1mUTC);
1057  }
1058  else
1059  {
1060  Exception e("IERS convention is not defined");
1061  GNSSTK_THROW(e);
1062  }
1063  }
1064  catch (Exception& e)
1065  {
1066  GNSSTK_RETHROW(e);
1067  }
1068  }
1069 
1070  //---------------------------------------------------------------------------------
1071  /* Compute Greenwich Apparent Sidereal Time, or the Greenwich hour angle of
1072  the true vernal equinox (radians), given the coordinate time of interest,
1073  and this object's UT1-UTC (sec), which comes from the IERS bulletin.
1074  param t EphTime epoch of the rotation.
1075  param reduced, bool true when UT1mUTC is 'reduced', meaning assumes
1076  'no tides', as is the case with the NGA EOPs (default=F).
1077  throw if convention is not defined */
1078  double EarthOrientation::GAST(const EphTime& t, bool reduced)
1079  {
1080  try
1081  {
1082  if (convention == IERSConvention::IERS1996)
1083  {
1084  return GAST1996(t, UT1mUTC, reduced);
1085  }
1086  else if (convention == IERSConvention::IERS2003)
1087  {
1088  return GAST2003(t, UT1mUTC);
1089  }
1090  else if (convention == IERSConvention::IERS2010)
1091  {
1092  return GAST2010(t, UT1mUTC);
1093  }
1094  else
1095  {
1096  Exception e("IERS convention is not defined");
1097  GNSSTK_THROW(e);
1098  }
1099  }
1100  catch (Exception& e)
1101  {
1102  GNSSTK_RETHROW(e);
1103  }
1104  }
1105 
1106  //---------------------------------------------------------------------------------
1107  /* Generate transformation matrix (3X3 rotation) due to the EOPs (polar
1108  motion angles xp and yp (arcseconds)), as found in the IERS bulletin; see
1109  class EarthOrientation. param t EphTime epoch of the rotation. throw if
1110  convention is not defined */
1111  Matrix<double> EarthOrientation::polarMotionMatrix(const EphTime& t)
1112  {
1113  if (convention == IERSConvention::IERS1996)
1114  {
1115  return polarMotionMatrix1996(xp, yp);
1116  }
1117  else if (convention == IERSConvention::IERS2003)
1118  {
1119  return polarMotionMatrix2003(t, xp, yp);
1120  }
1121  else if (convention == IERSConvention::IERS2010)
1122  {
1123  return polarMotionMatrix2003(t, xp, yp); // valid also for 2010
1124  }
1125  else
1126  {
1127  Exception e("IERS convention is not defined");
1128  GNSSTK_THROW(e);
1129  }
1130  }
1131 
1132  //---------------------------------------------------------------------------------
1133  /* Compute the precession matrix, a 3x3 rotation matrix, given T,
1134  the coordinate transformation time at the time of interest
1135  throw if convention is not defined */
1136  Matrix<double> EarthOrientation::precessionMatrix(const EphTime& t)
1137  {
1138  double T = coordTransTime(t);
1139  if (convention == IERSConvention::IERS1996)
1140  {
1141  return precessionMatrix1996(T);
1142  }
1143  else if (convention == IERSConvention::IERS2003)
1144  {
1145  return precessionMatrix2003(T);
1146  }
1147  else if (convention == IERSConvention::IERS2010)
1148  {
1149  return precessionMatrix2010(T);
1150  }
1151  else
1152  {
1153  Exception e("IERS convention is not defined");
1154  GNSSTK_THROW(e);
1155  }
1156  }
1157 
1158  //------------------------------------------------------------------------------
1159  /* Nutation of the obliquity (deps) and of the longitude (dpsi)
1160  param T, the coordinate transformation time at the time of interest
1161  param deps, nutation of the obliquity (output) in radians
1162  param dpsi, nutation of the longitude (output) in radians
1163  throw if convention is not defined
1164  NB this is never called internally
1165  void EarthOrientation::nutationAngles(double T, double& deps, double&
1166  dpsi)
1167  {
1168  if(convention == IERSConvention::IERS1996) {
1169  double om;
1170  nutationAngles1996(T,deps,dpsi,om);
1171  }
1172  if(convention == IERSConvention::IERS2003)
1173  nutationAngles2003(T,deps,dpsi);
1174  if(convention == IERSConvention::IERS2010)
1175  nutationAngles2010(T,deps,dpsi);
1176  else {
1177  Exception e("IERS convention is not defined");
1178  GNSSTK_THROW(e);
1179  }
1180  } */
1181 
1182  //---------------------------------------------------------------------------------
1183  // Compute the nutation matrix, given coordinate transformation time T
1184  Matrix<double> EarthOrientation::nutationMatrix(const EphTime& t)
1185  {
1186  try
1187  {
1188  double T = coordTransTime(t);
1189  if (convention == IERSConvention::IERS1996)
1190  {
1191  return nutationMatrix1996(T);
1192  }
1193  else if (convention == IERSConvention::IERS2003)
1194  {
1195  return nutationMatrix2003(T);
1196  }
1197  else if (convention == IERSConvention::IERS2010)
1198  {
1199  return nutationMatrix2010(T);
1200  }
1201  else
1202  {
1203  Exception e("IERS convention is not defined");
1204  GNSSTK_THROW(e);
1205  }
1206  }
1207  catch (Exception& e)
1208  {
1209  GNSSTK_RETHROW(e);
1210  }
1211  }
1212 
1213  //------------------------------------------------------------------------------
1214  /* Generate precise transformation matrix (3X3 rotation) for Earth motion due
1215  to precession, nutation and frame bias, at the given time of interest.
1216  param t EphTime epoch of the rotation.
1217  return 3x3 rotation matrix
1218  throw if the TimeSystem conversion fails (if
1219  t.system==TimeSystem::Unknown) throw if convention is not defined */
1220  Matrix<double> EarthOrientation::preciseEarthRotation(const EphTime& t)
1221  {
1222  if (convention == IERSConvention::IERS1996)
1223  {
1224  return Matrix<double>(3, 3); // TD
1225  }
1226  else if (convention == IERSConvention::IERS2003)
1227  {
1228  return preciseEarthRotation2003(coordTransTime(t));
1229  }
1230  else if (convention == IERSConvention::IERS2010)
1231  {
1232  return preciseEarthRotation2010(coordTransTime(t));
1233  }
1234  else
1235  {
1236  Exception e("IERS convention is not defined");
1237  GNSSTK_THROW(e);
1238  }
1239  }
1240 
1241  //---------------------------------------------------------------------------------
1242  /* Generate the full transformation matrix (3x3 rotation) relating the ECEF
1243  frame to the conventional inertial frame.
1244  Input is the time of interest; use this object's EOPs - the polar motion
1245  angles xp and yp (arcseconds), and UT1-UTC (seconds) (as found in the IERS
1246  bulletin). throw if convention is not defined */
1247  Matrix<double> EarthOrientation::ECEFtoInertial(const EphTime& t,
1248  bool reduced)
1249  {
1250  try
1251  {
1252  if (convention == IERSConvention::IERS1996)
1253  {
1254  return ECEFtoInertial1996(t, xp, yp, UT1mUTC, reduced);
1255  }
1256  else if (convention == IERSConvention::IERS2003)
1257  {
1258  return ECEFtoInertial2003(t, xp, yp, UT1mUTC);
1259  }
1260  else if (convention == IERSConvention::IERS2010)
1261  {
1262  return ECEFtoInertial2010(t, xp, yp, UT1mUTC);
1263  }
1264  else
1265  {
1266  Exception e("IERS convention is not defined");
1267  GNSSTK_THROW(e);
1268  }
1269  }
1270  catch (Exception& e)
1271  {
1272  GNSSTK_RETHROW(e);
1273  }
1274  }
1275 
1276  //------------------------------------------------------------------------------
1277  /* Compute the transformation from ECEF to the J2000 dynamical (inertial)
1278  frame. This differs from the ECEFtoInertial transformation only by the
1279  frame bias matrix. Only available in IERS2010.
1280 
1281  TD this is not right
1282  Matrix<double> EarthOrientation::ECEFtoJ2000(const EphTime& t, bool
1283  reduced)
1284  {
1285  try {
1286  if(convention != IERSConvention::IERS2010) {
1287  Exception e("ECEFtoJ2000 implemented only for IERS2010");
1288  GNSSTK_THROW(e);
1289  }
1290 
1291  // get the frame bias matrix
1292  //double T=0.0; // T at J2000
1293  double gamb, phib, psib, eps;
1294  //FukushimaWilliams(T, gamb, phib, psib, eps);
1295  //eps = 0.0;
1296  //Matrix<double> FB = FukushimaWilliams(gamb, phib, psib, eps);
1297  gamb = 0.052928 * ARCSEC_TO_RAD;
1298  phib = 0.006891 * ARCSEC_TO_RAD;
1299  psib = 0.041775 * ARCSEC_TO_RAD;
1300  Matrix<double> FB =
1301  rotation(-psib,3)*rotation(phib,1)*rotation(gamb,3);
1302 
1303  // get the full ECEFtoInertial, and back out the frame bias
1304  Matrix<double> T2C = ECEFtoInertial2010(t, xp, yp, UT1mUTC);
1305 
1306  return (FB*T2C);
1307  }
1308  catch(Exception& e) { GNSSTK_RETHROW(e); }
1309  }
1310  */
1311 
1312  //------------------------------------------------------------------------------
1313  // private functions
1314  //------------------------------------------------------------------------------
1315 
1316  //------------------------------------------------------------------------------
1317  /* Compute the locator s which gives the position of the CIO on the equator
1318  of the CIP, given the coordinate transformation time T. Consistent with
1319  IAU 2000A (IERS2003) precession-nutation, and P03 precession (IERS2010).
1320  and the coordinates X,Y of the CIP. Derived in part from SOFA routine
1321  s00.c for IERS2003 and s06.c for IERS2010. param T, the coordinate
1322  transformation time at the time of interest param X, the X coordinate of
1323  the CIP (input) param Y, the Y coordinate of the CIP (input) param which,
1324  the IERS convention to be used (input) return S, the parameter that
1325  positions the CIO on the CIP equator in radians.
1326  */
1327  double EarthOrientation::S(double T, double& X, double& Y,
1328  IERSConvention which)
1329  {
1330  int i, j;
1331  double s, st[6], arg;
1332 
1333  // Fundamental arguments: all in radians
1334  double farg[8];
1335  farg[0] = L(T); // mean anomaly of the moon
1336  farg[1] = Lp(T); // mean anomaly of the sun
1337  farg[2] = F(T); // mean longitude of the moon minus Omega
1338  farg[3] = D(T); // mean elongation of the moon from the sun
1339  // mean longitude of lunar ascending node
1340  if (which == IERSConvention::IERS2010)
1341  {
1342  farg[4] = Omega2003(T);
1343  }
1344  else
1345  {
1346  farg[4] = Omega(T);
1347  }
1348  farg[5] = LV(T); // mean longitude of Venus
1349  farg[6] = LE(T); // mean longitude of Earth
1350  farg[7] = Pa(T); // general precession in longitude
1351 
1352  // Based on table 5.2c IERS Tech Note 32 Chap 5 (tab5.2c.txt obtained from
1353  // IERS) also used SOFA routine s00.c
1354 
1355  // coefficients of polynomial in T for IERS2003
1356  static const double polycoeff[6] = {94.00e-6, 3808.35e-6, -119.94e-6,
1357  -72574.09e-6, 27.70e-6, 15.61e-6};
1358  // same for 2010
1359  static const double polycoeff2010[6] = {
1360  94.00e-6, 3808.65e-6, -122.68e-6, -72574.11e-6, 27.98e-6, 15.62e-6};
1361 
1362  typedef struct
1363  {
1364  int coeff[8]; // coefficients of l,lp,f,d,o,lv,le,pa
1365  double sincoeff, coscoeff; // coefficients of sine and cosine
1366  } Coeffs;
1367 
1368  // constant terms (T^0)
1369  static const Coeffs C0[] = {
1370  // indexes 1-10
1371  {{0, 0, 0, 0, 1, 0, 0, 0}, -2640.73e-6, 0.39e-6},
1372  {{0, 0, 0, 0, 2, 0, 0, 0}, -63.53e-6, 0.02e-6},
1373  {{0, 0, 2, -2, 3, 0, 0, 0}, -11.75e-6, -0.01e-6},
1374  {{0, 0, 2, -2, 1, 0, 0, 0}, -11.21e-6, -0.01e-6},
1375  {{0, 0, 2, -2, 2, 0, 0, 0}, 4.57e-6, 0.00e-6},
1376  {{0, 0, 2, 0, 3, 0, 0, 0}, -2.02e-6, 0.00e-6},
1377  {{0, 0, 2, 0, 1, 0, 0, 0}, -1.98e-6, 0.00e-6},
1378  {{0, 0, 0, 0, 3, 0, 0, 0}, 1.72e-6, 0.00e-6},
1379  {{0, 1, 0, 0, 1, 0, 0, 0}, 1.41e-6, 0.01e-6},
1380  {{0, 1, 0, 0, -1, 0, 0, 0}, 1.26e-6, 0.01e-6},
1381  // indexes 11-20
1382  {{1, 0, 0, 0, -1, 0, 0, 0}, 0.63e-6, 0.00e-6},
1383  {{1, 0, 0, 0, 1, 0, 0, 0}, 0.63e-6, 0.00e-6},
1384  {{0, 1, 2, -2, 3, 0, 0, 0}, -0.46e-6, 0.00e-6},
1385  {{0, 1, 2, -2, 1, 0, 0, 0}, -0.45e-6, 0.00e-6},
1386  {{0, 0, 4, -4, 4, 0, 0, 0}, -0.36e-6, 0.00e-6},
1387  {{0, 0, 1, -1, 1, -8, 12, 0}, 0.24e-6, 0.12e-6},
1388  {{0, 0, 2, 0, 0, 0, 0, 0}, -0.32e-6, 0.00e-6},
1389  {{0, 0, 2, 0, 2, 0, 0, 0}, -0.28e-6, 0.00e-6},
1390  {{1, 0, 2, 0, 3, 0, 0, 0}, -0.27e-6, 0.00e-6},
1391  {{1, 0, 2, 0, 1, 0, 0, 0}, -0.26e-6, 0.00e-6},
1392  // indexes 21-30
1393  {{0, 0, 2, -2, 0, 0, 0, 0}, 0.21e-6, 0.00e-6},
1394  {{0, 1, -2, 2, -3, 0, 0, 0}, -0.19e-6, 0.00e-6},
1395  {{0, 1, -2, 2, -1, 0, 0, 0}, -0.18e-6, 0.00e-6},
1396  {{0, 0, 0, 0, 0, 8, -13, -1}, 0.10e-6, -0.05e-6},
1397  {{0, 0, 0, 2, 0, 0, 0, 0}, -0.15e-6, 0.00e-6},
1398  {{2, 0, -2, 0, -1, 0, 0, 0}, 0.14e-6, 0.00e-6},
1399  {{0, 1, 2, -2, 2, 0, 0, 0}, 0.14e-6, 0.00e-6},
1400  {{1, 0, 0, -2, 1, 0, 0, 0}, -0.14e-6, 0.00e-6},
1401  {{1, 0, 0, -2, -1, 0, 0, 0}, -0.14e-6, 0.00e-6},
1402  {{0, 0, 4, -2, 4, 0, 0, 0}, -0.13e-6, 0.00e-6},
1403  // indexes 31-33
1404  {{0, 0, 2, -2, 4, 0, 0, 0}, 0.11e-6, 0.00e-6},
1405  {{1, 0, -2, 0, -3, 0, 0, 0}, -0.11e-6, 0.00e-6},
1406  {{1, 0, -2, 0, -1, 0, 0, 0}, -0.11e-6, 0.00e-6}};
1407 
1408  // First order terms (T)
1409  // NB C1[1].sincoeff=1.71e-6 in 2003 becomes 1.73e-6 in 2010 (2nd row)
1410  static const double C1_1_sincoeff2010 = 1.73e-6;
1411  static const Coeffs C1[] = {
1412  {{0, 0, 0, 0, 2, 0, 0, 0}, -0.07e-6, 3.57e-6},
1413  {{0, 0, 0, 0, 1, 0, 0, 0}, 1.71e-6, -0.03e-6},
1414  {{0, 0, 2, -2, 3, 0, 0, 0}, 0.00e-6, 0.48e-6}};
1415 
1416  // Second order terms (T^2)
1417  // NB C2[0].sincoeff=743.53e-6 in 2003 becomes 743.52e-6 in 2010 (1st row)
1418  static const double C2_0_sincoeff = 743.52e-6;
1419  static const Coeffs C2[] = {
1420  // indexes 1-10
1421  {{0, 0, 0, 0, 1, 0, 0, 0}, 743.53e-6, -0.17e-6},
1422  {{0, 0, 2, -2, 2, 0, 0, 0}, 56.91e-6, 0.06e-6},
1423  {{0, 0, 2, 0, 2, 0, 0, 0}, 9.84e-6, -0.01e-6},
1424  {{0, 0, 0, 0, 2, 0, 0, 0}, -8.85e-6, 0.01e-6},
1425  {{0, 1, 0, 0, 0, 0, 0, 0}, -6.38e-6, -0.05e-6},
1426  {{1, 0, 0, 0, 0, 0, 0, 0}, -3.07e-6, 0.00e-6},
1427  {{0, 1, 2, -2, 2, 0, 0, 0}, 2.23e-6, 0.00e-6},
1428  {{0, 0, 2, 0, 1, 0, 0, 0}, 1.67e-6, 0.00e-6},
1429  {{1, 0, 2, 0, 2, 0, 0, 0}, 1.30e-6, 0.00e-6},
1430  {{0, 1, -2, 2, -2, 0, 0, 0}, 0.93e-6, 0.00e-6},
1431  // indexes 11-20
1432  {{1, 0, 0, -2, 0, 0, 0, 0}, 0.68e-6, 0.00e-6},
1433  {{0, 0, 2, -2, 1, 0, 0, 0}, -0.55e-6, 0.00e-6},
1434  {{1, 0, -2, 0, -2, 0, 0, 0}, 0.53e-6, 0.00e-6},
1435  {{0, 0, 0, 2, 0, 0, 0, 0}, -0.27e-6, 0.00e-6},
1436  {{1, 0, 0, 0, 1, 0, 0, 0}, -0.27e-6, 0.00e-6},
1437  {{1, 0, -2, -2, -2, 0, 0, 0}, -0.26e-6, 0.00e-6},
1438  {{1, 0, 0, 0, -1, 0, 0, 0}, -0.25e-6, 0.00e-6},
1439  {{1, 0, 2, 0, 1, 0, 0, 0}, 0.22e-6, 0.00e-6},
1440  {{2, 0, 0, -2, 0, 0, 0, 0}, -0.21e-6, 0.00e-6},
1441  {{2, 0, -2, 0, -1, 0, 0, 0}, 0.20e-6, 0.00e-6},
1442  // indexes 21-25
1443  {{0, 0, 2, 2, 2, 0, 0, 0}, 0.17e-6, 0.00e-6},
1444  {{2, 0, 2, 0, 2, 0, 0, 0}, 0.13e-6, 0.00e-6},
1445  {{2, 0, 0, 0, 0, 0, 0, 0}, -0.13e-6, 0.00e-6},
1446  {{1, 0, 2, -2, 2, 0, 0, 0}, -0.12e-6, 0.00e-6},
1447  {{0, 0, 2, 0, 0, 0, 0, 0}, -0.11e-6, 0.00e-6}};
1448 
1449  // Third order terms (T^3)
1450  static const Coeffs C3[] = {
1451  {{0, 0, 0, 0, 1, 0, 0, 0}, 0.30e-6, -23.51e-6},
1452  {{0, 0, 2, -2, 2, 0, 0, 0}, -0.03e-6, -1.39e-6},
1453  {{0, 0, 2, 0, 2, 0, 0, 0}, -0.01e-6, -0.24e-6},
1454  {{0, 0, 0, 0, 2, 0, 0, 0}, 0.00e-6, 0.22e-6}};
1455  static const Coeffs C32010[] = {
1456  {{0, 0, 0, 0, 1, 0, 0, 0}, 0.30e-6, -23.42e-6},
1457  {{0, 0, 2, -2, 2, 0, 0, 0}, -0.03e-6, -1.46e-6},
1458  {{0, 0, 2, 0, 2, 0, 0, 0}, -0.01e-6, -0.25e-6},
1459  {{0, 0, 0, 0, 2, 0, 0, 0}, 0.00e-6, 0.23e-6}};
1460 
1461  // Fourth order terms (T^4)
1462  static const Coeffs C4[] = {
1463  {{0, 0, 0, 0, 1, 0, 0, 0}, -0.26e-6, -0.01e-6}};
1464 
1465  // number of terms in each order
1466  const int n0 = int(sizeof C0 / sizeof(Coeffs));
1467  const int n1 = int(sizeof C1 / sizeof(Coeffs));
1468  const int n2 = int(sizeof C2 / sizeof(Coeffs));
1469  const int n3 = int(sizeof C3 / sizeof(Coeffs));
1470  const int n4 = int(sizeof C4 / sizeof(Coeffs));
1471 
1472  // initialize with the polynomial coefficients
1473  if (which != IERSConvention::IERS2010)
1474  {
1475  for (i = 0; i < 6; ++i)
1476  {
1477  st[i] = polycoeff[i];
1478  }
1479  }
1480  else
1481  {
1482  for (i = 0; i < 6; ++i)
1483  {
1484  st[i] = polycoeff2010[i];
1485  }
1486  }
1487 
1488  // do the sums
1489  for (i = n0 - 1; i >= 0; --i)
1490  { // order 0
1491  arg = 0.0;
1492  for (j = 0; j < 8; ++j)
1493  {
1494  arg += C0[i].coeff[j] * farg[j];
1495  }
1496  st[0] += C0[i].sincoeff * ::sin(arg) + C0[i].coscoeff * ::cos(arg);
1497  }
1498 
1499  for (i = n1 - 1; i >= 0; --i)
1500  { // order 1
1501  arg = 0.0;
1502  for (j = 0; j < 8; ++j)
1503  {
1504  arg += C1[i].coeff[j] * farg[j];
1505  }
1506  if (which == IERSConvention::IERS2010 && i == 1)
1507  {
1508  st[1] +=
1509  C1_1_sincoeff2010 * ::sin(arg) + C1[i].coscoeff * ::cos(arg);
1510  }
1511  else
1512  {
1513  st[1] += C1[i].sincoeff * ::sin(arg) + C1[i].coscoeff * ::cos(arg);
1514  }
1515  }
1516 
1517  for (i = n2 - 1; i >= 0; --i)
1518  { // order 2
1519  arg = 0.0;
1520  for (j = 0; j < 8; ++j)
1521  {
1522  arg += C2[i].coeff[j] * farg[j];
1523  }
1524  if (which == IERSConvention::IERS2010 && i == 0)
1525  {
1526  st[2] += C2_0_sincoeff * ::sin(arg) + C2[i].coscoeff * ::cos(arg);
1527  }
1528  else
1529  {
1530  st[2] += C2[i].sincoeff * ::sin(arg) + C2[i].coscoeff * ::cos(arg);
1531  }
1532  }
1533 
1534  // order 3
1535  if (which != IERSConvention::IERS2010) // 2003
1536  {
1537  for (i = n3 - 1; i >= 0; --i)
1538  {
1539  arg = 0.0;
1540  for (j = 0; j < 8; ++j)
1541  {
1542  arg += C3[i].coeff[j] * farg[j];
1543  }
1544  st[3] += C3[i].sincoeff * ::sin(arg) + C3[i].coscoeff * ::cos(arg);
1545  }
1546  }
1547  else // 2010
1548  {
1549  for (i = n3 - 1; i >= 0; --i)
1550  {
1551  arg = 0.0;
1552  for (j = 0; j < 8; ++j)
1553  {
1554  arg += C32010[i].coeff[j] * farg[j];
1555  }
1556  st[3] += C32010[i].sincoeff * ::sin(arg) +
1557  C32010[i].coscoeff * ::cos(arg);
1558  }
1559  }
1560 
1561  for (i = n4 - 1; i >= 0; --i)
1562  { // order 4
1563  arg = 0.0;
1564  for (j = 0; j < 8; ++j)
1565  {
1566  arg += C4[i].coeff[j] * farg[j];
1567  }
1568  st[4] += C4[i].sincoeff * ::sin(arg) + C4[i].coscoeff * ::cos(arg);
1569  }
1570 
1571  // combine all the terms
1572  s = st[0] +
1573  (st[1] + (st[2] + (st[3] + (st[4] + st[5] * T) * T) * T) * T) * T;
1574  s *= ARCSEC_TO_RAD;
1575  s -= X * Y / 2.0;
1576 
1577  return s;
1578  }
1579 
1580  //------------------------------------------------------------------------------
1581  /* The position of the Terrestrial Ephemeris Origin (TEO) on the equator of
1582  the Celestial Intermediate Pole (CIP), as given by the quantity s'. Also
1583  called the Terrestrial Intermediate Origin (TIO). Valid 2003 and 2010 Ref.
1584  IERS Tech Note 32 Chap 5 Eqn 12 and IERS Tech Note 36 Chap 5 Eqn 5.13
1585  param T Coordinate transformation time T.
1586  return angle 's prime' in radians */
1587  double EarthOrientation::Sprime(double T)
1588  {
1589  double sp = -47.0e-6 * T * ARCSEC_TO_RAD;
1590  return sp;
1591  }
1592 
1593  //------------------------------------------------------------------------------
1594  /* Compute the coordinates X,Y of the Celestial Intermediate Origin (CIO)
1595  using a series based on IAU 2006 precession and IAU 2000A nutation (IERS
1596  2010). The coordinates form a unit vector that points towards the CIO;
1597  they include the effects of frame bias, precession and nutation. cf. sofa
1598  xy06 Reference IERS(2010) Section 5.5.4 param T, the coordinate
1599  transformation time at the time of interest param X, x coordinate of CIO
1600  param Y, y coordinate of CIO */
1601  void EarthOrientation::XYCIO(double& T, double& X, double& Y)
1602  {
1603 // include data arrays : defines MAXPT
1604 #include "IERS2010CIOSeriesData.hpp"
1605 
1606  int i, j;
1607  double t;
1608 
1609  // compute and store powers of T
1610  double powsT[MAXPT + 1];
1611  t = 1.0;
1612  for (i = 0; i <= MAXPT; i++)
1613  {
1614  powsT[i] = t;
1615  t *= T;
1616  }
1617 
1618  // fundamental arguments
1619  double fa[14];
1620  fa[0] = L(T); // mean anomaly of the moon
1621  fa[1] = Lp(T); // mean anomaly of the sun
1622  fa[2] = F(T); // mean longitude of the moon - Omega
1623  fa[3] = D(T); // mean elongation of the moon from the sun
1624  fa[4] = Omega2003(T); // mean longitude of lunar ascending node
1625  fa[5] = LMe(T); // mean longitude Mercury
1626  fa[6] = LV(T); // mean longitude of Venus
1627  fa[7] = LE(T); // mean longitude of Earth
1628  fa[8] = LMa(T); // mean longitude Mars
1629  fa[9] = LJ(T); // mean longitude Jupiter
1630  fa[10] = LS(T); // mean longitude Saturn
1631  fa[11] = LU(T); // mean longitude Uranus
1632  fa[12] = LN(T); // mean longitude Neptune
1633  fa[13] = Pa(T); // general precession in longitude
1634 
1635  // intermediate totals
1636  double xypoly[2] = {0.0, 0.0}, xylunarsolar[2] = {0.0, 0.0},
1637  xyplanet[2] = {0.0, 0.0};
1638 
1639  // polynomial
1640  for (i = 0; i < 2; i++)
1641  {
1642  for (j = MAXPT; j >= 0; j--)
1643  {
1644  xypoly[i] += XYcoeff[i][j] * powsT[j];
1645  }
1646  }
1647 
1648  // nutation planetary terms
1649  int ilast(NAmp), jlast;
1650  double sc[2];
1651  for (int ifreq = NFAP - 1; ifreq >= 0; ifreq--)
1652  {
1653  // build the argument
1654  double arg(0.0);
1655  for (i = 0; i < 14; i++)
1656  {
1657  if (nFAplanetary[ifreq][i]) // don't add zero
1658  {
1659  arg += double(nFAplanetary[ifreq][i]) * fa[i];
1660  }
1661  }
1662 
1663  // store sin and cos
1664  sc[0] = ::sin(arg);
1665  sc[1] = ::cos(arg);
1666 
1667  // amplitudes
1668  jlast = iamp[ifreq + NFALS];
1669  for (i = ilast; i >= jlast; i--)
1670  {
1671  j = i - jlast; // coeff number
1672  xyplanet[jaxy[j]] += amp[i - 1] * sc[jasc[j]] * powsT[japt[j]];
1673  }
1674  ilast = jlast - 1;
1675  }
1676 
1677  // NB ilast is maintained through this point
1678 
1679  // nutation lunar solar terms
1680  for (int ifreq = NFALS - 1; ifreq >= 0; ifreq--)
1681  {
1682  // build the argument
1683  double arg(0.0);
1684  for (i = 0; i < 5; i++)
1685  {
1686  if (nFAlunarsolar[ifreq][i]) // don't add zero
1687  {
1688  arg += double(nFAlunarsolar[ifreq][i]) * fa[i];
1689  }
1690  }
1691 
1692  // store sin and cos
1693  sc[0] = ::sin(arg);
1694  sc[1] = ::cos(arg);
1695 
1696  // amplitudes
1697  jlast = iamp[ifreq];
1698  for (i = ilast; i >= jlast; i--)
1699  {
1700  j = i - jlast; // coeff number
1701  xylunarsolar[jaxy[j]] += amp[i - 1] * sc[jasc[j]] * powsT[japt[j]];
1702  }
1703  ilast = jlast - 1;
1704  }
1705 
1706  X = xypoly[0] + (xylunarsolar[0] + xyplanet[0]) * 1.e-6;
1707  X *= ARCSEC_TO_RAD;
1708  Y = xypoly[1] + (xylunarsolar[1] + xyplanet[1]) * 1.e-6;
1709  Y *= ARCSEC_TO_RAD;
1710  }
1711 
1712  //---------------------------------------------------------------------------------
1713  /* Starting with 2003 (and valid for 2010) conventions a new method for
1714  computing the transformation fron ITRS to GCRS is provided by the
1715  Celestial Ephemeris Origin (CEO) which is based on the Earth Rotation
1716  Angle. cf. sofa ERA00.c param t EphTime time of interest param UT1mUTC
1717  offset UT1-UTC return Earth rotation angle in radians throw if the
1718  TimeSystem conversion fails (if t.system==TimeSystem::Unknown) */
1719  double EarthOrientation::EarthRotationAngle(const EphTime& t,
1720  double& UT1mUTC)
1721  {
1722  try
1723  {
1724  EphTime tUT1(t);
1725  tUT1.convertSystemTo(TimeSystem::UTC);
1726  tUT1 += UT1mUTC;
1727 
1728  // TN36 eqn 5.15
1729  // double days = tUT1.dMJD() - JulianEpoch; // days since JE
1730  int idays(int(tUT1.dMJD() - 0.5) -
1731  intJulianEpoch); // days = idays+frac
1732  double frac(0.5 +
1733  tUT1.secOfDay() / 86400.0); // fractional part of days
1734  if (frac > 1.0)
1735  {
1736  frac -= 1.0;
1737  }
1738 
1739  // mod the terms with 1 individually to avoid numerical error
1740  double term1 = frac + 0.7790572732640 + 0.00273781191135448 * frac;
1741  if (term1 > 1.0)
1742  {
1743  term1 -= 1.0;
1744  }
1745  double term2 = ::fmod(0.00273781191135448 * double(idays), 1.0);
1746  double term = ::fmod(term1 + term2, 1.0);
1747 
1748  double era = TWOPI * term;
1749  if (era > TWOPI)
1750  {
1751  era -= TWOPI;
1752  }
1753 
1754  return era;
1755  }
1756  catch (Exception& e)
1757  {
1758  GNSSTK_RETHROW(e);
1759  }
1760  }
1761 
1762  //------------------------------------------------------------------------------
1763  /* Equation of the equinoxes complementary terms, IAU 2000 (IERS 2003)
1764  Note that GAST = GMST + equationOfEquinoxes2003
1765  param t EphTime epoch of interest
1766  return the ee in radians
1767  throw if the TimeSystem conversion fails (if
1768  t.system==TimeSystem::Unknown) Based on IERS function EECT2000.f; all
1769  planets but Venus dropped b/c their contribution is zero. */
1770  double EarthOrientation::equationOfEquinoxes2003(EphTime t)
1771  {
1772  // first define the coefficients
1773 
1774  // number of integer coefficients and fundamental arguments
1775  const int N(8);
1776 
1777  // NB drop the lma,lju,lsa,lur,lne terms - all zero!
1778  typedef struct
1779  {
1780  // coefficients of l,lp,f,d,o,lme,lv,le,lma,lju,lsa,lur,lne,pa in arg
1781  // *these were dropped: * * * * * *
1782  int coeff[N];
1783  // coefficients of sin(arg) and cos(arg) in the series
1784  double sincoeff, coscoeff;
1785  } Coeffs;
1786 
1787  // T^0 ...not very efficient
1788  static const Coeffs Czero[] = {
1789  // l lp f d o lv le pa c(sin) c(cos)
1790  // 1-10
1791  {{0, 0, 0, 0, 1, 0, 0, 0}, 2640.96e-6, -0.39e-6},
1792  {{0, 0, 0, 0, 2, 0, 0, 0}, 63.52e-6, -0.02e-6},
1793  {{0, 0, 2, -2, 3, 0, 0, 0}, 11.75e-6, 0.01e-6},
1794  {{0, 0, 2, -2, 1, 0, 0, 0}, 11.21e-6, 0.01e-6},
1795  {{0, 0, 2, -2, 2, 0, 0, 0}, -4.55e-6, 0.00e-6},
1796  {{0, 0, 2, 0, 3, 0, 0, 0}, 2.02e-6, 0.00e-6},
1797  {{0, 0, 2, 0, 1, 0, 0, 0}, 1.98e-6, 0.00e-6},
1798  {{0, 0, 0, 0, 3, 0, 0, 0}, -1.72e-6, 0.00e-6},
1799  {{0, 1, 0, 0, 1, 0, 0, 0}, -1.41e-6, -0.01e-6},
1800  {{0, 1, 0, 0, -1, 0, 0, 0}, -1.26e-6, -0.01e-6},
1801  // 11-20
1802  {{1, 0, 0, 0, -1, 0, 0, 0}, -0.63e-6, 0.00e-6},
1803  {{1, 0, 0, 0, 1, 0, 0, 0}, -0.63e-6, 0.00e-6},
1804  {{0, 1, 2, -2, 3, 0, 0, 0}, 0.46e-6, 0.00e-6},
1805  {{0, 1, 2, -2, 1, 0, 0, 0}, 0.45e-6, 0.00e-6},
1806  {{0, 0, 4, -4, 4, 0, 0, 0}, 0.36e-6, 0.00e-6},
1807  {{0, 0, 1, -1, 1, -8, 12, 0}, -0.24e-6, -0.12e-6},
1808  {{0, 0, 2, 0, 0, 0, 0, 0}, 0.32e-6, 0.00e-6},
1809  {{0, 0, 2, 0, 2, 0, 0, 0}, 0.28e-6, 0.00e-6},
1810  {{1, 0, 2, 0, 3, 0, 0, 0}, 0.27e-6, 0.00e-6},
1811  {{1, 0, 2, 0, 1, 0, 0, 0}, 0.26e-6, 0.00e-6},
1812  // 21-30
1813  {{0, 0, 2, -2, 0, 0, 0, 0}, -0.21e-6, 0.00e-6},
1814  {{0, 1, -2, 2, -3, 0, 0, 0}, 0.19e-6, 0.00e-6},
1815  {{0, 1, -2, 2, -1, 0, 0, 0}, 0.18e-6, 0.00e-6},
1816  {{0, 0, 0, 0, 0, 8, -13, -1}, -0.10e-6, 0.05e-6},
1817  {{0, 0, 0, 2, 0, 0, 0, 0}, 0.15e-6, 0.00e-6},
1818  {{2, 0, -2, 0, -1, 0, 0, 0}, -0.14e-6, 0.00e-6},
1819  {{1, 0, 0, -2, 1, 0, 0, 0}, 0.14e-6, 0.00e-6},
1820  {{0, 1, 2, -2, 2, 0, 0, 0}, -0.14e-6, 0.00e-6},
1821  {{1, 0, 0, -2, -1, 0, 0, 0}, 0.14e-6, 0.00e-6},
1822  {{0, 0, 4, -2, 4, 0, 0, 0}, 0.13e-6, 0.00e-6},
1823  // 31-33
1824  {{0, 0, 2, -2, 4, 0, 0, 0}, -0.11e-6, 0.00e-6},
1825  {{1, 0, -2, 0, -3, 0, 0, 0}, 0.11e-6, 0.00e-6},
1826  {{1, 0, -2, 0, -1, 0, 0, 0}, 0.11e-6, 0.00e-6}};
1827 
1828  // T^1 - do this manually
1829  // static const Coeffs C1[] = {
1830  // {{ 0, 0, 0, 0, 1, 0, 0, 0, 0 }, -0.87e-6, +0.00e-6 }
1831  //};
1832 
1833  // number of terms in each order
1834  const int n0 = int(sizeof Czero / sizeof(Coeffs));
1835  // const int n1 = int(sizeof C1 / sizeof(Coeffs));
1836 
1837  try
1838  {
1839  // coordinate transformation time
1840  double T(coordTransTime(t));
1841  // fundamental arguments l lp f d o lme lv le pa
1842  double farg[N];
1843 
1844  LOG(DEBUG7) << "\nT = " << fixed << setprecision(15) << T;
1845  farg[0] = L(T); // mean anomaly of the moon
1846  LOG(DEBUG7) << "L(T) = " << fixed << setprecision(15) << farg[0];
1847  farg[1] = Lp(T); // mean anomaly of the sun
1848  LOG(DEBUG7) << "Lp(T) = " << fixed << setprecision(15) << farg[1];
1849  farg[2] = F(T); // mean longitude of the moon - Omega
1850  LOG(DEBUG7) << "F(T) = " << fixed << setprecision(15) << farg[2];
1851  farg[3] = D(T); // mean elongation of the moon from the sun
1852  LOG(DEBUG7) << "D(T) = " << fixed << setprecision(15) << farg[3];
1853  farg[4] = Omega2003(T); // mean longitude of lunar ascending node
1854  LOG(DEBUG7) << "Omega2003(T) = " << fixed << setprecision(15)
1855  << farg[4];
1856  farg[5] = LV(T); // mean longitude of Venus
1857  LOG(DEBUG7) << "LV(T) = " << fixed << setprecision(15) << farg[5];
1858  farg[6] = LE(T); // mean longitude of Earth
1859  LOG(DEBUG7) << "LE(T) = " << fixed << setprecision(15) << farg[6];
1860  farg[7] = Pa(T); // general precession in longitude
1861  LOG(DEBUG7) << "Pa(T) = " << fixed << setprecision(15) << farg[7];
1862 
1863  // do the sums
1864  double ee(0.0);
1865  for (int i = n0 - 1; i >= 0; --i)
1866  { // order 0
1867  double arg(0.0);
1868  for (int j = 0; j < N; ++j)
1869  {
1870  if (Czero[i].coeff[j])
1871  {
1872  arg += Czero[i].coeff[j] * farg[j];
1873  }
1874  }
1875  ee += Czero[i].sincoeff * ::sin(arg);
1876  if (Czero[i].coscoeff)
1877  {
1878  ee += Czero[i].coscoeff * ::cos(arg);
1879  }
1880  }
1881 
1882  // the T^1 term
1883  ee += -0.87e-6 * ::sin(farg[4]) * T;
1884 
1885  // convert to radians
1886  ee *= ARCSEC_TO_RAD;
1887 
1888  return ee;
1889  }
1890  catch (Exception& e)
1891  {
1892  GNSSTK_RETHROW(e);
1893  }
1894  }
1895 
1896  //------------------------------------------------------------------------------
1897  /* Zonal tide terms for corrections of UT1mUTC when that quantity does not
1898  include tides (e.g. NGA EOP), ref. IERS 1996 Ch. 8, table 8.1 pg 74.
1899  param T, the coordinate transformation time at the time of interest
1900  param UT1mUT1R, the correction to UT1mUTC (seconds)
1901  param dlodR, the correction to the length of day (seconds)
1902  param domegaR, the correction to the Earth rotation rate (radians/second) */
1903  void EarthOrientation::UT1mUTCTidalCorrections(double T, double& UT1mUT1R,
1904  double& dlodR,
1905  double& domegaR)
1906  {
1907  // Define (all doubles) and all in radians
1908  double o(Omega(T)); // mean longitude of lunar ascending node
1909  double f(F(T)); // mean longitude of the moon - Omega
1910  double d(D(T)); // mean elongation of the moon from the sun
1911  double l(L(T)); // mean anomaly of the moon
1912  double lp(Lp(T)); // mean anomaly of the sun
1913 
1914  //-----------------------------------------------------------------------
1915  // include code that forms UT1mUT1R dlodR domegaR
1916  #include "IERS1996UT1mUTCData.hpp"
1917  }
1918 
1919  //---------------------------------------------------------------------------------
1920  /* Compute eps, the obliquity of the ecliptic, in radians,
1921  given T, the coordTransTime at the time of interest. IAU76 IAU80 for
1922  IERS1996,03 */
1923  double EarthOrientation::obliquity1996(double T)
1924  {
1925  /* double eps; // seconds of arc
1926  eps = T*(-46.8150 + T*(-0.00059 + T*0.001813));
1927  eps /= 3600.0; // convert to degrees
1928  eps += 23.4392911111111111; // = 84381.448/3600.0
1929  eps *= DEG_TO_RAD;
1930  return eps; */
1931  return (84381.448 + T * (-46.8150 + T * (-0.00059 + T * (0.001813)))) *
1932  ARCSEC_TO_RAD;
1933  }
1934 
1935  //---------------------------------------------------------------------------------
1936  /* Compute eps, the obliquity of the ecliptic, in radians,
1937  given T, the coordTransTime at the time of interest. */
1938  double EarthOrientation::obliquity2010(double T)
1939  {
1940  /* double eps;
1941  mean obliquity cf. sofa obl06.c
1942  eps = (84381.406 + (-46.836769 + (-0.0001831 + (0.00200340 +
1943  (-0.000000576 - 0.0000000434*T)*T)*T)*T)*T)*ARCSEC_TO_RAD;
1944  return eps; */
1945  return (84381.406 +
1946  T * (-46.836769 +
1947  T * (-0.0001831 +
1948  T * (0.00200340 +
1949  T * (-0.000000576 + T * (-0.0000000434)))))) *
1950  ARCSEC_TO_RAD;
1951  }
1952 
1953  //------------------------------------------------------------------------------
1954  /* Compute Greenwich Mean Sidereal Time, or the Greenwich hour angle of
1955  the mean vernal equinox (radians), given the UT1 time of interest,
1956  and UT1-UTC (sec), from the IERS bulletin. For IERS1996. cf sofa gmst82.c
1957  param t EphTime epoch of the rotation.
1958  param UT1mUTC, UT1-UTC in seconds, as found in the IERS bulletin.
1959  return GMST in radians */
1960  double EarthOrientation::GMST1996(EphTime t, double UT1mUTC, bool reduced)
1961  {
1962  try
1963  {
1964  // convert to UTC first
1965  t.convertSystemTo(TimeSystem::UTC);
1966 
1967  // if reduced, compute tidal terms
1968  if (reduced)
1969  {
1970  double dlodR, domegaR, UT1mUT1R;
1971  UT1mUTCTidalCorrections(coordTransTime(t), UT1mUT1R, dlodR,
1972  domegaR);
1973  UT1mUTC = UT1mUT1R - UT1mUTC;
1974  }
1975  // convert to UT1
1976  t += UT1mUTC;
1977 
1978  // dont use coordTransTime() b/c UT1 is needed here, not TT
1979  double T((t.dMJD() - JulianEpoch) / 36525.0);
1980 
1981  // Compute GMST in radians
1982  double G = -19089.45159 // first term is 24110.54841-43200. seconds
1983  + T * (8640184.812866 + T * (0.093104 - T * 6.2e-6));
1984  // convert seconds to days
1985  G /= 86400.0;
1986 
1987  /* this should be the same ... seconds (24060s = 6h 41min)
1988  G = 24110.54841 + (8640184.812866 + (0.093104 - 6.2e-6*T)*T)*T; sec
1989  G /= 86400.0; // instead, divide the numbers above manually
1990  G = 0.279057273264 + 100.0021390378009*T // seconds/86400 =
1991  days
1992  + (0.093104 - 6.2e-6*T)*T*T/86400.0;
1993  G += (UT1mUTC + t.secOfDay())/86400.0; // days */
1994 
1995  // add fraction of day
1996  double frac(0.5 + t.secOfDay() / 86400); // fraction of day
1997  if (frac > 1.0)
1998  {
1999  frac -= 1.0;
2000  }
2001  G += frac;
2002 
2003  // convert to radians
2004  G *= TWOPI; // radians
2005  G = ::fmod(G, TWOPI);
2006  if (G < 0.0)
2007  {
2008  G += TWOPI;
2009  }
2010 
2011  return G;
2012  }
2013  catch (Exception& e)
2014  {
2015  GNSSTK_RETHROW(e);
2016  }
2017  }
2018 
2019  //------------------------------------------------------------------------------
2020  /* Compute Greenwich Mean Sidereal Time, or the Greenwich hour angle of
2021  the mean vernal equinox (radians), given the coordinate time of interest,
2022  and UT1-UTC (sec), which comes from the IERS bulletin. For IERS2003
2023  param t EphTime epoch of the rotation.
2024  param UT1mUTC, UT1-UTC in seconds, as found in the IERS bulletin. */
2025  double EarthOrientation::GMST2003(EphTime t, double UT1mUTC)
2026  {
2027  try
2028  {
2029  // TT days since epoch
2030  double T(coordTransTime(t)), G;
2031  double era = EarthRotationAngle(t, UT1mUTC);
2032  G = ::fmod(
2033  era + (0.014506 +
2034  (4612.15739966 +
2035  (1.39667721 + (-0.00009344 + 0.00001882 * T) * T) * T) *
2036  T) *
2037  ARCSEC_TO_RAD,
2038  TWOPI);
2039 
2040  return G;
2041  }
2042  catch (Exception& e)
2043  {
2044  GNSSTK_RETHROW(e);
2045  }
2046  }
2047 
2048  //------------------------------------------------------------------------------
2049  /* Compute Greenwich Mean Sidereal Time, or the Greenwich hour angle of
2050  the mean vernal equinox (radians), given the coordinate time of interest,
2051  and UT1-UTC (sec), which comes from the IERS bulletin. For IERS2010
2052  param t EphTime epoch of the rotation.
2053  param UT1mUTC, UT1-UTC in seconds, as found in the IERS bulletin. */
2054  double EarthOrientation::GMST2010(EphTime t, double UT1mUTC)
2055  {
2056  try
2057  {
2058  double era = EarthRotationAngle(t, UT1mUTC); // radians
2059 
2060  // IERS2010 - cf sofa gmst06.c and TN36 eqn 5.32
2061  double T = coordTransTime(t);
2062  /* double G = ::fmod(era+(0.014506 + (4612.156534 + (1.3915817 +
2063  (-0.00000044
2064  + (-0.000029956 + (-0.0000000368)*T)*T)*T)*T)*T)*ARCSEC_TO_RAD,
2065  TWOPI);
2066  return G; */
2067  return (era +
2068  (0.014506 +
2069  T * (4612.156534 // NB era in radians already
2070  + T * (1.3915817 +
2071  T * (-0.00000044 + T * (-0.000029956 +
2072  T * (-0.0000000368)))))) *
2073  ARCSEC_TO_RAD);
2074  }
2075  catch (Exception& e)
2076  {
2077  GNSSTK_RETHROW(e);
2078  }
2079  }
2080 
2081  //---------------------------------------------------------------------------------
2082  /* Helper to compute the Greenwich hour angle of the true vernal equinox, or
2083  Greenwich Apparent Sidereal Time (GAST) in radians, for IERS1996,
2084  given the (UT) time of interest t, and, where T = coordTransTime(t),
2085  o = Omega(T) = mean longitude of lunar ascending node, in radians,
2086  eps = the obliquity of the ecliptic, in radians,
2087  dpsi = nutation in longitude (counted in the ecliptic),
2088  in seconds of arc. */
2089  double EarthOrientation::gast1996(EphTime t, double om, double eps,
2090  double dpsi, double UT1mUTC)
2091  {
2092  try
2093  {
2094  double G = GMST1996(t, UT1mUTC, false);
2095 
2096  // add equation of equinoxes: dpsi, eps and Omega terms
2097  double ee =
2098  dpsi * ::cos(eps) +
2099  (0.00264 * ::sin(om) + 0.000063 * ::sin(2.0 * om)) * ARCSEC_TO_RAD;
2100 
2101  LOG(DEBUG7) << "\nequequinox = " << fixed << setprecision(15)
2102  << showpos << 1.e3 * ee * RAD_TO_DEG << " / 1.e3";
2103  LOG(DEBUG7) << "\nGMST = " << fixed << setprecision(15) << showpos
2104  << G * RAD_TO_DEG;
2105 
2106  return (G + ee);
2107  }
2108  catch (Exception& e)
2109  {
2110  GNSSTK_RETHROW(e);
2111  }
2112  }
2113 
2114  //------------------------------------------------------------------------------
2115  /* Compute Greenwich Apparent Sidereal Time, or the Greenwich hour angle of
2116  the true vernal equinox (radians), given the coordinate time of interest,
2117  and UT1-UTC (sec), which comes from the IERS bulletin.
2118  param t EphTime epoch of the rotation.
2119  param UT1mUTC, UT1-UTC in seconds, as found in the IERS bulletin.
2120 
2121  GAST = Greenwich hour angle of the true vernal equinox
2122  GAST = GMST + dpsi*cos(eps) + 0.00264" * sin(Omega) +0.000063" *
2123  sin(2*Omega)
2124  (these terms account for the accumulated precession and nutation in
2125  right ascension and minimize any discontinuity in UT1)
2126 
2127  GMST = Greenwich hour angle of the mean vernal equinox
2128  = Greenwich Mean Sideral Time
2129  = GMST0 + r*[UTC + (UT1-UTC)]
2130  r = is the ratio of universal to sidereal time
2131  = 1.002737909350795 + 5.9006E-11*T' - 5.9e-15*T'^2
2132  T' = days'/36525
2133  days'= number of days elapsed since the Julian Epoch t0 (J2000)
2134  = +/-(integer+0.5)
2135  and
2136  (UT1-UTC) (seconds) is taken from the IERS bulletin
2137 
2138  GMST0 = GMST at 0h UT1
2139  = 6h 41min (50.54841+8640184.812866*T'+0.093104*T'^2-6.2E-6*T'^3)sec
2140 
2141  see pg 21 of the Reference (IERS 1996). */
2142  double EarthOrientation::GAST1996(EphTime t, double UT1mUTC, bool reduced)
2143  {
2144  try
2145  {
2146  double T(coordTransTime(t));
2147  double omega, deps, dpsi, G;
2148  double eps(obliquity1996(T)), epsa;
2149 
2150  nutationAngles1996(T, deps, dpsi, omega); // deps is not used...
2151 
2152  // if reduced (NGA), correct for tides
2153  double UT1mUT1R, dlodR, domegaR;
2154  if (reduced)
2155  {
2156  UT1mUTCTidalCorrections(T, UT1mUT1R, dlodR, domegaR);
2157  G = gast1996(t, omega, eps, dpsi, UT1mUT1R - UT1mUTC);
2158  }
2159  else
2160  {
2161  G = gast1996(t, omega, eps, dpsi, UT1mUTC);
2162  }
2163 
2164  return G;
2165  }
2166  catch (Exception& e)
2167  {
2168  GNSSTK_RETHROW(e);
2169  }
2170  }
2171 
2172  //------------------------------------------------------------------------------
2173  /* Compute Greenwich Apparent Sidereal Time, or the Greenwich hour angle of
2174  the true vernal equinox (radians), given the coordinate time of interest,
2175  and UT1-UTC (sec), which comes from the IERS bulletin.
2176  param t EphTime epoch of the rotation.
2177  param UT1mUTC, UT1-UTC in seconds, as found in the IERS bulletin. */
2178  double EarthOrientation::GAST2003(EphTime t, double UT1mUTC)
2179  {
2180  try
2181  {
2182  double omega, deps, dpsi, G;
2183  double T(coordTransTime(t));
2184 
2185  // precession and obliquity corrections (rad/century)
2186  double dpsipr, depspr;
2187  precessionRateCorrections2003(T, dpsipr, depspr); // dpsipr not used
2188 
2189  // compute mean obliquity from IERS Tech Note 32 Chapter 5, eqn 32.
2190  double eps(obliquity1996(T));
2191  // mean obliquity consistent with IAU 2000 P-N models
2192  double epsa(eps + depspr);
2193 
2194  nutationAngles2003(T, deps, dpsi);
2195 
2196  // Equation of the equinoxes.
2197  double ee = equationOfEquinoxes2003(t);
2198  LOG(DEBUG7) << "\nee = " << fixed << setprecision(15) << showpos
2199  << 1.e6 * ee * RAD_TO_DEG << " / 1.e6";
2200  ee = dpsi * ::cos(epsa) + ee;
2201 
2202  G = GMST2003(t, UT1mUTC) + ee;
2203 
2204  return G;
2205  }
2206  catch (Exception& e)
2207  {
2208  GNSSTK_RETHROW(e);
2209  }
2210  }
2211 
2212  //------------------------------------------------------------------------------
2213  /* Compute Greenwich Apparent Sidereal Time, or the Greenwich hour angle of
2214  the true vernal equinox (radians), given the coordinate time of interest,
2215  and UT1-UTC (sec), which comes from the IERS bulletin.
2216  param t EphTime epoch of the rotation.
2217  param UT1mUTC, UT1-UTC in seconds, as found in the IERS bulletin. */
2218  double EarthOrientation::GAST2010(EphTime t, double UT1mUTC)
2219  {
2220  try
2221  {
2222  Matrix<double> NutPreBias =
2223  preciseEarthRotation2010(coordTransTime(t));
2224 
2225  // extract X and Y coords of the CIP from the matrix
2226  // cf. sofa bpn2xy.c
2227  double X(NutPreBias(2, 0)), Y(NutPreBias(2, 1));
2228 
2229  // get T and the CIO locator s
2230  double T(coordTransTime(t));
2231  double s(S(T, X, Y, IERSConvention::IERS2010));
2232 
2233  // get ERA(UT1)
2234  double era = EarthRotationAngle(t, UT1mUTC);
2235 
2236  // equation of the origins. cf. sofa eors.c
2237  double eo, ax, xs, ys, zs, p, q;
2238  ax = X / (1.0 + NutPreBias(2, 2));
2239  xs = 1.0 - ax * X;
2240  ys = -ax * Y;
2241  zs = -X;
2242  p = NutPreBias(0, 0) * xs + NutPreBias(0, 1) * ys +
2243  NutPreBias(0, 2) * zs;
2244  q = NutPreBias(1, 0) * xs + NutPreBias(1, 1) * ys +
2245  NutPreBias(1, 2) * zs;
2246  eo = ((p != 0) || (q != 0)) ? s - ::atan2(q, p) : s;
2247 
2248  /* double G(era-eo);
2249  G = ::fmod(G, TWOPI);
2250  return G; */
2251  return ::fmod(era - eo, TWOPI);
2252  }
2253  catch (Exception& e)
2254  {
2255  GNSSTK_RETHROW(e);
2256  }
2257  }
2258 
2259  //---------------------------------------------------------------------------------
2260  /* Generate transformation matrix (3X3 rotation) due to the polar motion
2261  angles xp and yp (arcseconds), as found in the IERS bulletin;
2262  see class EarthOrientation.
2263  param xp, Earth wobble in arcseconds, as found in the IERS bulletin.
2264  param yp, Earth wobble in arcseconds, as found in the IERS bulletin.
2265  return Matrix<double>(3,3) rotation matrix */
2266  Matrix<double> EarthOrientation::polarMotionMatrix1996(double xp,
2267  double yp)
2268  {
2269  xp *= ARCSEC_TO_RAD;
2270  yp *= ARCSEC_TO_RAD;
2271  Matrix<double> R1, R2;
2272  R1 = rotation(-yp, 1);
2273  R2 = rotation(-xp, 2);
2274  return (R2 * R1);
2275  }
2276 
2277  //---------------------------------------------------------------------------------
2278  /* Generate transformation matrix (3X3 rotation) due to the polar motion
2279  angles xp and yp (arcseconds), as found in the IERS bulletin;
2280  The returned matrix R transforms the CIP into TRS: V(TRS) = R * V(CIP).
2281  Also see class EarthOrientation, sofa pom00. Also valid for IERS2010.
2282  param t EphTime epoch of the rotation.
2283  param xp, Earth wobble in arcseconds, as found in the IERS bulletin.
2284  param yp, Earth wobble in arcseconds, as found in the IERS bulletin.
2285  return Matrix<double>(3,3) rotation matrix CIP -> TRS */
2286  Matrix<double> EarthOrientation::polarMotionMatrix2003(EphTime t, double xp,
2287  double yp)
2288  {
2289  double sp(Sprime(t)); // parameter s' provides position of TEO on CIP
2290  LOG(DEBUG7) << "\nsp = " << fixed << setprecision(15) << setw(18) << sp
2291  << " with T = " << coordTransTime(t);
2292  xp *= ARCSEC_TO_RAD;
2293  yp *= ARCSEC_TO_RAD;
2294  Matrix<double> R1, R2, R3;
2295  R3 = rotation(sp, 3);
2296  R2 = rotation(-xp, 2);
2297  R1 = rotation(-yp, 1);
2298  return (R1 * R2 * R3);
2299  }
2300 
2301  //------------------------------------------------------------------------------
2302  /* Compute Fukushima-Williams angles for computing nutation, frame bias and
2303  precession matrices in IERS2010; cf. FukushimaWilliams().
2304  NB. fourth angle is obliquity. */
2305  void EarthOrientation::FukushimaWilliams(double T, double& gamb,
2306  double& phib, double& psib,
2307  double& eps)
2308  {
2309  // P03 bias+precession angles cf. sofa pfw06.c
2310  gamb = (-0.052928 +
2311  (10.556378 +
2312  (0.4932044 +
2313  (-0.00031238 + (-0.000002788 + 0.0000000260 * T) * T) * T) *
2314  T) *
2315  T) *
2316  ARCSEC_TO_RAD;
2317  phib = (84381.412819 +
2318  (-46.811016 +
2319  (0.0511268 +
2320  (0.00053289 + (-0.000000440 - 0.0000000176 * T) * T) * T) *
2321  T) *
2322  T) *
2323  ARCSEC_TO_RAD;
2324  psib = (-0.041775 +
2325  (5038.481484 +
2326  (1.5584175 +
2327  (-0.00018522 + (-0.000026452 - 0.0000000148 * T) * T) * T) *
2328  T) *
2329  T) *
2330  ARCSEC_TO_RAD;
2331 
2332  // obliquity
2333  eps = obliquity2010(T);
2334  }
2335 
2336  //---------------------------------------------------------------------------------
2337  /* Generate any of B = frame bias matrix
2338  PB = Precession*Bias matrix
2339  NPB = Nutation*Precession*Bias matrix
2340  for IERS 2010, given the four F-W angles with caviats, as follows.
2341  Get B by passing the full F-W angles at J2000 (T=0).
2342  Get PB by passing the full F-W angles at time time of interest.
2343  Get NPB by passing the full F-W angles at time time of interest with
2344  nutation angle corrections (nutationAngles2010()). Specifically,
2345  FukushimaWilliams(T, gamb, phib, psib, eps);
2346  nutationAngles2010(T, deps, dpsi);
2347  NPB = FukushimaWilliams(gamb, phib, psib+dpsi, eps+deps);
2348  Thus the precession matrix is computed as PB * transpose(B), and
2349  the nutation matrix is computes as N = NPB * transpose(PB).
2350  param gamb F-W angle
2351  param phib F-W angle
2352  param psib F-W angle
2353  param epsa F-W angle, the obliquity
2354  return 3x3 rotation matrix */
2355  Matrix<double> EarthOrientation::FukushimaWilliams(double gamb, double phib,
2356  double psib,
2357  double epsa)
2358  {
2359  Matrix<double> R;
2360  R = rotation(-epsa, 1) * rotation(-psib, 3) * rotation(phib, 1) *
2361  rotation(gamb, 3);
2362  return R;
2363  }
2364 
2365  //---------------------------------------------------------------------------------
2366  /* Nutation of the obliquity (deps) and of the longitude (dpsi), IERS 1996
2367  model (ref SOFA nut80.c - not found in Tech Note 21!).
2368  param T, the coordinate transformation time at the time of interest
2369  param deps, nutation of the obliquity, in radians (output)
2370  param dpsi, nutation of the longitude, in radians (output)
2371  param om, longitude mean ascending node of lunar orbit, from mean equinox */
2372  void EarthOrientation::nutationAngles1996(double T, double& deps,
2373  double& dpsi, double& om)
2374  {
2375 // include coefficients array; defines coeff[Ncoeff]
2376 #include "IERS1996NutationData.hpp"
2377 
2378  // Define fundamental arguments in radians - these do not appear elsewhere
2379  // (!?) Mean longitude of Moon minus mean longitude of Moon's perigee
2380  double el =
2381  ::fmod((485866.733 + (715922.633 + (31.310 + 0.064 * T) * T) * T) *
2382  ARCSEC_TO_RAD +
2383  ::fmod(1325.0 * T, 1.0) * TWOPI,
2384  TWOPI);
2385 
2386  // Mean longitude of Sun minus mean longitude of Sun's perigee
2387  double elp =
2388  ::fmod((1287099.804 + (1292581.224 + (-0.577 - 0.012 * T) * T) * T) *
2389  ARCSEC_TO_RAD +
2390  ::fmod(99.0 * T, 1.0) * TWOPI,
2391  TWOPI);
2392 
2393  // Mean longitude of Moon minus mean longitude of Moon's node
2394  double f =
2395  ::fmod((335778.877 + (295263.137 + (-13.257 + 0.011 * T) * T) * T) *
2396  ARCSEC_TO_RAD +
2397  ::fmod(1342.0 * T, 1.0) * TWOPI,
2398  TWOPI);
2399 
2400  // Mean elongation of Moon from Sun
2401  double d =
2402  ::fmod((1072261.307 + (1105601.328 + (-6.891 + 0.019 * T) * T) * T) *
2403  ARCSEC_TO_RAD +
2404  ::fmod(1236.0 * T, 1.0) * TWOPI,
2405  TWOPI);
2406 
2407  // Longitude of the mean ascending node of the lunar orbit on the
2408  // ecliptic, measured from the mean equinox of date
2409  om = ::fmod((450160.280 + (-482890.539 + (7.455 + 0.008 * T) * T) * T) *
2410  ARCSEC_TO_RAD +
2411  ::fmod(-5.0 * T, 1.0) * TWOPI,
2412  TWOPI);
2413 
2414  // sum the series
2415  double arg, scoeff, ccoeff;
2416  deps = dpsi = 0.0;
2417  for (int i = Ncoeff - 1; i >= 0; i--)
2418  {
2419  // form argument
2420  double arg = coeff[i].nl * el + coeff[i].nlp * elp + coeff[i].nf * f +
2421  coeff[i].nd * d + coeff[i].nom * om;
2422 
2423  // sine and cosine terms
2424  scoeff = coeff[i].sp + coeff[i].spt * T;
2425  ccoeff = coeff[i].ce + coeff[i].cet * T;
2426  if (scoeff != 0.0)
2427  {
2428  dpsi += scoeff * ::sin(arg);
2429  }
2430  if (ccoeff != 0.0)
2431  {
2432  deps += ccoeff * ::cos(arg);
2433  }
2434  }
2435 
2436  // convert from 0.1 milliarcseconds to radians
2437  deps *= ARCSEC_TO_RAD * 1.e-4;
2438  dpsi *= ARCSEC_TO_RAD * 1.e-4;
2439  }
2440 
2441  //---------------------------------------------------------------------------------
2442  /* Nutation of the obliquity (deps) and of the longitude (dpsi), IERS 2003 or
2443  IAU 2000A model (MHB2000 luni-solar and planetary nutation without free
2444  core n). param T, the coordinate transformation time at the time of
2445  interest param deps, nutation of the obliquity, in radians (output) param
2446  dpsi, nutation of the longitude, in radians (output) */
2447  void EarthOrientation::nutationAngles2003(double T, double& deps,
2448  double& dpsi)
2449  {
2450  // sin and cos coefficients have units 0.1 microarcsec = 1e-7as
2451  const double COEFF_TO_RAD(ARCSEC_TO_RAD * 1.0e-7);
2452 
2453 // -----------------------------------------
2454 // include huge static arrays of coefficients
2455 #include "IERS2003NutationData.hpp"
2456 
2457  /* -----------------------------------------
2458  Lunar-Solar nutation
2459  fundamental arguments, in radians */
2460  double l(L(T)); // mean anomaly of the moon
2461 
2462  double lp(
2463  ::fmod(1287104.79305 // mean anomaly of the sun MHB2000 value
2464  + T * (129596581.0481 +
2465  T * (-0.5532 + T * (0.000136 + T * (-0.00001149)))),
2466  ARCSEC_PER_CIRCLE) *
2467  ARCSEC_TO_RAD);
2468 
2469  double f(
2470  ::fmod(335779.526232 // mean longitude of moon minus Omega MHB2000
2471  + T * (1739527262.8478 +
2472  T * (-12.7512 + T * (-0.001037 + T * (0.00000417)))),
2473  ARCSEC_PER_CIRCLE) *
2474  ARCSEC_TO_RAD);
2475 
2476  double d(
2477  ::fmod(1072260.70369 // mean elongation moon from sun MHB2000
2478  + T * (1602961601.2090 +
2479  T * (-6.3706 + T * (0.006593 + T * (-0.00003169)))),
2480  ARCSEC_PER_CIRCLE) *
2481  ARCSEC_TO_RAD);
2482 
2483  double Om(Omega2003(T)); // mean longitude of lunar ascending node
2484 
2485  // initialize
2486  deps = dpsi = 0.0;
2487  int i;
2488  double arg, sina, cosa;
2489 
2490  // form the LS series
2491  for (i = NLS - 1; i >= 0; --i)
2492  {
2493  // argument
2494  arg =
2495  ::fmod(LSCoeff[i].nl * l + LSCoeff[i].nlp * lp + LSCoeff[i].nf * f +
2496  LSCoeff[i].nd * d + LSCoeff[i].nom * Om,
2497  TWOPI);
2498  sina = ::sin(arg);
2499  cosa = ::cos(arg);
2500  // term
2501  deps +=
2502  (LSCoeff[i].ce + LSCoeff[i].cet * T) * cosa + LSCoeff[i].se * sina;
2503  dpsi +=
2504  (LSCoeff[i].sp + LSCoeff[i].spt * T) * sina + LSCoeff[i].cp * cosa;
2505  }
2506 
2507  // -----------------------------------------
2508  // Planetary nutation
2509  // fundamental arguments, in radians.
2510 
2511  // NB MHB2000 values are very close to IERS2003; follow SOFA here TD ??
2512  // mean anomaly of the moon MHB2000 value
2513  l = ::fmod(2.35555598 + 8328.6914269554 * T, TWOPI);
2514  // mean longitude of the moon minus Omega MHB2000 value
2515  f = ::fmod(1.627905234 + 8433.466158131 * T, TWOPI);
2516  // mean elongation of the Moon from the Sun MHB2000 value
2517  d = ::fmod(5.198466741 + 7771.3771468121 * T, TWOPI);
2518  // mean longitude of lunar ascending node MHB2000 value
2519  Om = ::fmod(2.18243920 - 33.757045 * T, TWOPI);
2520 
2521  // mean longitude Mercury
2522  double lme(LMe(T));
2523  // mean longitude of Venus
2524  double lve(LV(T));
2525  // mean longitude of Earth
2526  double lea(LE(T));
2527  // mean longitude Mars
2528  double lma(LMa(T));
2529  // mean longitude Jupiter
2530  double lju(LJ(T));
2531  // mean longitude Saturn
2532  double lsa(LS(T));
2533  // mean longitude Uranus
2534  double lur(LU(T));
2535  // mean longitude Neptune
2536  double lne(::fmod(5.321159000 + 3.8127774000 * T, TWOPI));
2537  // general precession in longitude
2538  double pa(Pa(T));
2539 
2540  // form the planetary series
2541  for (i = NP - 1; i >= 0; --i)
2542  {
2543  // argument
2544  arg = ::fmod(PCoeff[i].nl * l + PCoeff[i].nf * f + PCoeff[i].nd * d +
2545  PCoeff[i].nom * Om + PCoeff[i].nme * lme +
2546  PCoeff[i].nve * lve + PCoeff[i].nea * lea +
2547  PCoeff[i].nma * lma + PCoeff[i].nju * lju +
2548  PCoeff[i].nsa * lsa + PCoeff[i].nur * lur +
2549  PCoeff[i].nne * lne + PCoeff[i].npa * pa,
2550  TWOPI);
2551  sina = ::sin(arg);
2552  cosa = ::cos(arg);
2553  // term
2554  deps += PCoeff[i].ce * cosa + PCoeff[i].se * sina;
2555  dpsi += PCoeff[i].sp * sina + PCoeff[i].cp * cosa;
2556  }
2557 
2558  // convert 0.1microarcsec to radians
2559  deps *= COEFF_TO_RAD;
2560  dpsi *= COEFF_TO_RAD;
2561 
2562  return;
2563  }
2564 
2565  //---------------------------------------------------------------------------------
2566  /* Nutation of the obliquity (deps) and of the longitude (dpsi), IERS 2010 or
2567  IAU 2000A model (MHB2000 luni-solar and planetary nutation without free
2568  core n) with P03 adjustments. cf. sofa nut06a.c param T, the coordinate
2569  transformation time at the time of interest param deps, nutation of the
2570  obliquity, in radians (output) param dpsi, nutation of the longitude, in
2571  radians (output) */
2572  void EarthOrientation::nutationAngles2010(double T, double& deps,
2573  double& dpsi)
2574  {
2575  nutationAngles2003(T, deps, dpsi);
2576  double fj2(-2.7774e-6 * T);
2577  dpsi *= (1.0 + 0.4697e-6 + fj2);
2578  deps *= (1.0 + fj2);
2579  }
2580 
2581  //---------------------------------------------------------------------------------
2582  /* Compute the nutation matrix, given
2583  eps, the obliquity of the ecliptic, in radians,
2584  dpsi, the nutation in longitude (counted in the ecliptic), in radians
2585  deps, the nutation in obliquity, in radians. */
2586  Matrix<double> EarthOrientation::nutationMatrix(double eps, double dpsi,
2587  double deps)
2588  {
2589  Matrix<double> R1 = rotation(eps, 1);
2590  Matrix<double> R2 = rotation(-dpsi, 3);
2591  Matrix<double> R3 = rotation(-(eps + deps), 1);
2592  return (R3 * R2 * R1);
2593  }
2594 
2595  //------------------------------------------------------------------------------
2596  /* IERS1996 nutation matrix, a 3x3 rotation matrix, given
2597  param T, the coordinate transformation time at the time of interest
2598  return nutation matrix Matrix<double>(3,3) */
2599  Matrix<double> EarthOrientation::nutationMatrix1996(double T)
2600  {
2601  double eps(obliquity1996(T)), deps, dpsi, om;
2602 
2603  nutationAngles1996(T, deps, dpsi, om);
2604  return nutationMatrix(eps, dpsi, deps);
2605  }
2606 
2607  //------------------------------------------------------------------------------
2608  /* IERS2003 nutation matrix, a 3x3 rotation matrix
2609  (including the frame bias matrix), given
2610  param T, the coordinate transformation time at the time of interest
2611  return nutation matrix Matrix<double>(3,3) */
2612  Matrix<double> EarthOrientation::nutationMatrix2003(double T)
2613  {
2614  double eps(obliquity1996(T)), deps, dpsi; // same as obliquity2003
2615  nutationAngles2003(T, deps, dpsi);
2616 
2617  /* Precession rate contributions with respect to IAU 2000
2618  Precession and obliquity corrections (radians) */
2619  double depspr = -0.02524 * ARCSEC_TO_RAD * T;
2620  eps += depspr;
2621 
2622  return nutationMatrix(eps, dpsi, deps);
2623  }
2624 
2625  //------------------------------------------------------------------------------
2626  /* IERS2010 nutation matrix, a 3x3 rotation matrix, given
2627  param T, the coordinate transformation time at the time of interest;
2628  cf. FukushimaWilliams().
2629  return nutation matrix Matrix<double>(3,3) */
2630  Matrix<double> EarthOrientation::nutationMatrix2010(double T)
2631  {
2632  double deps, dpsi, eps;
2633 
2634  /* get the F-W angles at epoch
2635  double gamb,phib,psib;
2636  FukushimaWilliams(T, gamb, phib, psib, eps);
2637 
2638  get nutation angles
2639  nutationAngles2010(T,deps,dpsi);
2640 
2641  construct nutation x precession x frame bias matrix
2642  NB this is the same as preciseEarthRotation2010(T)
2643  Matrix<double> NPB = FukushimaWilliams(gamb,phib,psib+dpsi,eps+deps);
2644 
2645  now get PB alone - see FukushimaWilliams()
2646  Matrix<double> PB = FukushimaWilliams(gamb,phib,psib,eps);
2647 
2648  return (NPB * transpose(PB)); */
2649 
2650  // same result
2651  nutationAngles2010(T, deps, dpsi);
2652  return nutationMatrix(obliquity2010(T), dpsi, deps);
2653  }
2654 
2655  //------------------------------------------------------------------------------
2656  /* Compute the IERS1996 precession matrix, a 3x3 rotation matrix, given
2657  param T, the coordinate transformation time at the time of interest
2658  return precession matrix Matrix<double>(3,3) */
2659  Matrix<double> EarthOrientation::precessionMatrix1996(double T)
2660  {
2661  /* IAU76 - ref McCarthy - seconds of arc
2662  NB t0==0 in sofa prec76.c - TD why do they do these things? */
2663  double TAR(T * ARCSEC_TO_RAD); // convert to radians
2664  double zeta = TAR * (2306.2181 + T * (0.30188 + T * 0.017998));
2665  double theta = TAR * (2004.3109 - T * (0.42665 + T * 0.041833));
2666  double z = TAR * (2306.2181 + T * (1.09468 + T * 0.018203));
2667 
2668  Matrix<double> R1 = rotation(-zeta, 3);
2669  Matrix<double> R2 = rotation(theta, 2);
2670  Matrix<double> R3 = rotation(-z, 3);
2671  Matrix<double> P = R3 * R2 * R1;
2672 
2673  return P;
2674  }
2675 
2676  //------------------------------------------------------------------------------
2677  /* Compute the IERS2003 precession matrix, a 3x3 rotation matrix, given
2678  param T, the coordinate transformation time at the time of interest
2679  Includes the frame bias matrix. cf sofa bp00.c
2680  return precession matrix Matrix<double>(3,3) */
2681  Matrix<double> EarthOrientation::precessionMatrix2003(double T)
2682  {
2683  // obliquity at the J2000.0 epoch
2684  static const double eps0(84381.448 * ARCSEC_TO_RAD);
2685  LOG(DEBUG7) << "\nobliquity at J2000 eps0 = " << fixed << setprecision(15)
2686  << showpos << eps0;
2687 
2688  // frame bias corrections in longitude and obliquity
2689  static const double psibias = -0.041775 * ARCSEC_TO_RAD;
2690  static const double epsbias = -0.0068192 * ARCSEC_TO_RAD;
2691  // ICRS right ascension of the J2000.0 equinox
2692  static const double raeps0 = -0.0146 * ARCSEC_TO_RAD;
2693  LOG(DEBUG7) << "frame bias psi = " << fixed << setprecision(15) << showpos
2694  << psibias;
2695  LOG(DEBUG7) << "frame bias eps = " << fixed << setprecision(15) << showpos
2696  << epsbias;
2697  LOG(DEBUG7) << "frame bias dra = " << fixed << setprecision(15) << showpos
2698  << raeps0;
2699 
2700  // precession angles
2701  double psia((5038.7784 + (-1.07259 + (-0.001147) * T) * T) * T *
2702  ARCSEC_TO_RAD);
2703  double epsa(eps0 + ((0.05127 + (-0.007726) * T) * T) * T * ARCSEC_TO_RAD);
2704  double chia((10.5526 + (-2.38064 + (-0.001125) * T) * T) * T *
2705  ARCSEC_TO_RAD);
2706  LOG(DEBUG7) << "\nprecession angle psi = " << fixed << setprecision(15)
2707  << showpos << psia;
2708  LOG(DEBUG7) << "precession angle eps = " << fixed << setprecision(15)
2709  << showpos << epsa;
2710  LOG(DEBUG7) << "precession angle chi = " << fixed << setprecision(15)
2711  << showpos << chia;
2712 
2713  /* Precession rate contributions with respect to IAU 2000 p-n models
2714  Precession and obliquity corrections (radians) cf sofa pr00.c */
2715  double dpsipr, depspr;
2716  precessionRateCorrections2003(T, dpsipr, depspr);
2717 
2718  // Apply precession corrections
2719  LOG(DEBUG7) << "precession rate = " << fixed << setprecision(15)
2720  << showpos << dpsipr << " " << depspr;
2721  psia += dpsipr;
2722  epsa += depspr;
2723 
2724  // Frame bias matrix
2725  Matrix<double> R1 = rotation(raeps0, 3);
2726  Matrix<double> R2 = rotation(psibias * ::sin(eps0), 2);
2727  Matrix<double> R3 = rotation(-epsbias, 1);
2728  Matrix<double> FrameBias(R3 * R2 * R1);
2729  LOG(DEBUG7) << "\nframe bias matrix:\n"
2730  << fixed << setprecision(15) << showpos << FrameBias;
2731 
2732  // Precession matrix
2733  R1 = rotation(eps0, 1);
2734  R2 = rotation(-psia, 3);
2735  R3 = rotation(-epsa, 1);
2736  Matrix<double> R4 = rotation(chia, 3);
2737  Matrix<double> Precess(R4 * R3 * R2 * R1);
2738  LOG(DEBUG7) << "\nprecession matrix:\n"
2739  << fixed << setprecision(15) << showpos << Precess;
2740 
2741  LOG(DEBUG7) << "\nprecession*framebias matrix:\n"
2742  << fixed << setprecision(15) << showpos
2743  << (Precess * FrameBias);
2744 
2745  return (Precess * FrameBias);
2746  }
2747 
2748  //------------------------------------------------------------------------------
2749  /* Compute the IERS2003 precession and obliquity rate corrections, IAU 2000
2750  param T, the coordinate transformation time at the time of interest
2751  return precession, obliquity corrections in radians */
2752  void EarthOrientation::precessionRateCorrections2003(double T, double& dpsi,
2753  double& deps)
2754  {
2755  // Precession rate contributions with respect to IAU 2000
2756  // Precession and obliquity corrections (radians)
2757  dpsi = -0.29965 * ARCSEC_TO_RAD * T;
2758  deps = -0.02524 * ARCSEC_TO_RAD * T;
2759  }
2760 
2761  //------------------------------------------------------------------------------
2762  /* IERS2010 frame bias matrix, a 3x3 rotation matrix; cf.
2763  FukushimaWilliams(). return frame bias matrix Matrix<double>(3,3) */
2764  Matrix<double> EarthOrientation::biasMatrix2010()
2765  {
2766  // get F-W angles at J2000
2767  double gamb, phib, psib, epsa;
2768  FukushimaWilliams(0.0, gamb, phib, psib, epsa);
2769 
2770  // frame bias matrix
2771  return FukushimaWilliams(gamb, phib, psib, epsa);
2772  }
2773 
2774  //---------------------------------------------------------------------------------
2775  /* Compute the IERS2010 precession matrix, a 3x3 rotation matrix, given
2776  param T, the coordinate transformation time at the time of interest
2777  Does not include the frame bias matrix. Cf. FukushimaWilliams().
2778  return precession matrix Matrix<double>(3,3) */
2779  Matrix<double> EarthOrientation::precessionMatrix2010(double T)
2780  {
2781  // the F-W angles
2782  double gamb, phib, psib, epsa;
2783 
2784  // get frame bias matrix
2785  Matrix<double> B = biasMatrix2010();
2786 
2787  // get F-W angles at epoch
2788  FukushimaWilliams(T, gamb, phib, psib, epsa);
2789 
2790  // precession x frame bias matrix
2791  Matrix<double> PB = FukushimaWilliams(gamb, phib, psib, epsa);
2792 
2793  return (PB * transpose(B));
2794  }
2795 
2796  //------------------------------------------------------------------------------
2797  /* Generate precise transformation matrix (3X3 rotation) for Earth motion due
2798  to precession, nutation and frame bias (NPB matrix), at the given time of
2799  interest, for IERS 2003. param T coordTransTime of interest return 3x3
2800  rotation matrix */
2801  Matrix<double> EarthOrientation::preciseEarthRotation2003(double T)
2802  {
2803  try
2804  {
2805  Matrix<double> N = nutationMatrix2003(T);
2806  Matrix<double> P = precessionMatrix2003(T); // includes bias
2807 
2808  /* Matrix<double> NPB(N*P);
2809  LOG(DEBUG7) << "\nNPB matrix:\n" << fixed << setprecision(15) <<
2810  setw(18)
2811  << showpos << NPB; */
2812 
2813  return (N * P);
2814  }
2815  catch (Exception& e)
2816  {
2817  GNSSTK_RETHROW(e);
2818  }
2819  }
2820 
2821  //------------------------------------------------------------------------------
2822  /* Generate precise transformation matrix (3X3 rotation) for Earth motion due
2823  to precession, nutation and frame bias (NPB matrix), at the given time of
2824  interest, for IERS 2010. param T coordTransTime of interest return 3x3
2825  rotation matrix */
2826  Matrix<double> EarthOrientation::preciseEarthRotation2010(double T)
2827  {
2828  try
2829  {
2830  double deps, dpsi, epsa;
2831 
2832  // get the F-W angles
2833  double gamb, phib, psib;
2834  FukushimaWilliams(T, gamb, phib, psib, epsa);
2835 
2836  // get nutation angles
2837  nutationAngles2010(T, deps, dpsi);
2838 
2839  // construct nutation x precession x frame bias matrix
2840  return FukushimaWilliams(gamb, phib, psib + dpsi, epsa + deps);
2841  }
2842  catch (Exception& e)
2843  {
2844  GNSSTK_RETHROW(e);
2845  }
2846  }
2847 
2848  //---------------------------------------------------------------------------------
2849  /* Generate the full transformation matrix (3x3 rotation) relating the ECEF
2850  frame to the conventional inertial frame, using IERS 1996 conventions.
2851  Input is the time of interest, the polar motion angles xp and yp
2852  (arcsecs), and UT1-UTC (sec) (xp,yp and UT1-UTC are just as found in the
2853  IERS bulletin; see class EarthOrientation). */
2854  Matrix<double> EarthOrientation::ECEFtoInertial1996(EphTime t, double xp,
2855  double yp,
2856  double UT1mUTC,
2857  bool reduced)
2858  {
2859  try
2860  {
2861  Matrix<double> P, N, W, S;
2862 
2863  double T = coordTransTime(t);
2864 
2865  // precession
2866  P = precessionMatrix1996(T);
2867  LOG(DEBUG7) << "\nprecession matrix:\n"
2868  << fixed << setprecision(15) << setw(18) << showpos << P;
2869 
2870  // nutation
2871  double eps, deps, dpsi, om;
2872  // mean obliquity radians
2873  eps = obliquity1996(T);
2874  LOG(DEBUG7) << "\nmean obliquity " << fixed << setprecision(15)
2875  << showpos << eps;
2876  // nutation angles - om is used in gast
2877  nutationAngles1996(T, deps, dpsi, om);
2878  LOG(DEBUG7) << "\nnutation angles psi eps " << fixed
2879  << setprecision(15) << showpos << dpsi << " " << deps;
2880  // nutation matrix
2881  N = nutationMatrix(eps, dpsi, deps);
2882  LOG(DEBUG7) << "\nnutation matrix:\n"
2883  << fixed << setprecision(15) << setw(18) << showpos << N;
2884 
2885  LOG(DEBUG7) << "\nNPB matrix:\n"
2886  << fixed << setprecision(15) << setw(18) << showpos
2887  << N * P;
2888 
2889  // if reduced (NGA), correct UT1mUTC for tides
2890  double UT1mUT1R, dlodR, domegaR;
2891  if (reduced)
2892  {
2893  UT1mUTCTidalCorrections(T, UT1mUT1R, dlodR, domegaR);
2894  UT1mUTC = UT1mUT1R - UT1mUTC;
2895  }
2896 
2897  double g = gast1996(t, om, eps, dpsi, UT1mUTC);
2898  LOG(DEBUG7) << "\nGAST = " << fixed << setprecision(15) << showpos
2899  << g * RAD_TO_DEG;
2900 
2901  S = rotation(g, 3);
2902  LOG(DEBUG7) << "\ncelestial-to-terrestrial matrix (no polar motion):\n"
2903  << fixed << setprecision(15) << setw(18) << showpos
2904  << S * N * P;
2905 
2906  // Polar Motion
2907  W = polarMotionMatrix1996(xp, yp);
2908  LOG(DEBUG7) << "\npolar motion matrix:\n"
2909  << fixed << setprecision(15) << setw(18) << showpos << W;
2910 
2911  return transpose(W * S * N * P);
2912  }
2913  catch (Exception& e)
2914  {
2915  GNSSTK_RETHROW(e);
2916  }
2917  }
2918 
2919  //---------------------------------------------------------------------------------
2920  /* Generate the full transformation matrix (3x3 rotation) relating the ECEF
2921  frame to the conventional inertial frame, using IERS 2003 conventions.
2922  Input is the time of interest, the polar motion angles xp and yp
2923  (arcsecs), and UT1-UTC (sec) (xp,yp and UT1-UTC are just as found in the
2924  IERS bulletin; see class EarthOrientation). */
2925  Matrix<double> EarthOrientation::ECEFtoInertial2003(EphTime t, double xp,
2926  double yp,
2927  double UT1mUTC)
2928  {
2929  try
2930  {
2931  Matrix<double> P, N, R, W;
2932 
2933  double T(coordTransTime(t));
2934 
2935  if (LOGlevel >= DEBUG7)
2936  {
2937  double gmst = GMST2003(t, UT1mUTC);
2938  LOG(DEBUG7) << "\nGMST = " << fixed << setprecision(15) << showpos
2939  << gmst * RAD_TO_DEG;
2940 
2941  double gast = GAST2003(t, UT1mUTC);
2942  LOG(DEBUG7) << "\nGAST = " << fixed << setprecision(15) << showpos
2943  << gast * RAD_TO_DEG;
2944  }
2945 
2946  // nutation
2947  double deps, dpsi, dpsipr, depspr;
2948  nutationAngles2003(T, deps, dpsi);
2949  LOG(DEBUG7) << "\nnutation angles psi eps " << fixed
2950  << setprecision(15) << showpos << dpsi << " " << deps;
2951 
2952  /* Precession rate contributions with respect to IAU 2000
2953  Precession and obliquity corrections (radians) */
2954  precessionRateCorrections2003(T, dpsipr, depspr);
2955  LOG(DEBUG7) << "\nprecession-rate " << fixed << setprecision(15)
2956  << showpos << dpsipr << " " << depspr;
2957 
2958  double eps(obliquity1996(T)); // same as 2003
2959  LOG(DEBUG7) << "\nmean obliquity " << fixed << setprecision(15)
2960  << showpos << eps;
2961  eps += depspr;
2962 
2963  N = nutationMatrix(eps, dpsi, deps);
2964  LOG(DEBUG7) << "\nnutation matrix:\n"
2965  << fixed << setprecision(15) << setw(18) << showpos << N;
2966 
2967  // precession
2968  P = precessionMatrix2003(T);
2969 
2970  Matrix<double> NPB(N * P);
2971  LOG(DEBUG7) << "\nNPB matrix:\n"
2972  << fixed << setprecision(15) << setw(18) << showpos << NPB;
2973 
2974  // ERA replaces GAST in the Earth rotation matrix
2975  double era(EarthRotationAngle(t, UT1mUTC));
2976  LOG(DEBUG7) << "\nERA = " << fixed << setprecision(15) << showpos
2977  << era * RAD_TO_DEG;
2978  R = rotation(era, 3);
2979 
2980  /* double gast = GAST2003(t, UT1mUTC);
2981  R = rotation(gast,3); */
2982 
2983  // polar motion
2984  W = polarMotionMatrix2003(t, xp, yp);
2985  LOG(DEBUG7) << "\npolar motion matrix:\n"
2986  << fixed << setprecision(15) << setw(18) << showpos << W;
2987 
2988  return transpose(W * R * N * P);
2989  }
2990  catch (Exception& e)
2991  {
2992  GNSSTK_RETHROW(e);
2993  }
2994  }
2995 
2996  //---------------------------------------------------------------------------------
2997  /* Generate the full transformation matrix (3x3 rotation) relating the ECEF
2998  frame to the conventional inertial frame, using IERS 2010 conventions.
2999  Input is the time of interest, the polar motion angles xp and yp
3000  (arcsecs), and UT1-UTC (sec) (xp,yp and UT1-UTC are just as found in the
3001  IERS bulletin; see class EarthOrientation). */
3002  Matrix<double> EarthOrientation::ECEFtoInertial2010(EphTime t, double xp,
3003  double yp,
3004  double UT1mUTC)
3005  {
3006  try
3007  {
3008  double T(coordTransTime(t));
3009 
3010  /* get the CIO coordinates and s
3011  note that X,Y could also be obtained as (2,0),(2,1) components
3012  of FukushimaWilliams() */
3013  double X, Y, s;
3014  XYCIO(T, X, Y);
3015  s = S(T, X, Y, IERSConvention::IERS2010);
3016  LOG(DEBUG7) << "X = " << fixed << setprecision(15) << showpos << X;
3017  LOG(DEBUG7) << "Y = " << fixed << setprecision(15) << showpos << Y;
3018  LOG(DEBUG7) << "s\" = " << fixed << setprecision(15)
3019  << s / ARCSEC_TO_RAD;
3020 
3021  /* compute transformation GCRS-to-CIRS or
3022  inertial-to-intermediate-celestial cf. sofa c2ixys */
3023  double r2(X * X + Y * Y); // squared radius
3024  double e(r2 != 0.0 ? ::atan2(Y, X) : 0.0); // spherical angles
3025  double d(::atan(::sqrt(r2 / (1.0 - r2)))); //
3026  Matrix<double> GCRStoCIRS;
3027  GCRStoCIRS = rotation(-(e + s), 3) * rotation(d, 2) * rotation(e, 3);
3028  LOG(DEBUG7) << "\nNPB matrix:\n"
3029  << fixed << setprecision(15) << setw(18) << showpos
3030  << GCRStoCIRS;
3031 
3032  // note that we could have called preciseEarthRotation2010() instead
3033 
3034  // get ERA at UT1
3035  double era = EarthRotationAngle(t, UT1mUTC);
3036  LOG(DEBUG7) << "\nERA = " << fixed << setprecision(15) << showpos
3037  << era * RAD_TO_DEG;
3038 
3039  // compute transf. CIRS-to-TIRS or
3040  // intermediate-celestial-to-terrestrial
3041  Matrix<double> CIRStoTIRS;
3042  CIRStoTIRS = rotation(era, 3);
3043  LOG(DEBUG7) << "\ncelestial-to-terrestrial matrix (no polar motion):\n"
3044  << fixed << setprecision(15) << setw(18) << showpos
3045  << CIRStoTIRS * GCRStoCIRS;
3046 
3047  // compute the polar motion matrix, TIRS-to-ITRS
3048  // double sprime(Sprime(T));
3049  Matrix<double> polarMotion(
3050  polarMotionMatrix2003(t, xp, yp)); // 2010 == 2003
3051  LOG(DEBUG7) << "\npolar motion matrix:\n"
3052  << fixed << setprecision(15) << setw(18) << showpos
3053  << polarMotion;
3054 
3055  // combine to get GCRS-to-ITRS
3056  Matrix<double> GCRStoITRS;
3057  GCRStoITRS = polarMotion * CIRStoTIRS * GCRStoCIRS;
3058 
3059  // invert to get ITRS-to-GCRS or ECEFtoInertial
3060  return (transpose(GCRStoITRS));
3061  }
3062  catch (Exception& e)
3063  {
3064  GNSSTK_RETHROW(e);
3065  }
3066  }
3067 
3068 } // end namespace gnsstk
3069 
3070 /*=============================================================================
3071 ** SOME of these routines, as noted, are based on, but not simply copied from,
3072 *SOFA;
3073 ** SOFA has the following license.
3074 **
3075 ** Copyright (C) 2012
3076 ** Standards Of Fundamental Astronomy Board
3077 ** of the International Astronomical Union.
3078 **
3079 ** =====================
3080 ** SOFA Software License
3081 ** =====================
3082 **
3083 ** NOTICE TO USER:
3084 **
3085 ** BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING SIX TERMS AND
3086 ** CONDITIONS WHICH APPLY TO ITS USE.
3087 **
3088 ** 1. The Software is owned by the IAU SOFA Board ("SOFA").
3089 **
3090 ** 2. Permission is granted to anyone to use the SOFA software for any
3091 ** purpose, including commercial applications, free of charge and
3092 ** without payment of royalties, subject to the conditions and
3093 ** restrictions listed below.
3094 **
3095 ** 3. You (the user) may copy and distribute SOFA source code to others,
3096 ** and use and adapt its code and algorithms in your own software,
3097 ** on a world-wide, royalty-free basis. That portion of your
3098 ** distribution that does not consist of intact and unchanged copies
3099 ** of SOFA source code files is a "derived work" that must comply
3100 ** with the following requirements:
3101 **
3102 ** a) Your work shall be marked or carry a statement that it
3103 ** (i) uses routines and computations derived by you from
3104 ** software provided by SOFA under license to you; and
3105 ** (ii) does not itself constitute software provided by and/or
3106 ** endorsed by SOFA.
3107 **
3108 ** b) The source code of your derived work must contain descriptions
3109 ** of how the derived work is based upon, contains and/or differs
3110 ** from the original SOFA software.
3111 **
3112 ** c) The names of all routines in your derived work shall not
3113 ** include the prefix "iau" or "sofa" or trivial modifications
3114 ** thereof such as changes of case.
3115 **
3116 ** d) The origin of the SOFA components of your derived work must
3117 ** not be misrepresented; you must not claim that you wrote the
3118 ** original software, nor file a patent application for SOFA
3119 ** software or algorithms embedded in the SOFA software.
3120 **
3121 ** e) These requirements must be reproduced intact in any source
3122 ** distribution and shall apply to anyone to whom you have
3123 ** granted a further right to modify the source code of your
3124 ** derived work.
3125 **
3126 ** Note that, as originally distributed, the SOFA software is
3127 ** intended to be a definitive implementation of the IAU standards,
3128 ** and consequently third-party modifications are discouraged. All
3129 ** variations, no matter how minor, must be explicitly marked as
3130 ** such, as explained above.
3131 **
3132 ** 4. You shall not cause the SOFA software to be brought into
3133 ** disrepute, either by misuse, or use for inappropriate tasks, or
3134 ** by inappropriate modification.
3135 **
3136 ** 5. The SOFA software is provided "as is" and SOFA makes no warranty
3137 ** as to its use or performance. SOFA does not and cannot warrant
3138 ** the performance or results which the user may obtain by using the
3139 ** SOFA software. SOFA makes no warranties, express or implied, as
3140 ** to non-infringement of third party rights, merchantability, or
3141 ** fitness for any particular purpose. In no event will SOFA be
3142 ** liable to the user for any consequential, incidental, or special
3143 ** damages, including any lost profits or lost savings, even if a
3144 ** SOFA representative has been advised of such damages, or for any
3145 ** claim by any third party.
3146 **
3147 ** 6. The provision of any version of the SOFA software under the terms
3148 ** and conditions specified herein does not imply that future
3149 ** versions will also be made available under the same terms and
3150 ** conditions.
3151 *
3152 ** In any published work or commercial product which uses the SOFA
3153 ** software directly, acknowledgement (see www.iausofa.org) is
3154 ** appreciated.
3155 **
3156 ** Correspondence concerning SOFA software should be addressed as
3157 ** follows:
3158 **
3159 ** By email: sofa@ukho.gov.uk
3160 ** By post: IAU SOFA Center
3161 ** HM Nautical Almanac Office
3162 ** UK Hydrographic Office
3163 ** Admiralty Way, Taunton
3164 ** Somerset, TA1 2DN
3165 ** United Kingdom
3166 **
3167 **=============================================================================*/
Ncoeff
const int Ncoeff
Number of terms in the series.
Definition: IERS1996NutationData.hpp:169
IERS2010CIOSeriesData.hpp
gnsstk::EphTime::convertSystemTo
void convertSystemTo(const TimeSystem &ts)
Definition: EphTime.hpp:96
se
double se
obliquity cos, T*cos, sin coefficients
Definition: IERS2003NutationData.hpp:48
amp
static const double amp[]
Amplitudes (microarcsec); indexed using the iamp array below.
Definition: IERS2010CIOSeriesData.hpp:1630
nf
int nf
Definition: IERS1996NutationData.hpp:44
domegaR
domegaR
Definition: IERS1996UT1mUTCData.hpp:55
example6.mjd
mjd
Definition: example6.py:102
coeff
static const struct @1 coeff[]
Constants used in the nutation models, adapted from SOFA nut80.c.
npa
int npa
coefficient of general precession
Definition: IERS2003NutationData.hpp:810
NLS
const int NLS
Number of terms in the luni-solar nutation model.
Definition: IERS2003NutationData.hpp:801
NP
const int NP
Number of terms in the planetary nutation model.
Definition: IERS2003NutationData.hpp:1573
nFAlunarsolar
static const int nFAlunarsolar[][5]
multipliers of the fundamental arguments, lunar and solar terms
Definition: IERS2010CIOSeriesData.hpp:56
gnsstk::computeFundamentalArgs
static void computeFundamentalArgs(double T, double args[6])
Definition: EarthOrientation.cpp:274
gnsstk::correctEarthRotationZonalTides2003
static void correctEarthRotationZonalTides2003(const double args[6], double &dUT, double &dld, double &dom)
Definition: EarthOrientation.cpp:629
docstring_generator.args
args
Definition: docstring_generator.py:106
gnsstk::EphTime::setTimeSystem
void setTimeSystem(TimeSystem sys)
Definition: EphTime.hpp:141
nme
int nme
Definition: IERS2003NutationData.hpp:809
nju
int nju
Definition: IERS2003NutationData.hpp:809
gnsstk::rotation
Matrix< T > rotation(T angle, int axis)
Definition: MatrixOperators.hpp:495
nlp
int nlp
Definition: IERS1996NutationData.hpp:44
gnsstk::TrackingCode::Y
@ Y
Encrypted legacy GPS precise code.
IERS1996NutationData.hpp
logstream.hpp
gnsstk::correctEarthRotationZonalTides
static void correctEarthRotationZonalTides(const double args[6], double &dUT, double &dld, double &dom)
Definition: EarthOrientation.cpp:516
gnsstk::PI
const double PI
GPS value of PI; also specified by GAL.
Definition: GNSSconstants.hpp:62
gnsstk::GMST
static double GMST(const CommonTime &t)
Definition: SolarPosition.cpp:58
XYcoeff
static const double XYcoeff[2][MAXPT+1]
polynominal coeff for X, Y in arcseconds
Definition: IERS2010CIOSeriesData.hpp:51
gnsstk
For Sinex::InputHistory.
Definition: BasicFramework.cpp:50
jasc
static const int jasc[]
for use with amp[] : sin or cos
Definition: IERS2010CIOSeriesData.hpp:2107
LOGlevel
#define LOGlevel
Definition: logstream.hpp:323
std::sin
double sin(gnsstk::Angle x)
Definition: Angle.hpp:144
gnsstk::EarthOrientation::yp
double yp
Polar motion angle yp, in arcseconds.
Definition: EarthOrientation.hpp:100
IERS2003NutationData.hpp
MiscMath.hpp
nne
int nne
Definition: IERS2003NutationData.hpp:809
gnsstk::Exception
Definition: Exception.hpp:151
UT1mUT1R
UT1mUT1R
Definition: IERS1996UT1mUTCData.hpp:49
nur
int nur
Definition: IERS2003NutationData.hpp:809
gnsstk::RAD_TO_DEG
static const double RAD_TO_DEG
Conversion Factor from radians to degrees (units: degrees)
Definition: GNSSconstants.hpp:78
nsa
int nsa
Definition: IERS2003NutationData.hpp:809
NFALS
static const int NFALS
Number of frequencies: lunar solar fundamental args.
Definition: IERS2010CIOSeriesData.hpp:778
gnsstk::transpose
SparseMatrix< T > transpose(const SparseMatrix< T > &M)
transpose
Definition: SparseMatrix.hpp:829
gnsstk::Matrix< double >
example4.time
time
Definition: example4.py:103
example4.err
err
Definition: example4.py:126
nd
int nd
Definition: IERS1996NutationData.hpp:44
gnsstk::EphTime
Definition: EphTime.hpp:67
nl
int nl
Definition: IERS1996NutationData.hpp:44
nma
int nma
Definition: IERS2003NutationData.hpp:809
LSCoeff
static const struct @2 LSCoeff[]
arg
double arg
Definition: IERS1996UT1mUTCData.hpp:48
nFAplanetary
static const int nFAplanetary[][14]
multipliers of the fundamental arguments, planetary terms
Definition: IERS2010CIOSeriesData.hpp:781
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::TrackingCode::P
@ P
Legacy GPS precise code.
japt
static const int japt[]
for use with amp[] : power of T
Definition: IERS2010CIOSeriesData.hpp:2110
std::operator<<
std::ostream & operator<<(std::ostream &s, gnsstk::StringUtils::FFLead v)
Definition: FormattedDouble_T.cpp:44
cp
double cp
longitude sin, T*sin, cos coefficients
Definition: IERS2003NutationData.hpp:47
iamp
static const int iamp[]
indexes into the amplitudes array for each frequency; length = NFALS + NFAP
Definition: IERS2010CIOSeriesData.hpp:1509
LOG
#define LOG(level)
define the macro that is used to write to the log stream
Definition: logstream.hpp:315
nea
int nea
Definition: IERS2003NutationData.hpp:809
PCoeff
static const struct @3 PCoeff[]
Planetary coefficients, sin and cos coeff have units 0.1 microarcsec.
EarthOrientation.hpp
std
Definition: Angle.hpp:142
dlodR
dlodR
Definition: IERS1996UT1mUTCData.hpp:54
gnsstk::EarthOrientation::convention
IERSConvention convention
IERS convention appropriate for this data.
Definition: EarthOrientation.hpp:106
NFAP
static const int NFAP
Number of frequencies: planetary fundamental args.
Definition: IERS2010CIOSeriesData.hpp:1506
gnsstk::EphTime::secOfDay
double secOfDay() const
Definition: EphTime.hpp:177
gnsstk::correctEarthRotationLibrations
static void correctEarthRotationLibrations(const double args[6], double &dUT, double &dld)
Definition: EarthOrientation.cpp:742
gnsstk::correctEOPOceanTides
static void correctEOPOceanTides(double mjd, double &dxp, double &dyp, double &dUT)
Definition: EarthOrientation.cpp:327
nve
int nve
Definition: IERS2003NutationData.hpp:809
jaxy
static const int jaxy[]
for use with amp[] : X or Y
Definition: IERS2010CIOSeriesData.hpp:2104
GNSSTK_THROW
#define GNSSTK_THROW(exc)
Definition: Exception.hpp:366
nom
int nom
coefficients of l,l',F,D,Om
Definition: IERS1996NutationData.hpp:44
gnsstk::EphTime::dMJD
double dMJD() const
Definition: EphTime.hpp:171
gnsstk::DEBUG7
@ DEBUG7
Definition: logstream.hpp:58
gnsstk::EarthOrientation
Definition: EarthOrientation.hpp:93
gnsstk::LagrangeInterpolation
T LagrangeInterpolation(const std::vector< T > &X, const std::vector< T > &Y, const T &x, T &err)
Definition: MiscMath.hpp:98
gnsstk::EphTime::setMJD
void setMJD(long double mjd)
Definition: EphTime.hpp:155
MAXPT
static const int MAXPT
order of polynomials in T for X and Y
Definition: IERS2010CIOSeriesData.hpp:48
IERS1996UT1mUTCData.hpp
NAmp
static const int NAmp
Number of amplitude coefficients.
Definition: IERS2010CIOSeriesData.hpp:2101
gnsstk::EarthOrientation::xp
double xp
Polar motion angle xp, in arcseconds.
Definition: EarthOrientation.hpp:97
gnsstk::EarthOrientation::UT1mUTC
double UT1mUTC
Time offset UT1 minus UTC, in seconds.
Definition: EarthOrientation.hpp:103
sp
double sp
Definition: IERS1996NutationData.hpp:46
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:38