SolarSystemEphemeris.cpp
Go to the documentation of this file.
1 //==============================================================================
2 //
3 // This file is part of GNSSTk, the ARL:UT GNSS Toolkit.
4 //
5 // The GNSSTk is free software; you can redistribute it and/or modify
6 // it under the terms of the GNU Lesser General Public License as published
7 // by the Free Software Foundation; either version 3.0 of the License, or
8 // any later version.
9 //
10 // The GNSSTk is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU Lesser General Public License for more details.
14 //
15 // You should have received a copy of the GNU Lesser General Public
16 // License along with GNSSTk; if not, write to the Free Software Foundation,
17 // Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA
18 //
19 // This software was developed by Applied Research Laboratories at the
20 // University of Texas at Austin.
21 // Copyright 2004-2022, The Board of Regents of The University of Texas System
22 //
23 //==============================================================================
24 
25 //==============================================================================
26 //
27 // This software was developed by Applied Research Laboratories at the
28 // University of Texas at Austin, under contract to an agency or agencies
29 // within the U.S. Department of Defense. The U.S. Government retains all
30 // rights to use, duplicate, distribute, disclose, or release this software.
31 //
32 // Pursuant to DoD Directive 523024
33 //
34 // DISTRIBUTION STATEMENT A: This software has been approved for public
35 // release, distribution is unlimited.
36 //
37 //==============================================================================
38 
50 //------------------------------------------------------------------------------------
51 #include "SolarSystemEphemeris.hpp"
52 // GNSSTk
53 #include "FormattedDouble.hpp"
54 #include "StringUtils.hpp"
55 #include "TimeConverters.hpp"
56 #include "logstream.hpp"
57 
58 //------------------------------------------------------------------------------------
59 using namespace std;
60 using namespace gnsstk::StringUtils;
61 
62 namespace gnsstk
63 {
65  class SSEDouble : public FormattedDouble
66  {
67  public:
68  SSEDouble(double d = 0)
69  : FormattedDouble(d, StringUtils::FFLead::Decimal, 18, 2, 26, 'D',
70  StringUtils::FFSign::NegOnly,
71  StringUtils::FFAlign::Right)
72  {
73  }
75  SSEDouble(const string& str) : FormattedDouble(str, 26, 'D') {}
76 
78  SSEDouble& operator=(const string& s)
79  {
80  FormattedDouble::operator=(s);
81  return *this;
82  }
83  };
84 
85  //---------------------------------------------------------------------------------
86  void SolarSystemEphemeris::readASCIIheader(const string& filename)
87  {
88  try
89  {
90  // cout << "Reporting in SolarSystemEphemeris is "
91  // << ConfigureLOG::ToString(ConfigureLOG::ReportingLevel()) << endl;
92 
93  // open the file
94  ifstream strm;
95  strm.open(filename.c_str());
96  if (!strm.is_open())
97  {
98  Exception e("Failed to open input file " + filename + ". Abort.");
99  GNSSTK_THROW(e);
100  }
101 
102  // clear existing data
103  constants.clear();
104 
105  /* read the file one line at a time, process depending on the value of
106  group */
107  int group = 0, n = 0; // n will count lines/items within a group
108  string line, word;
109  vector<string>
110  const_names; // Temporary - name,value stored in map constants
111  while (1)
112  {
113  getline(strm, line);
114  stripTrailing(line, '\r');
115 
116  // catch new groups
117  if (line.substr(0, 5) == "GROUP")
118  {
119  word = stripFirstWord(line);
120  group = asInt(stripFirstWord(line));
121  LOG(DEBUG) << "Group is " << group; // temp
122  n = 0; // n will count lines/items within a group
123  continue;
124  }
125 
126  // skip blank lines
127  stripLeading(line, " ");
128  if (line.empty())
129  {
130  if (strm.eof() || !strm.good())
131  {
132  break; // if the last line is blank
133  }
134  else
135  {
136  continue;
137  }
138  }
139 
140  /* process entire line at once
141  first line (no GROUP) */
142  if (group == 0)
143  {
144  word = stripFirstWord(line);
145  word = stripFirstWord(line); // ignore KSIZE
146  word = stripFirstWord(line);
147  if (word == "NCOEFF=")
148  {
149  Ncoeff = asInt(stripFirstWord(line));
150  LOG(DEBUG) << "Ncoeff is " << Ncoeff; // temp
151  continue;
152  }
153  else
154  {
155  Exception e(
156  "Confused on the first line - 3rd word is not NCOEFF=");
157  GNSSTK_THROW(e);
158  }
159  }
160  // GROUP 1010
161  else if (group == 1010)
162  {
163  if (n > 2)
164  { // this should not happen
165  Exception e("Too many labels under GROUP 1010");
166  GNSSTK_THROW(e);
167  }
168  else
169  {
170  stripTrailing(line, " ");
171  label[n++] = line;
172  LOG(VERBOSE) << "label " << n << " is " << line; // temp
173  continue;
174  }
175  }
176  // GROUP 1030
177  else if (group == 1030)
178  {
179  /* start and stop times. These are meaningless here, because they
180  will be determined by the data that follows this header, and
181  so are meaningful only in the binary file. */
182  startJD = SSEDouble(stripFirstWord(line));
183  endJD = SSEDouble(stripFirstWord(line));
184  // interval in days covered by each block of coefficients
185  interval = SSEDouble(stripFirstWord(line));
186  LOG(VERBOSE)
187  << "Times JD " << fixed << setprecision(3) << startJD
188  << " to JD " << endJD << " = " << interval << " days";
189  }
190  // GROUP 1070 - end-of-header
191  else if (group == 1070)
192  {
193  break;
194  }
195 
196  // process the line one (whitespace-separated) word at a time
197  while (!line.empty())
198  {
199  word = stripFirstWord(line);
200 
201  if (group == 1040)
202  {
203  if (n++ == 0)
204  {
205  Nconst = asInt(word);
206  LOG(DEBUG) << "Nconst is " << Nconst;
207  }
208  else
209  {
210  const_names.push_back(word);
211  }
212  }
213  else if (group == 1041)
214  {
215  if (n++ == 0)
216  {
217  if (Nconst != asInt(word))
218  {
219  Exception e("Nconst does not match N in GROUP 1041 : " +
220  asString(Nconst) + " != " + word);
221  GNSSTK_THROW(e);
222  }
223  LOG(DEBUG)
224  << "Nconst matches: " << Nconst << " = " << word;
225  }
226  else
227  {
228  constants[const_names[n - 2]] = SSEDouble(word);
229  }
230  }
231  else if (group == 1050)
232  {
233  if (n < 13)
234  {
235  c_offset[n] = asInt(word);
236  LOG(DEBUG) << "c_offset[" << n << "] = " << c_offset[n];
237  }
238  else if (n < 26)
239  {
240  c_ncoeff[n - 13] = asInt(word);
241  LOG(DEBUG)
242  << "c_ncoeff[" << n - 13 << "] = " << c_ncoeff[n - 13];
243  }
244  else
245  {
246  c_nsets[n - 26] = asInt(word);
247  LOG(DEBUG)
248  << "c_nsets[" << n - 26 << "] = " << c_nsets[n - 26];
249  }
250  n++;
251  }
252  else
253  {
254  Exception e("Confused about GROUP : " + asString(group));
255  GNSSTK_THROW(e);
256  }
257  } // end loop over words
258 
259  if (strm.eof() || !strm.good())
260  {
261  break; // if the last line is not blank
262  }
263 
264  } // end read loop over lines
265 
266  strm.clear();
267  strm.close();
268 
269  // test that we got a full header
270  if (group != 1070)
271  {
272  Exception e("Premature end of header");
273  GNSSTK_THROW(e);
274  }
275 
276  // EphemerisNumber != -1 means the header is complete
277  EphemerisNumber = int(constants["DENUM"]);
278 
279  // clear the data arrays
280  store.clear();
281  }
282  catch (Exception& e)
283  {
284  GNSSTK_RETHROW(e);
285  }
286  catch (exception& e)
287  {
288  Exception E("std except: " + string(e.what()));
289  GNSSTK_THROW(E);
290  }
291  catch (...)
292  {
293  Exception e("Unknown exception");
294  GNSSTK_THROW(e);
295  }
296  } // end void SolarSystemEphemeris::readASCIIheader(const string& filename)
297 
298  //---------------------------------------------------------------------------------
299  int SolarSystemEphemeris::readASCIIdata(vector<string>& filenames)
300  {
301  try
302  {
303  int i, n;
304  double jd;
305 
306  if (filenames.size() == 0)
307  {
308  return 0;
309  }
310 
311  // read the files, in any order; the data map will be sorted on time
312  for (i = 0; i < filenames.size(); i++)
313  {
314  int iret = readASCIIdata(filenames[i]);
315  if (iret)
316  {
317  return iret;
318  }
319  }
320 
321  // set the start and stop times in the header
322  map<double, vector<double>>::iterator it = store.begin();
323  startJD = it->second[0];
324  it = store.end();
325  it--;
326  endJD = it->second[1];
327 
328  LOG(VERBOSE) << "After reading data files, store size is "
329  << store.size() << " and new start and stop times are JD "
330  << fixed << setprecision(9) << startJD << " and JD "
331  << endJD;
332 
333  // Mod the header labels to reflect the new time limits
334  ostringstream oss;
335  int yy, mm, dd;
336  convertJDtoCalendar(static_cast<long>(startJD + 0.5), yy, mm, dd);
337  oss << "Start Epoch: JED= " << fixed << setw(10) << setprecision(1)
338  << startJD << " = " << yy << "/" << mm << "/" << dd;
339  label[1] = leftJustify(oss.str(), 81);
340  oss.seekp(ios_base::beg);
341 
342  convertJDtoCalendar(static_cast<long>(endJD + 0.5), yy, mm, dd);
343  oss.str("");
344  oss << "Final Epoch: JED= " << fixed << setw(10) << setprecision(1)
345  << endJD << " = " << yy << "/" << mm << "/" << dd;
346  label[2] = leftJustify(oss.str(), 81);
347 
348  LOG(VERBOSE) << "New label 1 is " << stripTrailing(label[1]);
349  LOG(VERBOSE) << "New label 2 is " << stripTrailing(label[2]);
350 
351  return 0;
352  }
353  catch (Exception& e)
354  {
355  GNSSTK_RETHROW(e);
356  }
357  catch (exception& e)
358  {
359  Exception E("std except: " + string(e.what()));
360  GNSSTK_THROW(E);
361  }
362  catch (...)
363  {
364  Exception e("Unknown exception");
365  GNSSTK_THROW(e);
366  }
367  } // end int SolarSystemEphemeris::readASCIIdata(vector<string>& filenames)
368 
369  //---------------------------------------------------------------------------------
370  int SolarSystemEphemeris::readASCIIdata(const string& filename)
371  {
372  try
373  {
374  if (EphemerisNumber < 0)
375  {
376  Exception e("readASCIIdata called before header read");
377  GNSSTK_THROW(e);
378  }
379 
380  int iret = 0;
381  string line, word;
382  ifstream strm;
383 
384  // open the file
385  strm.open(filename.c_str());
386  if (!strm.is_open())
387  {
388  Exception e("Could not open file " + filename);
389  GNSSTK_THROW(e);
390  }
391 
392  // expect this many lines per record
393  int nmax = Ncoeff / 3 + (Ncoeff % 3 ? 1 : 0);
394 
395  // loop over lines in the file
396  int ntot = 0; // counts the total number of lines
397  int n = 0; // counts the lines within a set of coefficients
398  int nc = 0; // count coefficients within a record
399  int rec;
400  vector<double> dataVector;
401  while (1)
402  {
403  getline(strm, line);
404  stripTrailing(line, '\r');
405 
406  if (line.empty())
407  {
408  if (strm.eof())
409  {
410  break;
411  }
412  if (!strm.good())
413  {
414  iret = retEarly;
415  break;
416  }
417  continue;
418  }
419 
420  if (n == 0)
421  {
422  rec =
423  asInt(stripFirstWord(line)); // 1st word is the record number
424  if (rec % 25 == 0)
425  {
426  LOG(VERBOSE) << "Record number " << rec;
427  }
428  int ncc = asInt(stripFirstWord(line)); // 2nd word is ncoeff
429  if (ncc != Ncoeff)
430  {
431  Exception e(
432  "readASCIIdata finds conflicting sizes in header (" +
433  asString(Ncoeff) + ") and data (" + asString(ncc) +
434  ") in file " + filename + " at line #" + asString(ntot));
435  GNSSTK_THROW(e);
436  }
437  nc = 0;
438  }
439  else
440  {
441  for (int j = 0; j < 3; j++)
442  {
443  double coeff = SSEDouble(stripFirstWord(line));
444  nc++;
445  dataVector.push_back(coeff);
446  if (nc >= Ncoeff)
447  {
448  vector<double> dtemp = dataVector;
449  store[dataVector[0]] = dtemp;
450  dataVector.clear();
451  break;
452  }
453  }
454  }
455 
456  if (strm.eof())
457  {
458  break;
459  }
460  if (!strm.good())
461  {
462  iret = retStrm;
463  break;
464  }
465 
466  if (n == nmax)
467  {
468  n = 0;
469  }
470  else
471  {
472  n++;
473  }
474  ntot++;
475  }
476  LOG(INFO) << "Read " << rec << " records from file " << filename;
477 
478  strm.close();
479 
480  return iret;
481  }
482  catch (Exception& e)
483  {
484  GNSSTK_RETHROW(e);
485  }
486  catch (exception& e)
487  {
488  Exception E("std except: " + string(e.what()));
489  GNSSTK_THROW(E);
490  }
491  catch (...)
492  {
493  Exception e("Unknown exception");
494  GNSSTK_THROW(e);
495  }
496  } // end int SolarSystemEphemeris::readASCIIdata(string filename)
497 
498  //---------------------------------------------------------------------------------
499  int SolarSystemEphemeris::writeASCIIheader(ostream& os)
500  {
501  try
502  {
503  if (EphemerisNumber < 0)
504  {
505  return retEphN;
506  }
507 
508  int i;
509  string str;
510  string blank(81, ' ');
511  blank += string("\n");
512  ostringstream oss;
513 
514  oss << "KSIZE= 0000 NSIZE=" << setw(5) << Ncoeff << blank;
515  os << leftJustify(oss.str(), 81) << endl << blank;
516  oss.seekp(ios_base::beg);
517 
518  os << leftJustify("GROUP 1010", 81) << endl << blank;
519  for (i = 0; i < 3; i++)
520  {
521  str = label[i];
522  os << leftJustify(str, 81) << endl;
523  }
524  os << blank;
525 
526  os << leftJustify("GROUP 1030", 81) << endl << blank;
527  oss << fixed << setprecision(2) << setw(12) << startJD << setw(12)
528  << endJD << setw(12) << interval << blank;
529  os << leftJustify(oss.str(), 81) << endl << blank;
530  oss.seekp(ios_base::beg);
531 
532  os << leftJustify("GROUP 1040", 81) << endl << blank;
533  oss << setw(6) << Nconst << blank;
534  os << leftJustify(oss.str(), 81) << endl;
535  oss.seekp(ios_base::beg);
536 
537  map<string, double>::const_iterator it = constants.begin();
538  for (i = 0; it != constants.end(); it++, i++)
539  {
540  oss << leftJustify(" " + it->first, 8);
541  if ((i + 1) % 10 == 0)
542  {
543  oss << blank;
544  os << leftJustify(oss.str(), 81) << endl;
545  oss.seekp(ios_base::beg);
546  }
547  }
548  if (Nconst % 10 != 0)
549  {
550  oss << blank;
551  os << leftJustify(oss.str(), 81) << endl;
552  oss.seekp(ios_base::beg);
553  }
554  os << blank;
555 
556  os << leftJustify("GROUP 1041", 81) << endl << blank;
557  oss << setw(6) << Nconst << blank;
558  os << leftJustify(oss.str(), 81) << endl;
559  oss.seekp(ios_base::beg);
560 
561  for (i = 0, it = constants.begin(); it != constants.end(); it++, i++)
562  {
563  oss << SSEDouble(it->second);
564  if ((i + 1) % 3 == 0)
565  {
566  oss << blank;
567  os << leftJustify(oss.str(), 81) << endl;
568  oss.seekp(ios_base::beg);
569  }
570  }
571  if (Nconst % 3 != 0)
572  {
573  i--;
574  while ((i + 1) % 3 != 0)
575  {
576  oss << SSEDouble(0.0);
577  i++;
578  }
579  oss << blank;
580  os << leftJustify(oss.str(), 81) << endl;
581  oss.seekp(ios_base::beg);
582  }
583  os << blank;
584 
585  os << leftJustify("GROUP 1050", 81) << endl << blank;
586  for (i = 0; i < 13; i++)
587  oss << rightJustify(asString(c_offset[i]), 6);
588  oss << blank;
589  os << leftJustify(oss.str(), 81) << endl;
590  oss.seekp(ios_base::beg);
591  for (i = 0; i < 13; i++)
592  oss << rightJustify(asString(c_ncoeff[i]), 6);
593  oss << blank;
594  os << leftJustify(oss.str(), 81) << endl;
595  oss.seekp(ios_base::beg);
596  for (i = 0; i < 13; i++)
597  oss << rightJustify(asString(c_nsets[i]), 6);
598  oss << blank;
599  os << leftJustify(oss.str(), 81) << endl;
600  oss.seekp(ios_base::beg);
601  os << blank;
602 
603  os << leftJustify("GROUP 1070", 81) << endl << blank;
604  os << blank;
605 
606  return 0;
607  }
608  catch (Exception& e)
609  {
610  GNSSTK_RETHROW(e);
611  }
612  catch (exception& e)
613  {
614  Exception E("std except: " + string(e.what()));
615  GNSSTK_THROW(E);
616  }
617  catch (...)
618  {
619  Exception e("Unknown exception");
620  GNSSTK_THROW(e);
621  }
622  }
623 
624  //---------------------------------------------------------------------------------
625  int SolarSystemEphemeris::writeASCIIdata(ostream& os)
626  {
627  try
628  {
629  if (EphemerisNumber < 0)
630  {
631  return retEphN;
632  }
633 
634  string blank(81, ' ');
635  blank += string("\n");
636  int i, nrec;
637  ostringstream oss;
638 
639  map<double, vector<double>>::iterator jt;
640  for (nrec = 1, jt = store.begin(); jt != store.end(); jt++, nrec++)
641  {
642  os << setw(6) << nrec << setw(6) << Ncoeff << " " << endl;
643  for (i = 0; i < Ncoeff; i++)
644  {
645  oss << SSEDouble(jt->second[i]);
646  if ((i + 1) % 3 == 0)
647  {
648  oss << blank;
649  os << leftJustify(oss.str(), 81) << endl;
650  oss.seekp(ios_base::beg);
651  }
652  }
653  if (Ncoeff % 3 != 0)
654  {
655  i--;
656  while ((i + 1) % 3 != 0)
657  {
658  oss << SSEDouble(0.0);
659  i++;
660  }
661  oss << blank;
662  os << leftJustify(oss.str(), 81) << endl;
663  oss.seekp(ios_base::beg);
664  }
665  }
666 
667  // TD clear the array after writing
668  // store.clear();
669 
670  return 0;
671  }
672  catch (Exception& e)
673  {
674  GNSSTK_RETHROW(e);
675  }
676  catch (exception& e)
677  {
678  Exception E("std except: " + string(e.what()));
679  GNSSTK_THROW(E);
680  }
681  catch (...)
682  {
683  Exception e("Unknown exception");
684  GNSSTK_THROW(e);
685  }
686  }
687 
688  //---------------------------------------------------------------------------------
689  int SolarSystemEphemeris::writeBinaryFile(const string& filename)
690  {
691  try
692  {
693  int i, recLength;
694  string str;
695 
696  if (EphemerisNumber <= 0)
697  {
698  return retEphN;
699  }
700 
701  // open the output file
702  ofstream strm;
703  strm.open(filename.c_str(), ios::out | ios::binary);
704  if (!strm.is_open())
705  {
706  Exception e("Failed to open output file " + filename + ". Abort.");
707  GNSSTK_THROW(e);
708  }
709 
710  /* write the header records: two of them, both of length
711  Ncoeff*sizeof(double) the structure and ordering taken from the JPL
712  code */
713  recLength = 0;
714 
715  /* first header record
716  -------------------------------------------------
717  1. 3 labels each of length 84 */
718  for (i = 0; i < 3; i++)
719  {
720  str = label[i];
721  writeBinary(strm, leftJustify(str, 84).c_str(), 84);
722  recLength += 84;
723  LOG(VERBOSE) << "writeBinaryFile label " << i << " " << label[i];
724  }
725 
726  // 2. 400 keys from the const array, each of length 6
727  map<string, double>::const_iterator it = constants.begin();
728  for (i = 0; i < 400; i++)
729  {
730  if (it != constants.end())
731  {
732  str = it->first;
733  writeBinary(strm, leftJustify(str, 6).c_str(), 6);
734  it++;
735  }
736  else
737  {
738  writeBinary(strm, " ", 6);
739  }
740  recLength += 6;
741  }
742 
743  // 3. the three times
744  writeBinary(strm, (char *)&startJD, sizeof(double));
745  writeBinary(strm, (char *)&endJD, sizeof(double));
746  writeBinary(strm, (char *)&interval, sizeof(double));
747  recLength += 3 * sizeof(double);
748 
749  // 4. Ncoeff
750  writeBinary(strm, (char *)&Ncoeff, sizeof(int));
751  recLength += sizeof(int);
752 
753  // 5. AU and EMRAT
754  writeBinary(strm, (char *)&constants["AU"], sizeof(double));
755  writeBinary(strm, (char *)&constants["EMRAT"], sizeof(double));
756  recLength += 2 * sizeof(double);
757 
758  // 6. c_arrays for the first 12 planets
759  for (i = 0; i < 12; i++)
760  {
761  writeBinary(strm, (char *)&c_offset[i], sizeof(int));
762  writeBinary(strm, (char *)&c_ncoeff[i], sizeof(int));
763  writeBinary(strm, (char *)&c_nsets[i], sizeof(int));
764  recLength += 3 * sizeof(int);
765  }
766 
767  // 7. DENUM
768  writeBinary(strm, (char *)&constants["DENUM"], sizeof(double));
769  recLength += sizeof(double);
770  LOG(VERBOSE) << "WriteBinary outputs DENUM " << constants["DENUM"];
771 
772  // 8. c_arrays for libration
773  writeBinary(strm, (char *)&c_offset[12], sizeof(int));
774  writeBinary(strm, (char *)&c_ncoeff[12], sizeof(int));
775  writeBinary(strm, (char *)&c_nsets[12], sizeof(int));
776  recLength += 3 * sizeof(int);
777 
778  // 9. pad
779  LOG(DEBUG) << "Pad length 1 = " << Ncoeff * sizeof(double) - recLength;
780  char c = ' ';
781  for (i = 0; i < Ncoeff * sizeof(double) - recLength; i++)
782  writeBinary(strm, &c, 1);
783 
784  /* second header record
785  ------------------------------------------------- the second header
786  record: 400 values from the const array */
787  double z = 0.0;
788  it = constants.begin();
789  for (i = 0; i < 400; i++)
790  {
791  if (it != constants.end())
792  {
793  writeBinary(strm, (char *)&(it->second), sizeof(double));
794  it++;
795  }
796  else
797  {
798  writeBinary(strm, (char *)&z, sizeof(double));
799  }
800  }
801  LOG(DEBUG) << "Pad length 2 = " << (400 - Nconst) * sizeof(double);
802  for (i = 0; i < (400 - Nconst) * sizeof(double); i++)
803  writeBinary(strm, &c, 1);
804 
805  /* data records
806  --------------------------------------------------------- the data,
807  in time order */
808  int nrec = 1;
809  map<double, vector<double>>::iterator jt;
810  for (jt = store.begin(); jt != store.end(); jt++)
811  {
812  LOG(DEBUG) << "writeBinaryFile writes " << nrec << " " << fixed
813  << setprecision(6) << jt->second[0];
814  for (i = 0; i < jt->second.size(); i++)
815  writeBinary(strm, (char *)&jt->second[i], sizeof(double));
816  nrec++;
817  }
818 
819  // TD after writing it out, clear the store array
820  // store.clear();
821 
822  strm.clear();
823  strm.close();
824 
825  return 0;
826  }
827  catch (Exception& e)
828  {
829  GNSSTK_RETHROW(e);
830  }
831  catch (exception& e)
832  {
833  Exception E("std except: " + string(e.what()));
834  GNSSTK_THROW(E);
835  }
836  catch (...)
837  {
838  Exception e("Unknown exception");
839  GNSSTK_THROW(e);
840  }
841  }
842 
843  //---------------------------------------------------------------------------------
844  int SolarSystemEphemeris::readBinaryFile(const string& filename)
845  {
846  try
847  {
848  int iret;
849  readBinaryHeader(filename);
850  iret = readBinaryData(true); // true: store ALL the data in map
851 
852  istrm.clear();
853  istrm.close();
854 
855  return iret;
856  }
857  catch (Exception& e)
858  {
859  GNSSTK_RETHROW(e);
860  }
861  catch (exception& e)
862  {
863  Exception E("std except: " + string(e.what()));
864  GNSSTK_THROW(E);
865  }
866  catch (...)
867  {
868  Exception e("Unknown exception");
869  GNSSTK_THROW(e);
870  }
871  }
872 
873  //---------------------------------------------------------------------------------
874  int SolarSystemEphemeris::initializeWithBinaryFile(const string& filename)
875  {
876  try
877  {
878  int iret;
879 
880  // cout << "Reporting in SolarSystemEphemeris is "
881  // << ConfigureLOG::ToString(ConfigureLOG::ReportingLevel()) << endl;
882 
883  readBinaryHeader(filename);
884  iret = readBinaryData(false); // false: don't store data in map
885  if (iret == 0)
886  {
887  /* EphemerisNumber == -1 : the header has not been read
888  * EphemerisNumber == 0 : the fileposMap has not been read (binary)
889  * EphemerisNumber == constants["DENUM"] :
890  * object has been initialized (binary file),
891  * or header has been read (ASCII file) */
892  EphemerisNumber = int(constants["DENUM"]);
893  LOG(DEBUG) << "initialize sets EphemerisNumber " << EphemerisNumber;
894  }
895 
896  return iret;
897  }
898  catch (Exception& e)
899  {
900  GNSSTK_RETHROW(e);
901  }
902  catch (exception& e)
903  {
904  Exception E("std except: " + string(e.what()));
905  GNSSTK_THROW(E);
906  }
907  catch (...)
908  {
909  Exception e("Unknown exception");
910  GNSSTK_THROW(e);
911  }
912  }
913 
914  //---------------------------------------------------------------------------------
915  // get an inertial position of one body relative to another.
916  void SolarSystemEphemeris::relativeInertialPositionVelocity(
917  double MJD, SolarSystemEphemeris::Planet target,
918  SolarSystemEphemeris::Planet center, double pv[6], bool kilometers)
919  {
920  try
921  {
922  int iret, i;
923 
924  // initialize
925  for (i = 0; i < 6; i++)
926  pv[i] = 0.0;
927 
928  // trivial; return
929  if (target == center)
930  {
931  return;
932  }
933 
934  // get the right record from the file
935  double JD(MJD + MJD_TO_JD);
936  iret = seekToJD(JD);
937  /* -1 out of range : input time is before the first time in file
938  -2 out of range : input time is after the last time in file, or in gap
939  -3 stream is not open or not good, or EOF was found prematurely
940  -4 EphemerisNumber is not defined */
941  if (iret)
942  {
943  if (iret == retEarly || iret == retLate)
944  {
945  Exception e(string("Requested time is ") +
946  (iret == retEarly ? string("before") : string("after")) +
947  string(" the range spanned by the ephemeris."));
948  GNSSTK_THROW(e);
949  }
950  else if (iret == retStrm)
951  {
952  Exception e(string("Stream error on ephemeris binary file"));
953  GNSSTK_THROW(e);
954  }
955  else if (iret == retEphN)
956  {
957  Exception e(string("Ephemeris not initialized"));
958  GNSSTK_THROW(e);
959  }
960  else
961  {
962  Exception e(string("Unknown error on ephemeris binary file"));
963  GNSSTK_THROW(e);
964  }
965  }
966 
967  // compute Nutations or Librations
968  if (target == idNutations || target == idLibrations)
969  {
970  inertialPositionVelocity(
971  MJD, target == idNutations ? NUTATIONS : LIBRATIONS, pv);
972  return;
973  }
974 
975  // define computeID's for target and center
976  computeID TARGET, CENTER;
977 
978  if (target <= idSun)
979  {
980  TARGET = computeID(target - 1);
981  }
982  else if (target == idSolarSystemBarycenter)
983  {
984  TARGET = NONE;
985  }
986  else if (target == idEarthMoonBarycenter)
987  {
988  TARGET = EMBARY;
989  }
990  // (Nutations and Librations are done above)
991  if (center <= idSun)
992  {
993  CENTER = computeID(center - 1);
994  }
995  else if (center == idSolarSystemBarycenter)
996  {
997  CENTER = NONE;
998  }
999  else if (center == idEarthMoonBarycenter)
1000  {
1001  CENTER = EMBARY;
1002  }
1003 
1004  /* Earth and Moon need special treatment - get moon and Earth-moon
1005  barycenter */
1006  double pvmoon[6], pvembary[6], Eratio, Mratio;
1007 
1008  // special cases of Earth AND Moon: Moon result is always geocentric
1009  if (target == idEarth && center == idMoon)
1010  {
1011  TARGET = NONE;
1012  }
1013  if (center == idEarth && target == idMoon)
1014  {
1015  CENTER = NONE;
1016  }
1017 
1018  // special cases of Earth OR Moon, but not both:
1019  if ((target == idEarth && center != idMoon) ||
1020  (center == idEarth && target != idMoon))
1021  {
1022  Eratio = 1.0 / (1.0 + constants["EMRAT"]);
1023  inertialPositionVelocity(MJD, MOON, pvmoon);
1024  }
1025  if ((target == idMoon && center != idEarth) ||
1026  (center == idMoon && target != idEarth))
1027  {
1028  Mratio = constants["EMRAT"] / (1.0 + constants["EMRAT"]);
1029  inertialPositionVelocity(MJD, EMBARY, pvembary);
1030  }
1031 
1032  // compute states for target and center
1033  double pvtarget[6], pvcenter[6];
1034  inertialPositionVelocity(MJD, TARGET, pvtarget);
1035  inertialPositionVelocity(MJD, CENTER, pvcenter);
1036 
1037  /* handle the Earth/Moon special cases
1038  convert from E-M barycenter to Earth */
1039  if (target == idEarth && center != idMoon)
1040  {
1041  for (i = 0; i < 6; i++)
1042  {
1043  pvtarget[i] -= pvmoon[i] * Eratio;
1044  }
1045  }
1046  if (center == idEarth && target != idMoon)
1047  {
1048  for (i = 0; i < 6; i++)
1049  {
1050  pvcenter[i] -= pvmoon[i] * Eratio;
1051  }
1052  }
1053 
1054  if (target == idMoon && center != idEarth)
1055  {
1056  for (i = 0; i < 6; i++)
1057  {
1058  pvtarget[i] = pvembary[i] + pvtarget[i] * Mratio;
1059  }
1060  }
1061  if (center == idMoon && target != idEarth)
1062  {
1063  for (i = 0; i < 6; i++)
1064  {
1065  pvcenter[i] = pvembary[i] + pvcenter[i] * Mratio;
1066  }
1067  }
1068 
1069  // final relative result
1070  for (i = 0; i < 6; i++)
1071  {
1072  pv[i] = pvtarget[i] - pvcenter[i];
1073  }
1074 
1075  if (!kilometers)
1076  {
1077  double AU = constants["AU"];
1078  for (i = 0; i < 6; i++)
1079  pv[i] /= AU;
1080  }
1081  }
1082  catch (Exception& e)
1083  {
1084  GNSSTK_RETHROW(e);
1085  }
1086  catch (exception& e)
1087  {
1088  Exception E("std except: " + string(e.what()));
1089  GNSSTK_THROW(E);
1090  }
1091  catch (...)
1092  {
1093  Exception e("Unknown exception");
1094  GNSSTK_THROW(e);
1095  }
1096  }
1097 
1098  //---------------------------------------------------------------------------------
1099  //---------------------------------------------------------------------------------
1100  // private
1101  void SolarSystemEphemeris::writeBinary(ofstream& strm, const char *ptr,
1102  size_t size)
1103  {
1104  try
1105  {
1106  strm.write(ptr, size);
1107  if (!strm.good())
1108  {
1109  Exception e("Stream error");
1110  GNSSTK_THROW(e);
1111  }
1112  }
1113  catch (Exception& e)
1114  {
1115  GNSSTK_RETHROW(e);
1116  }
1117  catch (exception& e)
1118  {
1119  Exception E("std except: " + string(e.what()));
1120  GNSSTK_THROW(E);
1121  }
1122  catch (...)
1123  {
1124  Exception e("Unknown exception");
1125  GNSSTK_THROW(e);
1126  }
1127  }
1128 
1129  //---------------------------------------------------------------------------------
1130  // private
1131  void SolarSystemEphemeris::readBinary(char *ptr, size_t size)
1132  {
1133  try
1134  {
1135  istrm.read(ptr, size);
1136  if (istrm.eof() || !istrm.good())
1137  {
1138  Exception e("Stream error or premature EOF");
1139  GNSSTK_THROW(e);
1140  }
1141  }
1142  catch (Exception& e)
1143  {
1144  GNSSTK_RETHROW(e);
1145  }
1146  catch (exception& e)
1147  {
1148  Exception E("std except: " + string(e.what()));
1149  GNSSTK_THROW(E);
1150  }
1151  catch (...)
1152  {
1153  Exception e("Unknown exception");
1154  GNSSTK_THROW(e);
1155  }
1156  }
1157 
1158  //---------------------------------------------------------------------------------
1159  // private
1160  void SolarSystemEphemeris::readBinaryHeader(const string& filename)
1161  {
1162  try
1163  {
1164  int i, DENUM, recLength;
1165  char buffer[100];
1166  double AU, EMRAT;
1167  string word;
1168 
1169  // open the input binary file
1170  istrm.open(filename.c_str(), ios::in | ios::binary);
1171  if (!istrm.is_open())
1172  {
1173  Exception e("Failed to open input binary file " + filename +
1174  ". Abort.");
1175  GNSSTK_THROW(e);
1176  }
1177 
1178  // initialize
1179  EphemerisNumber = -1;
1180  constants.clear();
1181  store.clear();
1182  recLength = 0;
1183 
1184  /* ----------------------------------------------------------------
1185  read the first header record
1186  1. 3 labels each of length 84 */
1187  for (i = 0; i < 3; i++)
1188  {
1189  readBinary(buffer, 84); // istrm.read(buffer,84);
1190  recLength += 84;
1191  buffer[84] = '\0';
1192  label[i] = stripTrailing(stripLeading(buffer, " "), " ");
1193  LOG(DEBUG) << "readBinaryHeader reads label " << label[i];
1194  }
1195 
1196  // 2. 400 keys from the const array, each of length 6
1197  vector<string> consts_names;
1198  for (i = 0; i < 400; i++)
1199  {
1200  readBinary(buffer, 6);
1201  buffer[6] = '\0';
1202  recLength += 6;
1203  word = stripLeading(string(buffer));
1204  if (!word.empty())
1205  {
1206  consts_names.push_back(word);
1207  LOG(DEBUG) << "readBinaryHeader reads constant label " << word;
1208  }
1209  }
1210  Nconst = consts_names.size();
1211 
1212  // 3. the three times
1213  readBinary((char *)&startJD, sizeof(double));
1214  LOG(DEBUG) << "readBinaryHeader reads start JD " << fixed
1215  << setprecision(2) << startJD;
1216 
1217  readBinary((char *)&endJD, sizeof(double));
1218  LOG(DEBUG) << "readBinaryHeader reads end JD " << fixed
1219  << setprecision(2) << endJD;
1220 
1221  readBinary((char *)&interval, sizeof(double));
1222  LOG(DEBUG) << "readBinaryHeader reads interval " << interval;
1223 
1224  recLength += 3 * sizeof(double);
1225 
1226  // 4. Ncoeff
1227  readBinary((char *)&Ncoeff, sizeof(int));
1228  recLength += sizeof(int);
1229  LOG(DEBUG) << "readBinaryHeader reads number of coefficients "
1230  << Ncoeff;
1231 
1232  // 5. AU and EMRAT
1233  buffer[sizeof(double)] = '\0';
1234  readBinary((char *)&AU, sizeof(double));
1235  recLength += sizeof(double);
1236  LOG(DEBUG) << "readBinaryHeader reads AU " << fixed << setprecision(4)
1237  << AU;
1238 
1239  readBinary((char *)&EMRAT, sizeof(double));
1240  recLength += sizeof(double);
1241  LOG(DEBUG) << "readBinaryHeader reads EMRAT " << EMRAT;
1242 
1243  // 6. c_arrays for the first 12 planets
1244  for (i = 0; i < 12; i++)
1245  {
1246  readBinary((char *)&c_offset[i], sizeof(int));
1247  readBinary((char *)&c_ncoeff[i], sizeof(int));
1248  readBinary((char *)&c_nsets[i], sizeof(int));
1249  recLength += 3 * sizeof(int);
1250  LOG(DEBUG) << "readBinaryHeader reads " << i << " " << c_offset[i]
1251  << " " << c_ncoeff[i] << " " << c_nsets[i];
1252  }
1253 
1254  // 7. DENUM
1255  double denum;
1256  readBinary((char *)&denum, sizeof(double));
1257  recLength += sizeof(double);
1258  LOG(DEBUG) << "readBinaryHeader reads DENUM directly " << denum;
1259 
1260  // 8. c_arrays for libration
1261  readBinary((char *)&c_offset[12], sizeof(int));
1262  readBinary((char *)&c_ncoeff[12], sizeof(int));
1263  readBinary((char *)&c_nsets[12], sizeof(int));
1264  recLength += 3 * sizeof(int);
1265  LOG(DEBUG) << "readBinaryHeader reads " << 12 << " " << c_offset[12]
1266  << " " << c_ncoeff[12] << " " << c_nsets[12];
1267 
1268  /* 9. pad - records are padded to be the same length as the data
1269  records b/c
1270  JPL does it (for Fortran reasons) - not necessary */
1271  LOG(DEBUG) << "Pad length 1 = " << Ncoeff * sizeof(double) - recLength;
1272  for (i = 0; i < Ncoeff * sizeof(double) - recLength; i++)
1273  {
1274  readBinary(buffer, 1);
1275  }
1276 
1277  // ----------------------------------------------------------------
1278  // the second header record: 400 values from the const array
1279  double d;
1280  for (i = 0; i < 400; i++)
1281  {
1282  readBinary((char *)&d, sizeof(double));
1283  if (i < Nconst)
1284  {
1285  constants[stripTrailing(consts_names[i])] = d;
1286  LOG(DEBUG) << "readBinaryHeader reads " << consts_names[i]
1287  << " = " << fixed << setprecision(18) << setw(24)
1288  << d;
1289  }
1290  }
1291  // pad
1292  LOG(DEBUG) << "Pad length 2 = " << (400 - Nconst) * sizeof(double);
1293  for (i = 0; i < (400 - Nconst) * sizeof(double); i++)
1294  {
1295  readBinary(buffer, 1);
1296  }
1297 
1298  // ----------------------------------------------------------------
1299  // test the header
1300  if (denum == constants["DENUM"])
1301  {
1302  LOG(DEBUG) << "DENUM agrees " << denum;
1303 
1304  /* EphemerisNumber == -1 : the header has not been read
1305  * EphemerisNumber == 0 : the fileposMap has not been read (binary)
1306  * EphemerisNumber == constants["DENUM"] :
1307  * object has been initialized (binary file),
1308  * or header has been read (ASCII file) */
1309  EphemerisNumber = 0;
1310 
1311  // clear the data arrays
1312  store.clear();
1313  }
1314  else
1315  {
1316  LOG(WARNING) << "DENUM (" << denum
1317  << ") does not equal the array value ("
1318  << constants["DENUM"] << ")";
1319  }
1320  }
1321  catch (Exception& e)
1322  {
1323  GNSSTK_RETHROW(e);
1324  }
1325  catch (exception& e)
1326  {
1327  Exception E("std except: " + string(e.what()));
1328  GNSSTK_THROW(E);
1329  }
1330  catch (...)
1331  {
1332  Exception e("Unknown exception");
1333  GNSSTK_THROW(e);
1334  }
1335  }
1336 
1337  //---------------------------------------------------------------------------------
1338  /* private
1339  return 0 ok, or -4 EphemerisNumber is not defined */
1340  int SolarSystemEphemeris ::readBinaryData(bool save)
1341  {
1342  try
1343  {
1344  // has the header been read?
1345  if (EphemerisNumber == -1)
1346  {
1347  return retEphN;
1348  }
1349 
1350  // read the data, optionally storing it all; fill the file position map
1351  int iret = retEarly, nrec = 1;
1352  double prev = 0.0;
1353  vector<double> dataVector;
1354  while (!istrm.eof() && istrm.good())
1355  {
1356  long filepos = istrm.tellg();
1357  iret = readBinaryRecord(dataVector);
1358  if (iret == retLate)
1359  {
1360  iret = 0;
1361  break;
1362  } // EOF
1363  if (iret)
1364  {
1365  break;
1366  }
1367 
1368  // if saving all in store, add it here
1369  if (save)
1370  {
1371  store[dataVector[0]] = dataVector;
1372  }
1373 
1374  // put the first record in coefficients array
1375  if (nrec == 1)
1376  {
1377  coefficients = dataVector;
1378  }
1379 
1380  // build the positions map
1381  fileposMap[dataVector[0]] = filepos;
1382 
1383  if (nrec > 1 && dataVector[0] != prev)
1384  {
1385  ostringstream oss;
1386  oss << "ERROR: found gap in data at " << nrec << fixed
1387  << setprecision(6) << " : prev end = " << prev
1388  << " != new beg = " << dataVector[0];
1389  Exception e(oss.str());
1390  GNSSTK_THROW(e);
1391  }
1392 
1393  // remember the end time for next record
1394  prev = dataVector[1];
1395 
1396  // count records
1397  nrec++;
1398  }
1399 
1400  istrm.clear();
1401 
1402  return iret;
1403  }
1404  catch (Exception& e)
1405  {
1406  GNSSTK_RETHROW(e);
1407  }
1408  catch (exception& e)
1409  {
1410  Exception E("std except: " + string(e.what()));
1411  GNSSTK_THROW(E);
1412  }
1413  catch (...)
1414  {
1415  Exception e("Unknown exception");
1416  GNSSTK_THROW(e);
1417  }
1418  }
1419 
1420  //---------------------------------------------------------------------------------
1421  /* private
1422  return 0 ok, or
1423  -3 stream is not open or not good
1424  -4 EphemerisNumber is not defined
1425  For -3,-4 : initializeWithBinaryFile() has not been called, or reading
1426  failed. */
1427  int SolarSystemEphemeris::readBinaryRecord(vector<double>& dataVector)
1428  {
1429  try
1430  {
1431  if (!istrm)
1432  {
1433  return retStrm;
1434  }
1435  if (istrm.eof() || !istrm.good())
1436  {
1437  return retStrm;
1438  }
1439  if (EphemerisNumber <= -1)
1440  {
1441  return retEphN;
1442  }
1443 
1444  dataVector.clear();
1445 
1446  double d;
1447  for (int i = 0; i < Ncoeff; i++)
1448  {
1449  istrm.read((char *)&d,
1450  sizeof(double)); // don't use readBinary(), to catch EOF
1451  if (istrm.eof())
1452  {
1453  return retLate;
1454  }
1455  if (!istrm.good())
1456  {
1457  return retStrm;
1458  }
1459  dataVector.push_back(d);
1460  }
1461 
1462  return 0;
1463  }
1464  catch (Exception& e)
1465  {
1466  GNSSTK_RETHROW(e);
1467  }
1468  catch (exception& e)
1469  {
1470  Exception E("std except: " + string(e.what()));
1471  GNSSTK_THROW(E);
1472  }
1473  catch (...)
1474  {
1475  Exception e("Unknown exception");
1476  GNSSTK_THROW(e);
1477  }
1478  }
1479 
1480  //---------------------------------------------------------------------------------
1481  /* private
1482  return 0 ok, or
1483  -1 out of range : input time is before the first time in file
1484  -2 out of range : input time is after the last time in file, or in a gap
1485  -3 stream is not open or not good, or EOF was found prematurely
1486  -4 EphemerisNumber is not defined
1487  For -3,-4 : initializeWithBinaryFile() has not been called, or reading
1488  failed. */
1489  int SolarSystemEphemeris ::seekToJD(double JD)
1490  {
1491  try
1492  {
1493  if (!istrm)
1494  {
1495  return retStrm;
1496  }
1497  if (istrm.eof() || !istrm.good())
1498  {
1499  return retStrm;
1500  }
1501  if (EphemerisNumber != int(constants["DENUM"]))
1502  {
1503  return retEphN;
1504  }
1505 
1506  if (coefficients[0] <= JD && JD <= coefficients[1])
1507  {
1508  return 0;
1509  }
1510 
1511  map<double, long>::const_iterator it; // key >= input
1512  it = fileposMap.lower_bound(
1513  JD); // it points to first element with JD <= time
1514 
1515  if (it == fileposMap.begin() && JD < it->first)
1516  {
1517  return retEarly; // failure: JD is before the first record
1518  }
1519 
1520  if (it ==
1521  fileposMap.end() // if beyond the found record, go to previous;
1522  || JD < it->first)
1523  {
1524  it--; // but beware the "lower_bound found the =" case
1525  }
1526 
1527  istrm.seekg(it->second, ios_base::beg); // get the record
1528  int iret = readBinaryRecord(coefficients);
1529  if (iret == retLate)
1530  {
1531  iret = retStrm; // this means EOF during data read
1532  }
1533  if (iret)
1534  {
1535  return iret; // reading failed
1536  }
1537 
1538  if (JD > coefficients[1])
1539  {
1540  return retLate; // failure: JD is after the last record, or
1541  }
1542  // JD is in a gap between records
1543  return 0;
1544  }
1545  catch (Exception& e)
1546  {
1547  GNSSTK_RETHROW(e);
1548  }
1549  catch (exception& e)
1550  {
1551  Exception E("std except: " + string(e.what()));
1552  GNSSTK_THROW(E);
1553  }
1554  catch (...)
1555  {
1556  Exception e("Unknown exception");
1557  GNSSTK_THROW(e);
1558  }
1559  }
1560 
1561  //---------------------------------------------------------------------------------
1562  // private
1563  void SolarSystemEphemeris::inertialPositionVelocity(
1564  double MJD, SolarSystemEphemeris::computeID which, double PV[6])
1565  {
1566  try
1567  {
1568  int i, j, i0, ncomp, offset;
1569 
1570  for (i = 0; i < 6; i++)
1571  {
1572  PV[i] = 0.0;
1573  }
1574  if (which == NONE)
1575  {
1576  return;
1577  }
1578 
1579  /* coefficients[0,1] give span of JD's in which coefficients[2,...] are
1580  applicable coefficients[0,1] are even days JDs - 2452xxx.5 =>
1581  secOfDay() for these == 0. */
1582  double T, Tbeg, Tspan, Tspan0;
1583  Tbeg = coefficients[0];
1584  Tspan0 = Tspan = coefficients[1] - coefficients[0];
1585  i0 = c_offset[which] - 1; // index of first coefficient in array
1586  ncomp = (which == NUTATIONS ? 2 : 3); // number of components returned
1587 
1588  // if more than one set, find the right set
1589  if (c_nsets[which] > 1)
1590  {
1591  Tspan /= double(c_nsets[which]);
1592  for (j = c_nsets[which]; j > 0; j--)
1593  {
1594  Tbeg = coefficients[0] + double(j - 1) * Tspan;
1595  if (MJD > Tbeg - MJD_TO_JD)
1596  { // == with j==1 is the default
1597  i0 += (j - 1) * ncomp * c_ncoeff[which];
1598  break;
1599  }
1600  }
1601  }
1602 
1603  // normalized time
1604  T = 2.0 * (MJD - (Tbeg - MJD_TO_JD)) / Tspan - 1.0;
1605 
1606  // interpolate
1607  int N = c_ncoeff[which];
1608  vector<double> C(N, 0.0); // Chebyshev
1609  vector<double> U(N, 0.0); // derivative of Chebyshev
1610  for (i = 0; i < ncomp; i++)
1611  { // loop over components
1612 
1613  // seed the Chebyshev recursions
1614  C[0] = 1;
1615  C[1] = T; // C[2] = 2*T*T-1;
1616  U[0] = 0;
1617  U[1] = 1; // U[2] = 4*T;
1618 
1619  // generate the Chebyshevs
1620  for (j = 2; j < N; j++)
1621  {
1622  C[j] = 2 * T * C[j - 1] - C[j - 2];
1623  U[j] = 2 * T * U[j - 1] + 2 * C[j - 1] - U[j - 2];
1624  }
1625 
1626  // compute P and V
1627  // done above PV[i] = PV[i+3] = 0.0;
1628  for (j = N - 1; j > -1; j--) // POS
1629  {
1630  PV[i] += coefficients[i0 + j + i * N] * C[j];
1631  }
1632  for (j = N - 1; j > 0; j--) // j>0 b/c U[0]=0 // VEL
1633  {
1634  PV[i + ncomp] += coefficients[i0 + j + i * N] * U[j];
1635  }
1636 
1637  // convert velocity to 'per day'
1638  PV[i + ncomp] *= 2 * double(c_nsets[which]) / Tspan0;
1639  }
1640  }
1641  catch (Exception& e)
1642  {
1643  GNSSTK_RETHROW(e);
1644  }
1645  catch (exception& e)
1646  {
1647  Exception E("std except: " + string(e.what()));
1648  GNSSTK_THROW(E);
1649  }
1650  catch (...)
1651  {
1652  Exception e("Unknown exception");
1653  GNSSTK_THROW(e);
1654  }
1655  }
1656 
1657  //---------------------------------------------------------------------------------
1658 } // end namespace gnsstk
Ncoeff
const int Ncoeff
Number of terms in the series.
Definition: IERS1996NutationData.hpp:169
gnsstk::StringUtils::asInt
long asInt(const std::string &s)
Definition: StringUtils.hpp:713
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::SSEDouble::operator=
SSEDouble & operator=(const string &s)
Assign a value by decoding a string using existing formatting.
Definition: SolarSystemEphemeris.cpp:78
StringUtils.hpp
gnsstk::StringUtils::center
std::string & center(std::string &s, const std::string::size_type length, const char pad=' ')
Definition: StringUtils.hpp:1607
gnsstk::SolarSystemEphemeris::computeID
computeID
Definition: SolarSystemEphemeris.hpp:412
gnsstk::SSEDouble::SSEDouble
SSEDouble(const string &str)
Decode a string.
Definition: SolarSystemEphemeris.cpp:75
logstream.hpp
gnsstk::StringUtils::asString
std::string asString(IonexStoreStrategy e)
Convert a IonexStoreStrategy to a whitespace-free string name.
Definition: IonexStoreStrategy.cpp:46
gnsstk::SSEDouble::SSEDouble
SSEDouble(double d=0)
Definition: SolarSystemEphemeris.cpp:68
gnsstk::WARNING
@ WARNING
Definition: logstream.hpp:57
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
gnsstk::StringUtils::FFAlign::Right
@ Right
Formatted output will be right-aligned.
gnsstk::MJD_TO_JD
const double MJD_TO_JD
Add this offset to convert Modified Julian Date to Julian Date.
Definition: TimeConstants.hpp:51
gnsstk::StringUtils::FFLead
FFLead
Leading character for floatFormat(), after any whitespace or sign.
Definition: StringUtils.hpp:109
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::FormattedDouble
Definition: FormattedDouble.hpp:70
gnsstk::StringUtils::FFSign::NegOnly
@ NegOnly
Prefix output with a minus sign (neg) or nothing (pos)
SolarSystemEphemeris.hpp
gnsstk::convertJDtoCalendar
void convertJDtoCalendar(long jd, int &iyear, int &imonth, int &iday)
Definition: TimeConverters.cpp:52
gnsstk::SolarSystemEphemeris::Planet
Planet
These are indexes used by the caller of inertialPositionVelocity().
Definition: SolarSystemEphemeris.hpp:110
gnsstk::StringUtils::FFLead::Decimal
@ Decimal
Start with decimal, e.g. .12345.
gnsstk::StringUtils::stripFirstWord
std::string stripFirstWord(std::string &s, const char delimiter=' ')
Definition: StringUtils.hpp:2253
gnsstk::VERBOSE
@ VERBOSE
Definition: logstream.hpp:57
TimeConverters.hpp
FormattedDouble.hpp
GNSSTK_RETHROW
#define GNSSTK_RETHROW(exc)
Definition: Exception.hpp:369
example6.interval
interval
Definition: example6.py:75
gnsstk::StringUtils
Definition: IonexStoreStrategy.cpp:44
gnsstk::StringUtils::rightJustify
std::string & rightJustify(std::string &s, const std::string::size_type length, const char pad=' ')
Definition: StringUtils.hpp:1557
LOG
#define LOG(level)
define the macro that is used to write to the log stream
Definition: logstream.hpp:315
std
Definition: Angle.hpp:142
gnsstk::INFO
@ INFO
Definition: logstream.hpp:57
gnsstk::StringUtils::FFSign
FFSign
How to handle sign in floatFormat()
Definition: StringUtils.hpp:117
gnsstk::MJD
Definition: MJD.hpp:54
gnsstk::StringUtils::leftJustify
std::string & leftJustify(std::string &s, const std::string::size_type length, const char pad=' ')
Definition: StringUtils.hpp:1582
GNSSTK_THROW
#define GNSSTK_THROW(exc)
Definition: Exception.hpp:366
gnsstk::DEBUG
@ DEBUG
Definition: logstream.hpp:57
gnsstk::StringUtils::FFAlign
FFAlign
Alignment of data for floatFormat()
Definition: StringUtils.hpp:125
gnsstk::SSEDouble
Class for the format used in this code.
Definition: SolarSystemEphemeris.cpp:65

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