OceanLoadTides.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 
40 
41 //------------------------------------------------------------------------------------
42 // system includes
43 #include <algorithm>
44 #include <fstream>
45 #include <iostream>
46 // GNSSTk
47 #include "GNSSconstants.hpp"
48 #include "MiscMath.hpp"
49 #include "StringUtils.hpp"
50 // geomatics
51 #include "CubicSpline.hpp"
52 #include "OceanLoadTides.hpp"
53 #include "RobustStats.hpp" // for QSort
54 //#include "logstream.hpp" // TEMP
55 
56 using namespace std;
57 
58 namespace gnsstk
59 {
60 
61  using namespace StringUtils;
62 
63  // Number of standard (Schwiderski) tides read from BLQ file
64  const int OceanLoadTides::NSTD = 11;
65  // Number of derived tides computed by deriveTides()
66  const int OceanLoadTides::NDER = 342;
67 
68  //---------------------------------------------------------------------------------
69  /* Open and read the given file, containing ocean loading coefficients, and
70  initialize this object for the sites names in the input list that match a
71  name in the file (case sensitive). Return the number of successfully
72  initialized site names, and remove those sites from the input list.
73  Ocean loading files can be obtained from the web. For example all the ITRF
74  sites are found at
75  ftp://maia.usno.navy.mil/conventions/chapter7/olls25.blq Also, at
76  http://www.oso.chalmers.se/~loading one may submit site label and position
77  for one or more sites, and the resulting ocean loading file will be
78  emailed.
79  @param sites vector<string> On input contains site labels found in
80  the
81  file, on output contains only sites that were NOT found.
82  If sites is empty, all sites are loaded.
83  @param filename string Input ocean loading file name.
84  @return the number of sites successfully initialized.
85  @throw if the file could not be opened. */
86  int OceanLoadTides::initializeSites(vector<string>& sites, string filename)
87  {
88  try
89  {
90  bool allsites = false;
91  if (sites.size() == 0)
92  {
93  allsites = true; // return 0;
94  }
95  int i, n;
96 
97  ifstream infile(filename.c_str());
98  if (!infile || !infile.is_open())
99  {
100  Exception e("File " + filename + " could not be opened.");
101  GNSSTK_THROW(e);
102  }
103 
104  n = 0; // number of successes
105  bool looking = true; // true if looking for a site name
106  double lat, lon;
107  vector<double> coeff;
108  string site;
109  while (1)
110  { // read the file
111  int count;
112  string line, word;
113 
114  // get the next line
115  getline(infile, line);
116  stripTrailing(line, '\r');
117 
118  // process line
119  if (!line.empty())
120  {
121  word = firstWord(line);
122  // LOG(VERBOSE) << "Word is " << word << " and line is " << line;
123 
124  if (word == "$$")
125  { // NB ignore header - assume column order, etc.
126  // pick out the lat/lon
127  if (!looking)
128  {
129  while (!line.empty())
130  {
131  word = stripFirstWord(line);
132  if (word == string("lon/lat:"))
133  {
134  lon = asDouble(stripFirstWord(line));
135  lat = asDouble(stripFirstWord(line));
136  break;
137  }
138  }
139  }
140  }
141  /* should test be line length <= 21 ? ... what if site name =
142  number
143  else if(looking && !isDecimalString(word)) { */
144  else if (looking && line.length() <= 21)
145  {
146  // site name
147  site = line;
148  stripTrailing(site, string("\n"));
149  stripTrailing(site, string("\r"));
150  stripTrailing(site);
151  stripLeading(site);
152  // LOG(VERBOSE) << "Found site " << site;
153  if (allsites)
154  {
155  // LOG(VERBOSE) << "Push back " << site;
156  looking = false;
157  sites.push_back(site);
158  }
159  else
160  {
161  for (i = 0; i < sites.size(); i++)
162  {
163  // LOG(VERBOSE) << "Compare " << sites[i];
164  if (site == sites[i])
165  {
166  looking = false;
167  break;
168  }
169  }
170  }
171  if (!looking)
172  { // found a site
173  count = 0;
174  coeff.clear();
175  lat = lon = 0.0;
176  }
177  }
178  else if (!looking)
179  { // not comment and not looking - must be data
180  if (numWords(line) != 11)
181  {
182  Exception e("File " + filename +
183  " is corrupted for site " + site +
184  " - offending line follows\n" + line);
185  GNSSTK_THROW(e);
186  }
187  // LOG(VERBOSE) << "Push back line " << line;
188  for (i = 0; i < 11; i++)
189  coeff.push_back(asDouble(stripFirstWord(line)));
190  count++;
191  if (count == 6)
192  { // success
193  ostringstream oss;
194  oss << fixed;
195  for (i = 0; i < coeff.size(); i++)
196  {
197  if (i < 33)
198  {
199  oss << " " << setprecision(5) << setw(7) << coeff[i];
200  }
201  else
202  {
203  oss << " " << setprecision(1) << setw(7) << coeff[i];
204  }
205  if ((i + 1) % 11 == 0)
206  {
207  oss << "\n";
208  }
209  }
210  // LOG(VERBOSE) << " Found site " << site << " with
211  // coefficients:"; LOG(VERBOSE) << oss.str();
212 
213  // update coeff map
214  coefficientMap[site] = coeff;
215  n++;
216  // update position map
217  coeff.clear();
218  coeff.push_back(lat);
219  coeff.push_back(lon);
220  positionMap[site] = coeff;
221 
222  // erase a vector element
223  if (!allsites)
224  {
225  vector<string>::iterator pos;
226  pos = find(sites.begin(), sites.end(), site);
227  if (pos != sites.end())
228  {
229  sites.erase(pos);
230  }
231  }
232  looking = true;
233  }
234  }
235 
236  } // end if line not empty
237 
238  if (infile.eof() || !infile.good())
239  {
240  break;
241  }
242 
243  } // end loop over lines in the file
244 
245  return n;
246  }
247  catch (Exception& e)
248  {
249  GNSSTK_RETHROW(e);
250  }
251  catch (exception& e)
252  {
253  Exception E("std except: " + string(e.what()));
254  GNSSTK_THROW(E);
255  }
256  catch (...)
257  {
258  Exception e("Unknown exception");
259  GNSSTK_THROW(e);
260  }
261  }
262 
263  //---------------------------------------------------------------------------------
264  /* Compute the site displacement vector at the given time for the given site.
265  The site must have been successfully initialized; if not an exception is
266  thrown.
267  @param site string Input name of the site; must be the same as previously
268  successfully passed to initializeSites().
269  @param t EphTime Input time of interest.
270  @return Triple containing the North, East and Up components of the site
271  displacement in meters.
272  @throw if the site has not been initialized. */
273  Triple OceanLoadTides::computeDisplacement11(string site, EphTime time)
274  {
275  try
276  {
277  if (!isValid(site))
278  {
279  Exception e("Site " + site + " has not been initialized.");
280  GNSSTK_THROW(e);
281  }
282 
283  // get the coefficients for this site
284  vector<double> coeff = coefficientMap[site];
285 
286  // get the astronomical arguments in radians
287  double angles[11];
288  // inline this SchwiderskiArg(int(t.year())-1900, t.DOY(),
289  // t.secOfDay(), angles);
290  {
291  double fday(time.secOfDay());
292  long jday(
293  static_cast<long>(time.lMJD() + MJD_JDAY + fday / SEC_PER_DAY));
294  int iyear, imm, iday;
295  convertJDtoCalendar(jday, iyear, imm, iday);
296  iyear -= 1900;
297 
298  // ordering is: M2, S2, N2, K2, K1, O1, P1, Q1, Mf, Mm, Ssa
299  // which are : { semi-diurnal }{ diurnal }{long-period}
300  static const double speed[11] = {
301  1.40519E-4, 1.45444E-4, 1.37880E-4, 1.45842E-4,
302  0.72921E-4, 0.67598E-4, 0.72523E-4, 0.64959E-4,
303  0.053234E-4, 0.026392E-4, 0.003982E-4};
304  static const double angfac[44] = {
305  // sun
306  2.0, 0.0, 2.0, 2.0, // 4 : M2, S2, N2, K2
307  1.0, 1.0, -1.0, 1.0, // 8 : K1, O1, P1, Q1
308  0.0, 0.0, 2.0, // 11 : Mf, Mm, Ssa
309  // moon
310  -2.0, 0.0, -3.0, 0.0, // 15 : M2, S2, N2, K2
311  0.0, -2.0, 0.0, -3.0, // 19 : K1, O1, P1, Q1
312  2.0, 1.0, 0.0, // 22 : Mf, Mm, Ssa
313  // lunar perigee
314  0.0, 0.0, 1.0, 0.0, // 26 : M2, S2, N2, K2
315  0.0, 0.0, 0.0, 1.0, // 30 : K1, O1, P1, Q1
316  0.0, -1.0, 0.0, // 33 : Mf, Mm, Ssa
317  // two pi
318  0.0, 0.0, 0.0, 0.0, // 37 : M2, S2, N2, K2
319  0.25, -0.25, -0.25, -0.25, // 41 : K1, O1, P1, Q1
320  0.0, 0.0, 0.0 // 44 : Mf, Mm, Ssa
321  };
322 
323  int icapd = iday + 365 * (iyear - 75) + ((iyear - 73) / 4);
324 
325  // double capt = (27392.500528+1.000000035*double(icapd))/36525.0;
326  double capt =
327  0.74996579132101300 + 2.73785088295687885e-5 * double(icapd);
328 
329  // mean longitude of sun at beginning of day
330  double H0 = 279.69668 + (36000.768930485 + 0.000303 * capt) * capt;
331 
332  // mean longitude of moon at beginning of day
333  double S0 =
334  ((0.0000019 * capt - 0.001133) * capt + 481267.88314137) * capt +
335  270.434358;
336 
337  // mean longitude of lunar perigee at beginning of day
338  double P0 =
339  ((-0.000012 * capt - 0.010325) * capt + 4069.0340329577) * capt +
340  334.329653;
341 
342  // convert to radians
343  // static const double dtr = 0.0174532925199;
344  H0 *= DEG_TO_RAD;
345  S0 *= DEG_TO_RAD;
346  P0 *= DEG_TO_RAD;
347 
348  // LOG(INFO) << "Schwiderski " << iday << " " << fixed <<
349  // setprecision(5)
350  //<< setw(11) << fday << " " << icapd << " " << capt
351  //<< " " << H0 << " " << S0 << " " << P0;
352 
353  static const double twopi = 6.28318530718;
354  for (int k = 0; k < 11; k++)
355  {
356  angles[k] = speed[k] * fday + angfac[k] * H0 +
357  angfac[11 + k] * S0 + angfac[22 + k] * P0 +
358  angfac[33 + k] * twopi;
359  angles[k] = ::fmod(angles[k], twopi);
360  if (angles[k] < 0.0)
361  {
362  angles[k] += twopi;
363  }
364  }
365  } // end SchwiderskiArg()
366 
367  /* compute the radial, west and south components
368  coefficients are stored by rows: radial, west, south; first amp,
369  then phase column order same as in SchwiderskiArg() [ as in the file
370  ] */
371  Triple dc;
372  for (int i = 0; i < 3; i++)
373  { // components
374  dc[i] = 0.0;
375  for (int j = 0; j < 11; j++) // tidal modes
376  dc[i] += coeff[i * 11 + j] *
377  ::cos(angles[j] - coeff[33 + i * 11 + j] * DEG_TO_RAD);
378  }
379 
380  // convert radial,west,south to north,east,up
381  double temp = dc[0];
382  dc[0] = -dc[2]; // N = -S
383  dc[1] = -dc[1]; // E = -W
384  dc[2] = temp; // U = rad
385 
386  return dc;
387  }
388  catch (Exception& e)
389  {
390  GNSSTK_RETHROW(e);
391  }
392  catch (exception& e)
393  {
394  Exception E("std except: " + string(e.what()));
395  GNSSTK_THROW(E);
396  }
397  catch (...)
398  {
399  Exception e("Unknown exception");
400  GNSSTK_THROW(e);
401  }
402  }
403 
404  //---------------------------------------------------------------------------------
405  /* Compute the site displacement vector at the given time for the given site.
406  The site must have been successfully initialized; if not an exception is
407  thrown.
408  @param site string Input name of the site; must be the same as previously
409  successfully passed to initializeSites().
410  @param t EphTime Input time of interest.
411  @return Triple containing the North, East and Up components of the site
412  displacement in meters.
413  @throw if the site has not been initialized, if the time system is
414  unknown,
415  if there is corruption in the static arrays, or . */
416  Triple OceanLoadTides::computeDisplacement(string site, EphTime time)
417  {
418  try
419  {
420  ostringstream oss; // TEMP
421  int i;
422 
423  if (!isValid(site))
424  {
425  Exception e("Site " + site + " has not been initialized.");
426  GNSSTK_THROW(e);
427  }
428 
429  // get the coefficients for this site
430  vector<double> coeff = coefficientMap[site];
431 
432  // Cartwright-Tayler numbers of Scherneck tides
433  // ordering is: M2, S2, N2, K2, K1, O1, P1, Q1, Mf, Mm, Ssa
434 
435  // standard 11 Scherneck tides:
436  static const NVector SchInd[] = {
437  {2, 0, 0, 0, 0, 0}, // M2
438  {2, 2, -2, 0, 0, 0}, // S2
439  {2, -1, 0, 1, 0, 0}, // N2
440  {2, 2, 0, 0, 0, 0}, // K2
441  {1, 1, 0, 0, 0, 0}, // K1
442  {1, -1, 0, 0, 0, 0}, // O1
443  {1, 1, -2, 0, 0, 0}, // P1
444  {1, -2, 0, 1, 0, 0}, // Q1
445  {0, 2, 0, 0, 0, 0}, // Mf
446  {0, 1, 0, -1, 0, 0}, // Mm
447  {0, 0, 2, 0, 0, 0}, // Ssa
448  };
449 
450  // NB there must be 11 std tides in SchInd[]
451  if ((int)(sizeof(SchInd) / sizeof(NVector)) != NSTD)
452  {
453  Exception e("Static SchInd array is corrupted");
454  GNSSTK_THROW(e);
455  }
456 
457  // compute time argument
458  EphTime ttag(time);
459  ttag.convertSystemTo(TimeSystem::UTC);
460  double dayfr(ttag.secOfDay() / 86400.0);
461  ttag.convertSystemTo(TimeSystem::TT);
462  // T = EarthOrientation::CoordTransTime()
463  double T((ttag.dMJD() - 51544.5) / 36525.0);
464 
465  // get the Delauney arguments and frequencies at t
466  double Del[5], freqDel[5]; // degrees and cycles/day
467  Del[0] =
468  134.9634025100 + // EarthOrientation::L()
469  T * (477198.8675605000 +
470  T * (0.0088553333 + T * (0.0000143431 + T * (-0.0000000680))));
471  Del[1] =
472  357.5291091806 + // EarthOrientation::Lp()
473  T *
474  (35999.0502911389 +
475  T * (-0.0001536667 + T * (0.0000000378 + T * (-0.0000000032))));
476  Del[2] =
477  93.2720906200 + // EarthOrientation::F()
478  T *
479  (483202.0174577222 +
480  T * (-0.0035420000 + T * (-0.0000002881 + T * (0.0000000012))));
481  Del[3] =
482  297.8501954694 + // EarthOrientation::D()
483  T *
484  (445267.1114469445 +
485  T * (-0.0017696111 + T * (0.0000018314 + T * (-0.0000000088))));
486  Del[4] =
487  125.0445550100 + // EarthOrientation::Omega2003()
488  T * (-1934.1362619722 +
489  T * (0.0020756111 + T * (0.0000021394 + T * (-0.0000000165))));
490  for (i = 0; i < 5; i++)
491  Del[i] = ::fmod(Del[i], 360.0);
492  freqDel[0] = 0.0362916471 + 0.0000000013 * T;
493  freqDel[1] = 0.0027377786;
494  freqDel[2] = 0.0367481951 - 0.0000000005 * T;
495  freqDel[3] = 0.0338631920 - 0.0000000003 * T;
496  freqDel[4] = -0.0001470938 + 0.0000000003 * T;
497 
498  // convert to Doodson (Darwin) variables
499  double Dood[6], freqDood[6];
500  Dood[0] = 360.0 * dayfr - Del[3];
501  Dood[1] = Del[2] + Del[4];
502  Dood[2] = Dood[1] - Del[3];
503  Dood[3] = Dood[1] - Del[0];
504  Dood[4] = -Del[4];
505  Dood[5] = Dood[2] - Del[1];
506  for (i = 0; i < 6; i++)
507  Dood[i] = ::fmod(Dood[i], 360.0);
508 
509  freqDood[0] = 1.0 - freqDel[3];
510  freqDood[1] = freqDel[2] + freqDel[4];
511  freqDood[2] = freqDood[1] - freqDel[3];
512  freqDood[3] = freqDood[1] - freqDel[0];
513  freqDood[4] = -freqDel[4];
514  freqDood[5] = freqDood[2] - freqDel[1];
515 
516  // find amplitudes and phases for vertical, west and south components,
517  // for all 342 derived tides, from standard tides
518  double amp[NSTD], phs[NSTD];
519  double ampS[NDER], ampW[NDER],
520  ampU[NDER]; // south,west,up component amp.s
521  double phsS[NDER], phsW[NDER],
522  phsU[NDER]; // south,west,up component phs.s
523  double freq[NDER]; // frequencies (same for S,W,U)
524 
525  // vertical
526  int nder; // number returned, may be < NDER
527  for (i = 0; i < NSTD; i++)
528  {
529  amp[i] = coeff[i];
530  phs[i] = -coeff[33 + i];
531  }
532  // oss.str(""); oss << fixed << setprecision(5) << "TEST Amp 1 ";
533  // for(i=0; i<NSTD; i++) oss << " " << setw(8) << amp[i];
534  // LOG(INFO) << oss.str();
535  // oss.str(""); oss << fixed << setprecision(1) << "TEST Phs 1 ";
536  // for(i=0; i<NSTD; i++) oss << " " << setw(8) << phs[i];
537  // LOG(INFO) << oss.str();
538  // LOG(INFO) << "TEST T,DAYFR,DELTA" << fixed << setprecision(15)
539  // << setw(25) << T << setw(25) << dayfr << setw(25) <<
540  // ttag.secOfDay();
541  // LOG(INFO) << "TEST Delauneys " << fixed << setprecision(15)
542  // << setw(22) << Del[0] << setw(22) << Del[1] << setw(22) << Del[2]
543  // << setw(22) << Del[3] << setw(22) << Del[4];
544  // LOG(INFO) << "TEST Del freqs " << fixed << setprecision(15)
545  // << setw(22) << freqDel[0] << setw(22) << freqDel[1] << setw(22)
546  // << freqDel[2] << setw(22) << freqDel[3] << setw(22) << freqDel[4];
547  // LOG(INFO) << "TEST Doods " << fixed << setprecision(15)
548  // << setw(22) << Dood[0] << setw(22) << Dood[1] << setw(22)
549  // << Dood[2] << setw(22) << Dood[3] << setw(22) << Dood[4]
550  // << setw(22) << Dood[5];
551  // LOG(INFO) << "TEST Dood freqs" << fixed << setprecision(15)
552  // << setw(22) << freqDood[0] << setw(22) << freqDood[1] << setw(22)
553  // << freqDood[2] << setw(22) << freqDood[3] << setw(22) <<
554  // freqDood[4]
555  // << setw(22) << freqDood[5];
556  nder = deriveTides(SchInd, amp, phs, Dood, freqDood, ampU, phsU, freq,
557  NSTD);
558  // LOG(INFO) << "Vertical returned " << nder << " derived tides";
559 
560  // west
561  for (i = 0; i < NSTD; i++)
562  {
563  amp[i] = coeff[11 + i];
564  phs[i] = -coeff[44 + i];
565  }
566  // oss.str(""); oss << fixed << setprecision(5) << "TEST Amp 2 ";
567  // for(i=0; i<NSTD; i++) oss << " " << setw(8) << amp[i];
568  // LOG(INFO) << oss.str();
569  // oss.str(""); oss << fixed << setprecision(1) << "TEST Phs 2 ";
570  // for(i=0; i<NSTD; i++) oss << " " << setw(8) << phs[i];
571  // LOG(INFO) << oss.str();
572  nder = deriveTides(SchInd, amp, phs, Dood, freqDood, ampW, phsW, freq,
573  NSTD);
574  // LOG(INFO) << "West returned " << nder << " derived tides";
575 
576  // south
577  for (i = 0; i < NSTD; i++)
578  {
579  amp[i] = coeff[22 + i];
580  phs[i] = -coeff[55 + i];
581  }
582  // oss.str(""); oss << fixed << setprecision(5) << "TEST Amp 3 ";
583  // for(i=0; i<NSTD; i++) oss << " " << setw(8) << amp[i];
584  // LOG(INFO) << oss.str();
585  // oss.str(""); oss << fixed << setprecision(1) << "TEST Phs 3 ";
586  // for(i=0; i<NSTD; i++) oss << " " << setw(8) << phs[i];
587  // LOG(INFO) << oss.str();
588  nder = deriveTides(SchInd, amp, phs, Dood, freqDood, ampS, phsS, freq,
589  NSTD);
590  // LOG(INFO) << "TEST First 40 South amp, phase";
591  // for(i=0; i<40; i++) LOG(INFO) << "TEST " << setw(2) << i+1 << fixed
592  // << setprecision(15) << setw(22) << ampS[i] << setw(22) << phsS[i];
593 
594  // sum up
595  Triple dc(0.0, 0.0, 0.0); // U S W
596  for (i = 0; i < nder; i++)
597  {
598  dc[0] += ampU[i] * ::cos(phsU[i] * DEG_TO_RAD);
599  // LOG(INFO)<<"TEST LOOP U " << setw(3) << i+1 << fixed <<
600  // setprecision(15)
601  // << setw(22) << ampU[i]*::cos(phsU[i]*DEG_TO_RAD) << setw(22) <<
602  // dc[0];
603  }
604  // LOG(INFO) << "TEST RECURS result U " << fixed << setprecision(15)
605  // << setw(22) << dc[0];
606 
607  for (i = 0; i < nder; i++)
608  {
609  dc[1] += ampS[i] * ::cos(phsS[i] * DEG_TO_RAD);
610  // LOG(INFO)<<"TEST LOOP S " << setw(3) << i+1 << fixed <<
611  // setprecision(15)
612  // << setw(22) << ampS[i]*::cos(phsS[i]*DEG_TO_RAD) << setw(22) <<
613  // dc[1];
614  }
615  // LOG(INFO) << "TEST RECURS result S " << fixed << setprecision(15)
616  // << setw(22) << dc[1];
617 
618  for (i = 0; i < nder; i++)
619  {
620  dc[2] += ampW[i] * ::cos(phsW[i] * DEG_TO_RAD);
621  // LOG(INFO)<<"TEST LOOP W " << setw(3) << i+1 << fixed <<
622  // setprecision(15)
623  // << setw(22) << ampW[i]*::cos(phsW[i]*DEG_TO_RAD) << setw(22) <<
624  // dc[2];
625  }
626  // LOG(INFO) << "TEST RECURS result W " << fixed << setprecision(15)
627  // << setw(22) << dc[2];
628 
629  // LOG(INFO) << "TEST " << fixed << setprecision(6)
630  // << " " << dc[0] << " " << dc[1] << " " << dc[2];
631 
632  // convert vertical,south,west to north,east,up
633  double temp = dc[0];
634  dc[0] = -dc[1]; // N = -S
635  dc[1] = -dc[2]; // E = -W
636  dc[2] = temp; // U = U
637 
638  return dc;
639  }
640  catch (Exception& e)
641  {
642  GNSSTK_RETHROW(e);
643  }
644  catch (exception& e)
645  {
646  Exception E("std except: " + string(e.what()));
647  GNSSTK_THROW(E);
648  }
649  catch (...)
650  {
651  Exception e("Unknown exception");
652  GNSSTK_THROW(e);
653  }
654 
655  } // end Triple OceanLoadTides::computeDisplacement
656 
657  //---------------------------------------------------------------------------------
658  int OceanLoadTides::deriveTides(const NVector SchInd[], const double amp[],
659  const double phs[], const double Dood[],
660  const double freqDood[], double ampDer[],
661  double phsDer[], double freqDer[],
662  const int Nin)
663  {
664  // indexes for std tides: M2, S2, N2, K2, K1, O1, P1, Q1, Mf, Mm, Ssa
665  static const int stdindex[] = {0, 1, 2, 3, 109, 110,
666  111, 112, 263, 264, 265};
667 
668  static const double DerAmp[] = {
669  .632208, .294107, .121046, .079915, .023818, -.023589, .022994,
670  .019333, -.017871, .017192, .016018, .004671, -.004662, -.004519,
671  .004470, .004467, .002589, -.002455, -.002172, .001972, .001947,
672  .001914, -.001898, .001802, .001304, .001170, .001130, .001061,
673  -.001022, -.001017, .001014, .000901, -.000857, .000855, .000855,
674  .000772, .000741, .000741, -.000721, .000698, .000658, .000654,
675  -.000653, .000633, .000626, -.000598, .000590, .000544, .000479,
676  -.000464, .000413, -.000390, .000373, .000366, .000366, -.000360,
677  -.000355, .000354, .000329, .000328, .000319, .000302, .000279,
678  -.000274, -.000272, .000248, -.000225, .000224, -.000223, -.000216,
679  .000211, .000209, .000194, .000185, -.000174, -.000171, .000159,
680  .000131, .000127, .000120, .000118, .000117, .000108, .000107,
681  .000105, -.000102, .000102, .000099, -.000096, .000095, -.000089,
682  -.000085, -.000084, -.000081, -.000077, -.000072, -.000067, .000066,
683  .000064, .000063, .000063, .000063, .000062, .000062, -.000060,
684  .000056, .000053, .000051, .000050, .368645, -.262232, -.121995,
685  -.050208, .050031, -.049470, .020620, .020613, .011279, -.009530,
686  -.009469, -.008012, .007414, -.007300, .007227, -.007131, -.006644,
687  .005249, .004137, .004087, .003944, .003943, .003420, .003418,
688  .002885, .002884, .002160, -.001936, .001934, -.001798, .001690,
689  .001689, .001516, .001514, -.001511, .001383, .001372, .001371,
690  -.001253, -.001075, .001020, .000901, .000865, -.000794, .000788,
691  .000782, -.000747, -.000745, .000670, -.000603, -.000597, .000542,
692  .000542, -.000541, -.000469, -.000440, .000438, .000422, .000410,
693  -.000374, -.000365, .000345, .000335, -.000321, -.000319, .000307,
694  .000291, .000290, -.000289, .000286, .000275, .000271, .000263,
695  -.000245, .000225, .000225, .000221, -.000202, -.000200, -.000199,
696  .000192, .000183, .000183, .000183, -.000170, .000169, .000168,
697  .000162, .000149, -.000147, -.000141, .000138, .000136, .000136,
698  .000127, .000127, -.000126, -.000121, -.000121, .000117, -.000116,
699  -.000114, -.000114, -.000114, .000114, .000113, .000109, .000108,
700  .000106, -.000106, -.000106, .000105, .000104, -.000103, -.000100,
701  -.000100, -.000100, .000099, -.000098, .000093, .000093, .000090,
702  -.000088, .000083, -.000083, -.000082, -.000081, -.000079, -.000077,
703  -.000075, -.000075, -.000075, .000071, .000071, -.000071, .000068,
704  .000068, .000065, .000065, .000064, .000064, .000064, -.000064,
705  -.000060, .000056, .000056, .000053, .000053, .000053, -.000053,
706  .000053, .000053, .000052, .000050, -.066607, -.035184, -.030988,
707  .027929, -.027616, -.012753, -.006728, -.005837, -.005286, -.004921,
708  -.002884, -.002583, -.002422, .002310, .002283, -.002037, .001883,
709  -.001811, -.001687, -.001004, -.000925, -.000844, .000766, .000766,
710  -.000700, -.000495, -.000492, .000491, .000483, .000437, -.000416,
711  -.000384, .000374, -.000312, -.000288, -.000273, .000259, .000245,
712  -.000232, .000229, -.000216, .000206, -.000204, -.000202, .000200,
713  .000195, -.000190, .000187, .000180, -.000179, .000170, .000153,
714  -.000137, -.000119, -.000119, -.000112, -.000110, -.000110, .000107,
715  -.000095, -.000095, -.000091, -.000090, -.000081, -.000079, -.000079,
716  .000077, -.000073, .000069, -.000067, -.000066, .000065, .000064,
717  -.000062, .000060, .000059, -.000056, .000055, -.000051};
718 
719  static const NVector DerInd[] = {
720  {2, 0, 0, 0, 0, 0}, {2, 2, -2, 0, 0, 0},
721  {2, -1, 0, 1, 0, 0}, // M2,S2,N2
722  {2, 2, 0, 0, 0, 0}, {2, 2, 0, 0, 1, 0},
723  {2, 0, 0, 0, -1, 0}, // K2,x,x
724  {2, -1, 2, -1, 0, 0}, {2, -2, 2, 0, 0, 0},
725  {2, 1, 0, -1, 0, 0}, {2, 2, -3, 0, 0, 1},
726  {2, -2, 0, 2, 0, 0}, {2, -3, 2, 1, 0, 0},
727  {2, 1, -2, 1, 0, 0}, {2, -1, 0, 1, -1, 0},
728  {2, 3, 0, -1, 0, 0}, {2, 1, 0, 1, 0, 0},
729  {2, 2, 0, 0, 2, 0}, {2, 2, -1, 0, 0, -1},
730  {2, 0, -1, 0, 0, 1}, {2, 1, 0, 1, 1, 0},
731  {2, 3, 0, -1, 1, 0}, {2, 0, 1, 0, 0, -1},
732  {2, 0, -2, 2, 0, 0}, {2, -3, 0, 3, 0, 0},
733  {2, -2, 3, 0, 0, -1}, {2, 4, 0, 0, 0, 0},
734  {2, -1, 1, 1, 0, -1}, {2, -1, 3, -1, 0, -1},
735  {2, 2, 0, 0, -1, 0}, {2, -1, -1, 1, 0, 1},
736  {2, 4, 0, 0, 1, 0}, {2, -3, 4, -1, 0, 0},
737  {2, -1, 2, -1, -1, 0}, {2, 3, -2, 1, 0, 0},
738  {2, 1, 2, -1, 0, 0}, {2, -4, 2, 2, 0, 0},
739  {2, 4, -2, 0, 0, 0}, {2, 0, 2, 0, 0, 0},
740  {2, -2, 2, 0, -1, 0}, {2, 2, -4, 0, 0, 2},
741  {2, 2, -2, 0, -1, 0}, {2, 1, 0, -1, -1, 0},
742  {2, -1, 1, 0, 0, 0}, {2, 2, -1, 0, 0, 1},
743  {2, 2, 1, 0, 0, -1}, {2, -2, 0, 2, -1, 0},
744  {2, -2, 4, -2, 0, 0}, {2, 2, 2, 0, 0, 0},
745  {2, -4, 4, 0, 0, 0}, {2, -1, 0, -1, -2, 0},
746  {2, 1, 2, -1, 1, 0}, {2, -1, -2, 3, 0, 0},
747  {2, 3, -2, 1, 1, 0}, {2, 4, 0, -2, 0, 0},
748  {2, 0, 0, 2, 0, 0}, {2, 0, 2, -2, 0, 0},
749  {2, 0, 2, 0, 1, 0}, {2, -3, 3, 1, 0, -1},
750  {2, 0, 0, 0, -2, 0}, {2, 4, 0, 0, 2, 0},
751  {2, 4, -2, 0, 1, 0}, {2, 0, 0, 0, 0, 2},
752  {2, 1, 0, 1, 2, 0}, {2, 0, -2, 0, -2, 0},
753  {2, -2, 1, 0, 0, 1}, {2, -2, 1, 2, 0, -1},
754  {2, -1, 1, -1, 0, 1}, {2, 5, 0, -1, 0, 0},
755  {2, 1, -3, 1, 0, 1}, {2, -2, -1, 2, 0, 1},
756  {2, 3, 0, -1, 2, 0}, {2, 1, -2, 1, -1, 0},
757  {2, 5, 0, -1, 1, 0}, {2, -4, 0, 4, 0, 0},
758  {2, -3, 2, 1, -1, 0}, {2, -2, 1, 1, 0, 0},
759  {2, 4, 0, -2, 1, 0}, {2, 0, 0, 2, 1, 0},
760  {2, -5, 4, 1, 0, 0}, {2, 0, 2, 0, 2, 0},
761  {2, -1, 2, 1, 0, 0}, {2, 5, -2, -1, 0, 0},
762  {2, 1, -1, 0, 0, 0}, {2, 2, -2, 0, 0, 2},
763  {2, -5, 2, 3, 0, 0}, {2, -1, -2, 1, -2, 0},
764  {2, -3, 5, -1, 0, -1}, {2, -1, 0, 0, 0, 1},
765  {2, -2, 0, 0, -2, 0}, {2, 0, -1, 1, 0, 0},
766  {2, -3, 1, 1, 0, 1}, {2, 3, 0, -1, -1, 0},
767  {2, 1, 0, 1, -1, 0}, {2, -1, 2, 1, 1, 0},
768  {2, 0, -3, 2, 0, 1}, {2, 1, -1, -1, 0, 1},
769  {2, -3, 0, 3, -1, 0}, {2, 0, -2, 2, -1, 0},
770  {2, -4, 3, 2, 0, -1}, {2, -1, 0, 1, -2, 0},
771  {2, 5, 0, -1, 2, 0}, {2, -4, 5, 0, 0, -1},
772  {2, -2, 4, 0, 0, -2}, {2, -1, 0, 1, 0, 2},
773  {2, -2, -2, 4, 0, 0}, {2, 3, -2, -1, -1, 0},
774  {2, -2, 5, -2, 0, -1}, {2, 0, -1, 0, -1, 1},
775  {2, 5, -2, -1, 1, 0}, {1, 1, 0, 0, 0, 0},
776  {1, -1, 0, 0, 0, 0}, // x,K1,O1
777  {1, 1, -2, 0, 0, 0}, {1, -2, 0, 1, 0, 0},
778  {1, 1, 0, 0, 1, 0}, // P1,Q1,x
779  {1, -1, 0, 0, -1, 0}, {1, 2, 0, -1, 0, 0},
780  {1, 0, 0, 1, 0, 0}, {1, 3, 0, 0, 0, 0},
781  {1, -2, 2, -1, 0, 0}, {1, -2, 0, 1, -1, 0},
782  {1, -3, 2, 0, 0, 0}, {1, 0, 0, -1, 0, 0},
783  {1, 1, 0, 0, -1, 0}, {1, 3, 0, 0, 1, 0},
784  {1, 1, -3, 0, 0, 1}, {1, -3, 0, 2, 0, 0},
785  {1, 1, 2, 0, 0, 0}, {1, 0, 0, 1, 1, 0},
786  {1, 2, 0, -1, 1, 0}, {1, 0, 2, -1, 0, 0},
787  {1, 2, -2, 1, 0, 0}, {1, 3, -2, 0, 0, 0},
788  {1, -1, 2, 0, 0, 0}, {1, 1, 1, 0, 0, -1},
789  {1, 1, -1, 0, 0, 1}, {1, 4, 0, -1, 0, 0},
790  {1, -4, 2, 1, 0, 0}, {1, 0, -2, 1, 0, 0},
791  {1, -2, 2, -1, -1, 0}, {1, 3, 0, -2, 0, 0},
792  {1, -1, 0, 2, 0, 0}, {1, -1, 0, 0, -2, 0},
793  {1, 3, 0, 0, 2, 0}, {1, -3, 2, 0, -1, 0},
794  {1, 4, 0, -1, 1, 0}, {1, 0, 0, -1, -1, 0},
795  {1, 1, -2, 0, -1, 0}, {1, -3, 0, 2, -1, 0},
796  {1, 1, 0, 0, 2, 0}, {1, 1, -1, 0, 0, -1},
797  {1, -1, -1, 0, 0, 1}, {1, 0, 2, -1, 1, 0},
798  {1, -1, 1, 0, 0, -1}, {1, -1, -2, 2, 0, 0},
799  {1, 2, -2, 1, 1, 0}, {1, -4, 0, 3, 0, 0},
800  {1, -1, 2, 0, 1, 0}, {1, 3, -2, 0, 1, 0},
801  {1, 2, 0, -1, -1, 0}, {1, 0, 0, 1, -1, 0},
802  {1, -2, 2, 1, 0, 0}, {1, 4, -2, -1, 0, 0},
803  {1, -3, 3, 0, 0, -1}, {1, -2, 1, 1, 0, -1},
804  {1, -2, 3, -1, 0, -1}, {1, 0, -2, 1, -1, 0},
805  {1, -2, -1, 1, 0, 1}, {1, 4, -2, 1, 0, 0},
806  {1, -4, 4, -1, 0, 0}, {1, -4, 2, 1, -1, 0},
807  {1, 5, -2, 0, 0, 0}, {1, 3, 0, -2, 1, 0},
808  {1, -5, 2, 2, 0, 0}, {1, 2, 0, 1, 0, 0},
809  {1, 1, 3, 0, 0, -1}, {1, -2, 0, 1, -2, 0},
810  {1, 4, 0, -1, 2, 0}, {1, 1, -4, 0, 0, 2},
811  {1, 5, 0, -2, 0, 0}, {1, -1, 0, 2, 1, 0},
812  {1, -2, 1, 0, 0, 0}, {1, 4, -2, 1, 1, 0},
813  {1, -3, 4, -2, 0, 0}, {1, -1, 3, 0, 0, -1},
814  {1, 3, -3, 0, 0, 1}, {1, 5, -2, 0, 1, 0},
815  {1, 1, 2, 0, 1, 0}, {1, 2, 0, 1, 1, 0},
816  {1, -5, 4, 0, 0, 0}, {1, -2, 0, -1, -2, 0},
817  {1, 5, 0, -2, 1, 0}, {1, 1, 2, -2, 0, 0},
818  {1, 1, -2, 2, 0, 0}, {1, -2, 2, 1, 1, 0},
819  {1, 0, 3, -1, 0, -1}, {1, 2, -3, 1, 0, 1},
820  {1, -2, -2, 3, 0, 0}, {1, -1, 2, -2, 0, 0},
821  {1, -4, 3, 1, 0, -1}, {1, -4, 0, 3, -1, 0},
822  {1, -1, -2, 2, -1, 0}, {1, -2, 0, 3, 0, 0},
823  {1, 4, 0, -3, 0, 0}, {1, 0, 1, 1, 0, -1},
824  {1, 2, -1, -1, 0, 1}, {1, 2, -2, 1, -1, 0},
825  {1, 0, 0, -1, -2, 0}, {1, 2, 0, 1, 2, 0},
826  {1, 2, -2, -1, -1, 0}, {1, 0, 0, 1, 2, 0},
827  {1, 0, 1, 0, 0, 0}, {1, 2, -1, 0, 0, 0},
828  {1, 0, 2, -1, -1, 0}, {1, -1, -2, 0, -2, 0},
829  {1, -3, 1, 0, 0, 1}, {1, 3, -2, 0, -1, 0},
830  {1, -1, -1, 0, -1, 1}, {1, 4, -2, -1, 1, 0},
831  {1, 2, 1, -1, 0, -1}, {1, 0, -1, 1, 0, 1},
832  {1, -2, 4, -1, 0, 0}, {1, 4, -4, 1, 0, 0},
833  {1, -3, 1, 2, 0, -1}, {1, -3, 3, 0, -1, -1},
834  {1, 1, 2, 0, 2, 0}, {1, 1, -2, 0, -2, 0},
835  {1, 3, 0, 0, 3, 0}, {1, -1, 2, 0, -1, 0},
836  {1, -2, 1, -1, 0, 1}, {1, 0, -3, 1, 0, 1},
837  {1, -3, -1, 2, 0, 1}, {1, 2, 0, -1, 2, 0},
838  {1, 6, -2, -1, 0, 0}, {1, 2, 2, -1, 0, 0},
839  {1, -1, 1, 0, -1, -1}, {1, -2, 3, -1, -1, -1},
840  {1, -1, 0, 0, 0, 2}, {1, -5, 0, 4, 0, 0},
841  {1, 1, 0, 0, 0, -2}, {1, -2, 1, 1, -1, -1},
842  {1, 1, -1, 0, 1, 1}, {1, 1, 2, 0, 0, -2},
843  {1, -3, 1, 1, 0, 0}, {1, -4, 4, -1, -1, 0},
844  {1, 1, 0, -2, -1, 0}, {1, -2, -1, 1, -1, 1},
845  {1, -3, 2, 2, 0, 0}, {1, 5, -2, -2, 0, 0},
846  {1, 3, -4, 2, 0, 0}, {1, 1, -2, 0, 0, 2},
847  {1, -1, 4, -2, 0, 0}, {1, 2, 2, -1, 1, 0},
848  {1, -5, 2, 2, -1, 0}, {1, 1, -3, 0, -1, 1},
849  {1, 1, 1, 0, 1, -1}, {1, 6, -2, -1, 1, 0},
850  {1, -2, 2, -1, -2, 0}, {1, 4, -2, 1, 2, 0},
851  {1, -6, 4, 1, 0, 0}, {1, 5, -4, 0, 0, 0},
852  {1, -3, 4, 0, 0, 0}, {1, 1, 2, -2, 1, 0},
853  {1, -2, 1, 0, -1, 0}, {0, 2, 0, 0, 0, 0}, // x,x,Mf
854  {0, 1, 0, -1, 0, 0}, {0, 0, 2, 0, 0, 0},
855  {0, 0, 0, 0, 1, 0}, // Mm,SSa
856  {0, 2, 0, 0, 1, 0}, {0, 3, 0, -1, 0, 0},
857  {0, 1, -2, 1, 0, 0}, {0, 2, -2, 0, 0, 0},
858  {0, 3, 0, -1, 1, 0}, {0, 0, 1, 0, 0, -1},
859  {0, 2, 0, -2, 0, 0}, {0, 2, 0, 0, 2, 0},
860  {0, 3, -2, 1, 0, 0}, {0, 1, 0, -1, -1, 0},
861  {0, 1, 0, -1, 1, 0}, {0, 4, -2, 0, 0, 0},
862  {0, 1, 0, 1, 0, 0}, {0, 0, 3, 0, 0, -1},
863  {0, 4, 0, -2, 0, 0}, {0, 3, -2, 1, 1, 0},
864  {0, 3, -2, -1, 0, 0}, {0, 4, -2, 0, 1, 0},
865  {0, 0, 2, 0, 1, 0}, {0, 1, 0, 1, 1, 0},
866  {0, 4, 0, -2, 1, 0}, {0, 3, 0, -1, 2, 0},
867  {0, 5, -2, -1, 0, 0}, {0, 1, 2, -1, 0, 0},
868  {0, 1, -2, 1, -1, 0}, {0, 1, -2, 1, 1, 0},
869  {0, 2, -2, 0, -1, 0}, {0, 2, -3, 0, 0, 1},
870  {0, 2, -2, 0, 1, 0}, {0, 0, 2, -2, 0, 0},
871  {0, 1, -3, 1, 0, 1}, {0, 0, 0, 0, 2, 0},
872  {0, 0, 1, 0, 0, 1}, {0, 1, 2, -1, 1, 0},
873  {0, 3, 0, -3, 0, 0}, {0, 2, 1, 0, 0, -1},
874  {0, 1, -1, -1, 0, 1}, {0, 1, 0, 1, 2, 0},
875  {0, 5, -2, -1, 1, 0}, {0, 2, -1, 0, 0, 1},
876  {0, 2, 2, -2, 0, 0}, {0, 1, -1, 0, 0, 0},
877  {0, 5, 0, -3, 0, 0}, {0, 2, 0, -2, 1, 0},
878  {0, 1, 1, -1, 0, -1}, {0, 3, -4, 1, 0, 0},
879  {0, 0, 2, 0, 2, 0}, {0, 2, 0, -2, -1, 0},
880  {0, 4, -3, 0, 0, 1}, {0, 3, -1, -1, 0, 1},
881  {0, 0, 2, 0, 0, -2}, {0, 3, -3, 1, 0, 1},
882  {0, 2, -4, 2, 0, 0}, {0, 4, -2, -2, 0, 0},
883  {0, 3, 1, -1, 0, -1}, {0, 5, -4, 1, 0, 0},
884  {0, 3, -2, -1, -1, 0}, {0, 3, -2, 1, 2, 0},
885  {0, 4, -4, 0, 0, 0}, {0, 6, -2, -2, 0, 0},
886  {0, 5, 0, -3, 1, 0}, {0, 4, -2, 0, 2, 0},
887  {0, 2, 2, -2, 1, 0}, {0, 0, 4, 0, 0, -2},
888  {0, 3, -1, 0, 0, 0}, {0, 3, -3, -1, 0, 1},
889  {0, 4, 0, -2, 2, 0}, {0, 1, -2, -1, -1, 0},
890  {0, 2, -1, 0, 0, -1}, {0, 4, -4, 2, 0, 0},
891  {0, 2, 1, 0, 1, -1}, {0, 3, -2, -1, 1, 0},
892  {0, 4, -3, 0, 1, 1}, {0, 2, 0, 0, 3, 0},
893  {0, 6, -4, 0, 0, 0},
894  };
895 
896  if ((int)(sizeof(DerAmp) / sizeof(double)) != NDER ||
897  (int)(sizeof(DerInd) / sizeof(NVector)) != NDER)
898  {
899  Exception e("Static arrays are corrupted");
900  GNSSTK_THROW(e);
901  }
902 
903  int i, j, k, kk;
904  static const double dtr(0.01745329252);
905 
906  // get amplitude, phase and frequency for each of the standard tides
907  double RealAmp[NSTD], ImagAmp[NSTD], Freq[NSTD];
908  double phsrad, freq, phas;
909  for (i = 0; i < Nin; i++)
910  { // Nin is NSTD
911  // first find the index for this tide
912  j = stdindex[i];
913 
914  // amplitudes
915  phsrad = phs[i] * dtr; // DEG_TO_RAD; // phase in radians
916  RealAmp[i] = amp[i] * ::cos(phsrad) / ::fabs(DerAmp[j]);
917  ImagAmp[i] = amp[i] * ::sin(phsrad) / ::fabs(DerAmp[j]);
918  // LOG(INFO) << "TEST " << setw(2) << i+1 << fixed << setprecision(15)
919  // << setw(19) << amp[i] << setw(19) << phsrad << setw(19) <<
920  // DerAmp[j];
921 
922  // phase and freq
923  freq = phas = 0.0;
924  for (k = 0; k < 6; k++)
925  {
926  freq += DerInd[j].n[k] * freqDood[k];
927  // not used phas += DerInd[j].n[k] * Dood[k];
928  }
929  Freq[i] = freq;
930 
931  /* make 0 <= phas < 360 -- why?
932  not used phas = ::fmod(phas,360.0);
933  not used if(phas < 0.0) phas += 360.0;
934 
935  LOG(INFO) << "Dood " << setw(2) << i << " at index " << setw(3) << j
936  << " (" << DerInd[j].n[0] << "," << setw(2) << DerInd[j].n[1]
937  << "," << setw(2) << DerInd[j].n[2] << "," << setw(2) <<
938  DerInd[j].n[3]
939  << "," << setw(2) << DerInd[j].n[4] << "," << setw(2) <<
940  DerInd[j].n[5]
941  << ") rA iA F " << fixed << setprecision(10) << setw(13)
942  << RealAmp[i] << " " << setw(13) << ImagAmp[i]
943  << " " << setw(12) << Freq[i]; */
944  }
945 
946  // sort the frequency, and keep the key
947  int key[NSTD];
948  for (i = 0; i < Nin; ++i)
949  key[i] = i;
950  QSort(Freq, key, Nin);
951 
952  // use key to sort amplitudes
953  double tmpR[NSTD], tmpI[NSTD];
954  for (i = 0; i < Nin; ++i)
955  {
956  tmpR[i] = RealAmp[i];
957  tmpI[i] = ImagAmp[i];
958  }
959  for (i = 0; i < Nin; ++i)
960  {
961  RealAmp[i] = tmpR[key[i]];
962  ImagAmp[i] = tmpI[key[i]];
963  }
964 
965  // count the shells
966  int nl(0), nm(0), nh(0);
967  // LOG(INFO) << "TEST Sorted reamp, imamp, freq\n";
968  for (i = 0; i < Nin; i++)
969  { // Nin is NSTD
970  /* LOG(INFO) << "Sorted Dood " << setw(2) << key[i] << " rA iA F P "
971  << fixed << setprecision(10) << setw(13) << RealAmp[i]
972  << " " << setw(13) << ImagAmp[i] << " " << setw(12) << Freq[i];
973  LOG(INFO) << "TEST " << setw(2) << i+1 << fixed << setprecision(15)
974  << setw(19)<< RealAmp[i] << setw(19) << ImagAmp[i] << setw(19) <<
975  Freq[i]; */
976  if (Freq[i] < 0.5)
977  {
978  nl++;
979  }
980  else if (Freq[i] < 1.5)
981  {
982  nm++;
983  }
984  else if (Freq[i] < 2.5)
985  {
986  nh++;
987  }
988  // so freq cannot be >= 2.5??
989  }
990  // LOG(INFO) << "Shells contain " << nl << " " << nm << " " << nh;
991 
992  // split arrays into vector<double> for each shell
993  vector<double> Flow, Rlow, Ilow, Fmed, Rmed, Imed, Fhi, Rhi, Ihi;
994  for (i = 0; i < nl; i++)
995  {
996  Flow.push_back(Freq[i]);
997  Rlow.push_back(RealAmp[i]);
998  Ilow.push_back(ImagAmp[i]);
999  // LOG(INFO) << "Low shell Dood " << setw(2) << key[i] << " rA iA F "
1000  // << fixed << setprecision(10) << setw(13) << RealAmp[i]
1001  // << " " << setw(13) << ImagAmp[i] << " " << setw(12) << Freq[i];
1002  }
1003  for (i = nl; i < nl + nm; i++)
1004  {
1005  Fmed.push_back(Freq[i]);
1006  Rmed.push_back(RealAmp[i]);
1007  Imed.push_back(ImagAmp[i]);
1008  // LOG(INFO) << "Med shell Dood " << setw(2) << key[i] << " rA iA F "
1009  // << fixed << setprecision(10) << setw(13) << RealAmp[i]
1010  // << " " << setw(13) << ImagAmp[i] << " " << setw(12) << Freq[i];
1011  }
1012  for (i = nl + nm; i < nl + nm + nh; i++)
1013  {
1014  Fhi.push_back(Freq[i]);
1015  Rhi.push_back(RealAmp[i]);
1016  Ihi.push_back(ImagAmp[i]);
1017  // LOG(INFO) << "Hi shell Dood " << setw(2) << key[i] << " rA iA F "
1018  // << fixed << setprecision(10) << setw(13) << RealAmp[i]
1019  // << " " << setw(13) << ImagAmp[i] << " " << setw(12) << Freq[i];
1020  }
1021 
1022  // find splines of amp vs frequency in each shell
1023  CubicSpline<double> csRlow, csIlow, csRmed, csImed, csRhi, csIhi;
1024  if (nl > 0)
1025  {
1026  csRlow.initialize(Flow, Rlow);
1027  csIlow.initialize(Flow, Ilow);
1028  }
1029  csRmed.initialize(Fmed, Rmed);
1030  csImed.initialize(Fmed, Imed);
1031  csRhi.initialize(Fhi, Rhi);
1032  csIhi.initialize(Fhi, Ihi);
1033 
1034  // evaluate splines at each of the NDER waves; not all will contribute
1035  int nout(0);
1036  for (j = 0; j < NDER; j++)
1037  { // loop over 342 derived tides
1038  // this is why nout may be < NDER
1039  if (DerInd[j].n[0] == 0 && nl == 0)
1040  {
1041  continue;
1042  }
1043 
1044  // get phase and freq for this tide
1045  freqDer[nout] = phsDer[nout] = 0.0;
1046  for (k = 0; k < 6; k++)
1047  {
1048  freqDer[nout] += DerInd[j].n[k] * freqDood[k];
1049  phsDer[nout] += DerInd[j].n[k] * Dood[k];
1050  }
1051  phsDer[nout] = ::fmod(phsDer[nout], 360.0);
1052  if (phsDer[nout] < 0.0)
1053  {
1054  phsDer[nout] += 360.0;
1055  }
1056 
1057  // LOG(INFO) << "TEST TDFRPH "
1058  // << setw(3) << j+1 << fixed << setprecision(15)
1059  // << setw(22) << freqDer[nout] << setw(22) << phsDer[nout]
1060  // << setw(3) << DerInd[j].n[0] << setw(3) << DerInd[j].n[1]
1061  // << setw(3) << DerInd[j].n[2] << setw(3) << DerInd[j].n[3]
1062  // << setw(3) << DerInd[j].n[4] << setw(3) << DerInd[j].n[5];
1063 
1064  if (DerInd[j].n[0] == 0)
1065  {
1066  phsDer[nout] += 180.0;
1067  }
1068  else if (DerInd[j].n[0] == 1)
1069  {
1070  phsDer[nout] += 90.0;
1071  }
1072 
1073  // get amplitudes at freq
1074  freq = freqDer[nout];
1075  double ramp, iamp;
1076  if (DerInd[j].n[0] == 0)
1077  {
1078  if (csRlow.testLimits(freq, ramp))
1079  {
1080  ramp = csRlow.evaluate(freq);
1081  }
1082  if (csIlow.testLimits(freq, iamp))
1083  {
1084  iamp = csIlow.evaluate(freq);
1085  }
1086  }
1087  else if (DerInd[j].n[0] == 1)
1088  {
1089  if (csRmed.testLimits(freq, ramp))
1090  {
1091  ramp = csRmed.evaluate(freq);
1092  }
1093  if (csImed.testLimits(freq, iamp))
1094  {
1095  iamp = csImed.evaluate(freq);
1096  }
1097  }
1098  else if (DerInd[j].n[0] == 2)
1099  {
1100  if (csRhi.testLimits(freq, ramp))
1101  {
1102  ramp = csRhi.evaluate(freq);
1103  }
1104  if (csIhi.testLimits(freq, iamp))
1105  {
1106  iamp = csIhi.evaluate(freq);
1107  }
1108  }
1109 
1110  ampDer[nout] = DerAmp[j] * RSS(ramp, iamp);
1111  phsDer[nout] += ::atan2(iamp, ramp) / dtr; //*RAD_TO_DEG; // TEMP
1112  if (phsDer[nout] > 180.0)
1113  phsDer[nout] -= 360.0;
1114 
1115  // LOG(INFO) << "TEST RE AM " << setw(3) << j+1 << fixed
1116  // << setprecision(15) << setw(22) << ramp << setw(22) << iamp
1117  // << setw(22) << ampDer[nout] << setw(22) << phsDer[nout];
1118 
1119  nout++;
1120  }
1121 
1122  return nout;
1123 
1124  } // end int OceanLoadTides::deriveTides()
1125 
1126 } // end namespace gnsstk
1127 //------------------------------------------------------------------------------------
1128 //------------------------------------------------------------------------------------
gnsstk::EphTime::convertSystemTo
void convertSystemTo(const TimeSystem &ts)
Definition: EphTime.hpp:96
amp
static const double amp[]
Amplitudes (microarcsec); indexed using the iamp array below.
Definition: IERS2010CIOSeriesData.hpp:1630
gnsstk::OceanLoadTides::NVector::n
int n[6]
Definition: OceanLoadTides.hpp:187
CubicSpline.hpp
coeff
static const struct @1 coeff[]
Constants used in the nutation models, adapted from SOFA nut80.c.
gnsstk::StringUtils::word
std::string word(const std::string &s, const std::string::size_type wordNum=0, const char delimiter=' ')
Definition: StringUtils.hpp:1112
gnsstk::CubicSpline
Cubic spline interpolation.
Definition: CubicSpline.hpp:56
StringUtils.hpp
gnsstk::SEC_PER_DAY
const long SEC_PER_DAY
Seconds per day.
Definition: TimeConstants.hpp:63
gnsstk::MJD_JDAY
const long MJD_JDAY
'Julian day' offset from MJD
Definition: TimeConstants.hpp:53
gnsstk::CubicSpline::evaluate
T evaluate(const T &x)
Definition: CubicSpline.hpp:140
gnsstk::CubicSpline::testLimits
bool testLimits(const T &x, T &y)
Definition: CubicSpline.hpp:109
GNSSconstants.hpp
gnsstk::Triple
Definition: Triple.hpp:68
gnsstk::OceanLoadTides::NVector
Used for convenience by computeDisplacements.
Definition: OceanLoadTides.hpp:185
gnsstk
For Sinex::InputHistory.
Definition: BasicFramework.cpp:50
gnsstk::StringUtils::stripLeading
std::string & stripLeading(std::string &s, const std::string &aString, std::string::size_type num=std::string::npos)
Definition: StringUtils.hpp:1426
std::sin
double sin(gnsstk::Angle x)
Definition: Angle.hpp:144
gnsstk::CubicSpline::initialize
void initialize(const std::vector< T > &X, const std::vector< T > &Y)
Definition: CubicSpline.hpp:78
MiscMath.hpp
example4.temp
temp
Definition: example4.py:35
gnsstk::Exception
Definition: Exception.hpp:151
gnsstk::StringUtils::stripTrailing
std::string & stripTrailing(std::string &s, const std::string &aString, std::string::size_type num=std::string::npos)
Definition: StringUtils.hpp:1453
gnsstk::StringUtils::numWords
int numWords(const std::string &s, const char delimiter=' ')
Definition: StringUtils.hpp:2171
gnsstk::convertJDtoCalendar
void convertJDtoCalendar(long jd, int &iyear, int &imonth, int &iday)
Definition: TimeConverters.cpp:52
example4.time
time
Definition: example4.py:103
gnsstk::EphTime
Definition: EphTime.hpp:67
nl
int nl
Definition: IERS1996NutationData.hpp:44
gnsstk::StringUtils::stripFirstWord
std::string stripFirstWord(std::string &s, const char delimiter=' ')
Definition: StringUtils.hpp:2253
std::cos
double cos(gnsstk::Angle x)
Definition: Angle.hpp:146
gnsstk::StringUtils::asDouble
double asDouble(const std::string &s)
Definition: StringUtils.hpp:705
gnsstk::RSS
T RSS(T aa, T bb, T cc)
Perform the root sum square of aa, bb and cc.
Definition: MiscMath.hpp:246
GNSSTK_RETHROW
#define GNSSTK_RETHROW(exc)
Definition: Exception.hpp:369
example4.pos
pos
Definition: example4.py:125
OceanLoadTides.hpp
RobustStats.hpp
iamp
static const int iamp[]
indexes into the amplitudes array for each frequency; length = NFALS + NFAP
Definition: IERS2010CIOSeriesData.hpp:1509
std
Definition: Angle.hpp:142
gnsstk::QSort
void QSort(T *sa, int na, int(*comp)(const T &, const T &)=gnsstk::Qsort_compare)
Definition: RobustStats.hpp:139
gnsstk::EphTime::secOfDay
double secOfDay() const
Definition: EphTime.hpp:177
gnsstk::StringUtils::firstWord
std::string firstWord(const std::string &s, const char delimiter=' ')
Definition: StringUtils.hpp:2138
GNSSTK_THROW
#define GNSSTK_THROW(exc)
Definition: Exception.hpp:366
gnsstk::EphTime::dMJD
double dMJD() const
Definition: EphTime.hpp:171
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:40