PCL
JPLEphemeris.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.1.19
6 // ----------------------------------------------------------------------------
7 // pcl/JPLEphemeris.h - Released 2019-11-07T10:59:34Z
8 // ----------------------------------------------------------------------------
9 // This file is part of the PixInsight Class Library (PCL).
10 // PCL is a multiplatform C++ framework for development of PixInsight modules.
11 //
12 // Copyright (c) 2003-2019 Pleiades Astrophoto S.L. All Rights Reserved.
13 //
14 // Redistribution and use in both source and binary forms, with or without
15 // modification, is permitted provided that the following conditions are met:
16 //
17 // 1. All redistributions of source code must retain the above copyright
18 // notice, this list of conditions and the following disclaimer.
19 //
20 // 2. All redistributions in binary form must reproduce the above copyright
21 // notice, this list of conditions and the following disclaimer in the
22 // documentation and/or other materials provided with the distribution.
23 //
24 // 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
25 // of their contributors, may be used to endorse or promote products derived
26 // from this software without specific prior written permission. For written
27 // permission, please contact info@pixinsight.com.
28 //
29 // 4. All products derived from this software, in any form whatsoever, must
30 // reproduce the following acknowledgment in the end-user documentation
31 // and/or other materials provided with the product:
32 //
33 // "This product is based on software from the PixInsight project, developed
34 // by Pleiades Astrophoto and its contributors (http://pixinsight.com/)."
35 //
36 // Alternatively, if that is where third-party acknowledgments normally
37 // appear, this acknowledgment must be reproduced in the product itself.
38 //
39 // THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
40 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
41 // TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
42 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
43 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
44 // EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
45 // INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
46 // DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
47 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
48 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
49 // POSSIBILITY OF SUCH DAMAGE.
50 // ----------------------------------------------------------------------------
51 
52 #ifndef __PCL_JPLEphemeris_h
53 #define __PCL_JPLEphemeris_h
54 
56 
57 #include <pcl/Defs.h>
58 
59 #include <pcl/Arguments.h>
60 #include <pcl/EphemerisFile.h>
61 #include <pcl/ErrorHandler.h>
62 #include <pcl/KeyValue.h>
63 #include <pcl/MultiVector.h>
64 #include <pcl/StringList.h>
65 #include <pcl/TimePoint.h>
66 
67 namespace pcl
68 {
69 
70 // ----------------------------------------------------------------------------
71 
108 namespace JPLEphemerisItem
109 {
110  enum value_type
111  {
112  Unknown = -1,
113 
114  // Native JPL DE/LE items, as of 438t
115  Mercury = 0,
116  Venus,
117  EarthMoonBarycenter,
118  Mars,
119  Jupiter,
120  Saturn,
121  Uranus,
122  Neptune,
123  Pluto,
124  Moon,
125  Sun,
126  Nutations,
127  LunarLibration,
128  LunarMantleVelocity,
129  TT_TDB,
130 
131  NumberOfNativeItems,
132  NumberOfRequiredItems = Sun + 1,
133 
134  // Special items synthesized from EMB and Moon
135  Earth = 100, // Geocenter with respect to solar system barycenter
136  SSBMoon = 101 // Moon with respect to solar system barycenter
137  };
138 }
139 
140 // ----------------------------------------------------------------------------
141 
161 {
162 public:
163 
168  typedef JPLEphemerisItem::value_type item_index;
169 
175 
180 
188  JPLEphemeris( const String& filePath );
189 
210  void AddData( const String& filePath );
211 
212 #ifndef __PCL_NO_JPL_EPHEMERIS_TEST_ROUTINES
213 
236  bool Test( const String& filePath, bool verbose = true ) const;
237 
238 #endif // !__PCL_NO_JPL_EPHEMERIS_TEST_ROUTINES
239 
245  int DENumber() const
246  {
247  return TruncInt( ConstantValue( "DENUM" ) );
248  }
249 
255  int LENumber() const
256  {
257  return TruncInt( ConstantValue( "LENUM" ) );
258  }
259 
269  double EphemerisStartJD() const
270  {
271  return m_ephStartJD;
272  }
273 
283  double EphemerisEndJD() const
284  {
285  return m_ephEndJD;
286  }
287 
296  double StartJD() const
297  {
298  return m_startJD;
299  }
300 
309  double EndJD() const
310  {
311  return m_endJD;
312  }
313 
321  const constant_list& Constants() const
322  {
323  return m_constants;
324  }
325 
332  double ConstantValue( const IsoString& name ) const
333  {
334  constant_list::const_iterator i =
335  BinarySearch( m_constants.Begin(), m_constants.End(), EphemerisConstant( name ) );
336  if ( i == m_constants.End() )
337  throw Error( "Undefined integration constant '" + name + '\'' );
338  return i->value;
339  }
340 
346  bool IsItemAvailable( int i ) const
347  {
348  if ( i == JPLEphemerisItem::Earth || i == JPLEphemerisItem::SSBMoon )
349  return IsItemAvailable( JPLEphemerisItem::EarthMoonBarycenter ) && IsItemAvailable( JPLEphemerisItem::Moon );
350  return i >= 0 && i < JPLEphemerisItem::NumberOfNativeItems && m_blockIndex[i].subblocks > 0;
351  }
352 
366  static int ComponentsForItem( int i )
367  {
368  return (i == JPLEphemerisItem::Nutations) ? 2 : (i == JPLEphemerisItem::TT_TDB ? 1 : 3);
369  }
370 
389  void ComputeState( Vector& r, double jd1, double jd2, item_index i ) const
390  {
391  if ( i == JPLEphemerisItem::Earth )
392  {
393  Vector rm( 3 );
394  Interpolate( r.Begin(), nullptr, jd1, jd2, JPLEphemerisItem::EarthMoonBarycenter );
395  Interpolate( rm.Begin(), nullptr, jd1, jd2, JPLEphemerisItem::Moon );
396  rm *= m_emb2Earth;
397  r -= rm;
398  }
399  else if ( i == JPLEphemerisItem::SSBMoon )
400  {
401  Vector re( 3 );
402  Interpolate( re.Begin(), nullptr, jd1, jd2, JPLEphemerisItem::EarthMoonBarycenter );
403  Interpolate( r.Begin(), nullptr, jd1, jd2, JPLEphemerisItem::Moon );
404  r += re - m_emb2Earth*r;
405  }
406  else
407  Interpolate( r.Begin(), nullptr, jd1, jd2, i );
408 
409  if ( i != JPLEphemerisItem::Moon )
410  if ( i != JPLEphemerisItem::LunarLibration )
411  if ( i != JPLEphemerisItem::Nutations )
412  r *= m_km2au;
413  }
414 
438  void ComputeState( Vector& r, Vector& v, double jd1, double jd2, item_index i ) const
439  {
440  if ( i == JPLEphemerisItem::Earth )
441  {
442  Vector rm( 3 ), vm( 3 );
443  Interpolate( r.Begin(), v.Begin(), jd1, jd2, JPLEphemerisItem::EarthMoonBarycenter );
444  Interpolate( rm.Begin(), vm.Begin(), jd1, jd2, JPLEphemerisItem::Moon );
445  rm *= m_emb2Earth;
446  vm *= m_emb2Earth;
447  r -= rm;
448  v -= vm;
449  }
450  else if ( i == JPLEphemerisItem::SSBMoon )
451  {
452  Vector re( 3 ), ve( 3 );
453  Interpolate( re.Begin(), ve.Begin(), jd1, jd2, JPLEphemerisItem::EarthMoonBarycenter );
454  Interpolate( r.Begin(), v.Begin(), jd1, jd2, JPLEphemerisItem::Moon );
455  r += re - m_emb2Earth*r;
456  v += ve - m_emb2Earth*v;
457  }
458  else
459  Interpolate( r.Begin(), v.Begin(), jd1, jd2, i );
460 
461  if ( i != JPLEphemerisItem::Moon )
462  if ( i != JPLEphemerisItem::LunarLibration )
463  if ( i != JPLEphemerisItem::Nutations )
464  {
465  r *= m_km2au;
466  v *= m_km2au;
467  }
468  }
469 
477  void ComputeState( Vector& r, TimePoint t, item_index i ) const
478  {
479  ComputeState( r, t.JDI(), t.JDF(), i );
480  }
481 
490  void ComputeState( Vector& r, Vector& v, TimePoint t, item_index i ) const
491  {
492  ComputeState( r, v, t.JDI(), t.JDF(), i );
493  }
494 
503  IsoString Summary() const;
504 
505 private:
506 
507  struct block_index_item
508  {
509  int offset = 0; // offset of this ephemeris item in a block
510  int coefficients = 0; // number of coefficients per subblock
511  int subblocks = 0; // number of subblocks
512  };
513 
515 
516  typedef MultiVector block_list;
517 
518  double m_ephStartJD; // ephemeris start time
519  double m_ephEndJD; // ephemeris end time
520  int m_blockDelta; // block time span in days
521  constant_list m_constants; // integration constants
522  block_index m_blockIndex; // coefficients block index
523  block_list m_blocks; // coefficients blocks
524  double m_startJD; // start time of the first coefficients block
525  double m_endJD; // end time of the last coefficients block
526  double m_km2au; // astronomical unit in kilometers
527  double m_emb2Earth; // scale factor for translation Earth-Moon barycenter -> Earth
528 
534  void Interpolate( double* r, double* v, double jd1, double jd2, item_index n ) const
535  {
536  // block index
537  int ix = TruncInt( ((jd1 - m_startJD) + jd2)/m_blockDelta );
538  if ( ix < 0 || ix > int( m_blocks.Length() ) )
539  throw Error( String().Format( "Time point out of range: %.15g", jd1+jd2 ) );
540  if ( ix == int( m_blocks.Length() ) )
541  --ix;
542 
543  // block coefficients
544  const double* block = m_blocks[ix].Begin();
545 
546  // block index item
547  const block_index_item& index = m_blockIndex[n];
548 
549  // block starting time
550  double jd0 = *block;
551 
552  // number of components
553  int nv = ComponentsForItem( n );
554 
555  // number of coefficients per subblock
556  int nk = index.coefficients;
557 
558  // subblock time span in days
559  int dx = m_blockDelta/index.subblocks;
560 
561  // subblock index
562  ix = ((jd1 - jd0) + jd2)/dx;
563  if ( ix == index.subblocks )
564  --ix;
565 
566  // subblock coefficients
567  const double* k = block + index.offset + ix*nv*nk;
568 
569  // compute state variables
570  Vector pc( nk );
571  pc[0] = 1;
572  pc[1] = (2*((jd1 - (jd0 + ix*dx)) + jd2) - dx)/dx;
573  double y2 = 2*pc[1];
574  for ( int i = 2; i < nk; ++i )
575  pc[i] = y2*pc[i-1] - pc[i-2];
576  for ( int i = 0; i < nv; ++i )
577  {
578  r[i] = 0;
579  for ( int j = 0; j < nk; ++j )
580  r[i] += pc[j] * *k++;
581  }
582 
583  // if requested, compute first derivatives.
584  if ( v != nullptr )
585  {
586  double vf = 2.0*index.subblocks/m_blockDelta;
587  Vector vc( nk );
588  vc[0] = 0;
589  vc[1] = 1;
590  vc[2] = 2*y2;
591  for ( int i = 3; i < nk; ++i )
592  vc[i] = y2*vc[i-1] + 2*pc[i-1] - vc[i-2];
593  for ( int i = nv; --i >= 0; )
594  {
595  v[i] = 0;
596  for ( int j = nk; --j >= 0; )
597  v[i] += vc[j] * *--k;
598  v[i] *= vf;
599  }
600  }
601  }
602 };
603 
604 // ----------------------------------------------------------------------------
605 
606 } // pcl
607 
608 #endif // __PCL_JPLEphemeris_h
609 
610 // ----------------------------------------------------------------------------
611 // EOF pcl/JPLEphemeris.h - Released 2019-11-07T10:59:34Z
An instant in any timescale.
Definition: TimePoint.h:102
constexpr int JDI() const
Definition: TimePoint.h:556
bool IsItemAvailable(int i) const
Definition: JPLEphemeris.h:346
constexpr double JDF() const
Definition: TimePoint.h:565
double EphemerisStartJD() const
Definition: JPLEphemeris.h:269
JPLEphemerisItem::value_type item_index
Definition: JPLEphemeris.h:168
void ComputeState(Vector &r, TimePoint t, item_index i) const
Definition: JPLEphemeris.h:477
PCL root namespace.
Definition: AbstractImage.h:76
EphemerisConstant constant
Definition: JPLEphemeris.h:174
int LENumber() const
Definition: JPLEphemeris.h:255
double StartJD() const
Definition: JPLEphemeris.h:296
int DENumber() const
Definition: JPLEphemeris.h:245
static int ComponentsForItem(int i)
Definition: JPLEphemeris.h:366
void ComputeState(Vector &r, double jd1, double jd2, item_index i) const
Definition: JPLEphemeris.h:389
const constant_list & Constants() const
Definition: JPLEphemeris.h:321
64-bit floating point real multivector.
FI BinarySearch(FI i, FI j, const T &v)
Definition: Search.h:170
Unicode (UTF-16) string.
Definition: String.h:7911
int TruncInt(T x)
Definition: Math.h:1035
void ComputeState(Vector &r, Vector &v, TimePoint t, item_index i) const
Definition: JPLEphemeris.h:490
Dynamic list of ephemeris numerical constants.
A simple exception with an associated error message.
Definition: Exception.h:213
JPL planetary ephemeris.
Definition: JPLEphemeris.h:160
64-bit floating point real vector.
double EndJD() const
Definition: JPLEphemeris.h:309
void ComputeState(Vector &r, Vector &v, double jd1, double jd2, item_index i) const
Definition: JPLEphemeris.h:438
A numerical constant defined in an ephemeris file (XEPH format).
double EphemerisEndJD() const
Definition: JPLEphemeris.h:283
EphemerisConstantList constant_list
Definition: JPLEphemeris.h:179
double ConstantValue(const IsoString &name) const
Definition: JPLEphemeris.h:332
Eight-bit string (ISO/IEC-8859-1 or UTF-8 string)
Definition: String.h:5387