JPLEphemeris.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.7.0
// ----------------------------------------------------------------------------
// pcl/JPLEphemeris.h - Released 2024-06-18T15:48:54Z
// ----------------------------------------------------------------------------
// This file is part of the PixInsight Class Library (PCL).
// PCL is a multiplatform C++ framework for development of PixInsight modules.
//
// Copyright (c) 2003-2024 Pleiades Astrophoto S.L. All Rights Reserved.
//
// Redistribution and use in both source and binary forms, with or without
// modification, is permitted provided that the following conditions are met:
//
// 1. All redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//
// 2. All redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//
// 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
//    of their contributors, may be used to endorse or promote products derived
//    from this software without specific prior written permission. For written
//    permission, please contact info@pixinsight.com.
//
// 4. All products derived from this software, in any form whatsoever, must
//    reproduce the following acknowledgment in the end-user documentation
//    and/or other materials provided with the product:
//
//    "This product is based on software from the PixInsight project, developed
//    by Pleiades Astrophoto and its contributors (https://pixinsight.com/)."
//
//    Alternatively, if that is where third-party acknowledgments normally
//    appear, this acknowledgment must be reproduced in the product itself.
//
// THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
// INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
// DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
// ----------------------------------------------------------------------------
 
#ifndef __PCL_JPLEphemeris_h
#define __PCL_JPLEphemeris_h
 
 
#include <pcl/Defs.h>
 
#include <pcl/Arguments.h>
#include <pcl/EphemerisFile.h>
#include <pcl/ErrorHandler.h>
#include <pcl/KeyValue.h>
#include <pcl/MultiVector.h>
#include <pcl/StringList.h>
#include <pcl/TimePoint.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
namespace JPLEphemerisItem
{
   enum value_type
   {
      Unknown = -1,
 
      // Native JPL DE/LE items, as of 438t
      Mercury = 0,
      Venus,
      EarthMoonBarycenter,
      Mars,
      Jupiter,
      Saturn,
      Uranus,
      Neptune,
      Pluto,
      Moon,
      Sun,
      Nutations,
      LunarLibration,
      LunarMantleVelocity,
      TT_TDB,
 
      NumberOfNativeItems,
      NumberOfRequiredItems = Sun + 1,
 
      // Special items synthesized from EMB and Moon
      Earth = 100,  // Geocenter with respect to solar system barycenter
      SSBMoon = 101 // Moon with respect to solar system barycenter
   };
}
 
// ----------------------------------------------------------------------------
 
class JPLEphemeris
{
public:
 
   using item_index = JPLEphemerisItem::value_type;
 
   using constant = EphemerisConstant;
 
   using constant_list = EphemerisConstantList;
 
   JPLEphemeris( const String& filePath );
 
   void AddData( const String& filePath );
 
#ifndef __PCL_NO_JPL_EPHEMERIS_TEST_ROUTINES
 
   bool Test( const String& filePath, bool verbose = true ) const;
 
#endif   // !__PCL_NO_JPL_EPHEMERIS_TEST_ROUTINES
 
   int DENumber() const
   {
      return TruncInt( ConstantValue( "DENUM" ) );
   }
 
   int LENumber() const
   {
      return TruncInt( ConstantValue( "LENUM" ) );
   }
 
   double EphemerisStartJD() const
   {
      return m_ephStartJD;
   }
 
   double EphemerisEndJD() const
   {
      return m_ephEndJD;
   }
 
   double StartJD() const
   {
      return m_startJD;
   }
 
   double EndJD() const
   {
      return m_endJD;
   }
 
   const constant_list& Constants() const
   {
      return m_constants;
   }
 
   double ConstantValue( const IsoString& name ) const
   {
      constant_list::const_iterator i =
      BinarySearch( m_constants.Begin(), m_constants.End(), EphemerisConstant( name ) );
      if ( i == m_constants.End() )
         throw Error( "Undefined integration constant '" + name + '\'' );
      return i->value;
   }
 
   bool IsItemAvailable( int i ) const
   {
      if ( i == JPLEphemerisItem::Earth || i == JPLEphemerisItem::SSBMoon )
         return IsItemAvailable( JPLEphemerisItem::EarthMoonBarycenter ) && IsItemAvailable( JPLEphemerisItem::Moon );
      return i >= 0 && i < JPLEphemerisItem::NumberOfNativeItems && m_blockIndex[i].subblocks > 0;
   }
 
   static int ComponentsForItem( int i )
   {
      return (i == JPLEphemerisItem::Nutations) ? 2 : (i == JPLEphemerisItem::TT_TDB ? 1 : 3);
   }
 
   void ComputeState( Vector& r, double jd1, double jd2, item_index i ) const
   {
      if ( i == JPLEphemerisItem::Earth )
      {
         Vector rm( 3 );
         Interpolate( r.Begin(), nullptr, jd1, jd2, JPLEphemerisItem::EarthMoonBarycenter );
         Interpolate( rm.Begin(), nullptr, jd1, jd2, JPLEphemerisItem::Moon );
         rm *= m_emb2Earth;
         r -= rm;
      }
      else if ( i == JPLEphemerisItem::SSBMoon )
      {
         Vector re( 3 );
         Interpolate( re.Begin(), nullptr, jd1, jd2, JPLEphemerisItem::EarthMoonBarycenter );
         Interpolate( r.Begin(), nullptr, jd1, jd2, JPLEphemerisItem::Moon );
         r += re - m_emb2Earth*r;
      }
      else
         Interpolate( r.Begin(), nullptr, jd1, jd2, i );
 
      if ( i != JPLEphemerisItem::Moon )
         if ( i != JPLEphemerisItem::LunarLibration )
            if ( i != JPLEphemerisItem::Nutations )
               r *= m_km2au;
   }
 
   void ComputeState( Vector& r, Vector& v, double jd1, double jd2, item_index i ) const
   {
      if ( i == JPLEphemerisItem::Earth )
      {
         Vector rm( 3 ), vm( 3 );
         Interpolate( r.Begin(), v.Begin(), jd1, jd2, JPLEphemerisItem::EarthMoonBarycenter );
         Interpolate( rm.Begin(), vm.Begin(), jd1, jd2, JPLEphemerisItem::Moon );
         rm *= m_emb2Earth;
         vm *= m_emb2Earth;
         r -= rm;
         v -= vm;
      }
      else if ( i == JPLEphemerisItem::SSBMoon )
      {
         Vector re( 3 ), ve( 3 );
         Interpolate( re.Begin(), ve.Begin(), jd1, jd2, JPLEphemerisItem::EarthMoonBarycenter );
         Interpolate( r.Begin(), v.Begin(), jd1, jd2, JPLEphemerisItem::Moon );
         r += re - m_emb2Earth*r;
         v += ve - m_emb2Earth*v;
      }
      else
         Interpolate( r.Begin(), v.Begin(), jd1, jd2, i );
 
      if ( i != JPLEphemerisItem::Moon )
         if ( i != JPLEphemerisItem::LunarLibration )
            if ( i != JPLEphemerisItem::Nutations )
            {
               r *= m_km2au;
               v *= m_km2au;
            }
   }
 
   void ComputeState( Vector& r, TimePoint t, item_index i ) const
   {
      ComputeState( r, t.JDI(), t.JDF(), i );
   }
 
   void ComputeState( Vector& r, Vector& v, TimePoint t, item_index i ) const
   {
      ComputeState( r, v, t.JDI(), t.JDF(), i );
   }
 
   IsoString Summary() const;
 
private:
 
   struct block_index_item
   {
      int offset = 0;            // offset of this ephemeris item in a block
      int coefficients = 0;      // number of coefficients per subblock
      int subblocks = 0;         // number of subblocks
   };
 
   using block_index = Array<block_index_item>;
   using block_list = MultiVector;
 
   double        m_ephStartJD = 0;  // ephemeris start time
   double        m_ephEndJD = 0;    // ephemeris end time
   int           m_blockDelta = 0;  // block time span in days
   constant_list m_constants;       // integration constants
   block_index   m_blockIndex;      // coefficients block index
   block_list    m_blocks;          // coefficients blocks
   double        m_startJD = 0;     // start time of the first coefficients block
   double        m_endJD = 0;       // end time of the last coefficients block
   double        m_km2au = 0;       // astronomical unit in kilometers
   double        m_emb2Earth = 0;   // scale factor for translation Earth-Moon barycenter -> Earth
 
   void Interpolate( double* r, double* v, double jd1, double jd2, item_index n ) const
   {
      // block index
      int ix = TruncInt( ((jd1 - m_startJD) + jd2)/m_blockDelta );
      if ( ix < 0 || ix > int( m_blocks.Length() ) )
         throw Error( String().Format( "Time point out of range: %.15g", jd1+jd2 ) );
      if ( ix == int( m_blocks.Length() ) )
         --ix;
 
      // block coefficients
      const double* block = m_blocks[ix].Begin();
 
      // block index item
      const block_index_item& index = m_blockIndex[n];
 
      // block starting time
      double jd0 = *block;
 
      // number of components
      int nv = ComponentsForItem( n );
 
      // number of coefficients per subblock
      int nk = index.coefficients;
 
      // subblock time span in days
      int dx = m_blockDelta/index.subblocks;
 
      // subblock index
      ix = ((jd1 - jd0) + jd2)/dx;
      if ( ix == index.subblocks )
         --ix;
 
      // subblock coefficients
      const double* k = block + index.offset + ix*nv*nk;
 
      // compute state variables
      Vector pc( nk );
      pc[0] = 1;
      pc[1] = (2*((jd1 - (jd0 + ix*dx)) + jd2) - dx)/dx;
      double y2 = 2*pc[1];
      for ( int i = 2; i < nk; ++i )
         pc[i] = y2*pc[i-1] - pc[i-2];
      for ( int i = 0; i < nv; ++i )
      {
         r[i] = 0;
         for ( int j = 0; j < nk; ++j )
            r[i] += pc[j] * *k++;
      }
 
      // if requested, compute first derivatives.
      if ( v != nullptr )
      {
         double vf = 2.0*index.subblocks/m_blockDelta;
         Vector vc( nk );
         vc[0] = 0;
         vc[1] = 1;
         vc[2] = 2*y2;
         for ( int i = 3; i < nk; ++i )
            vc[i] = y2*vc[i-1] + 2*pc[i-1] - vc[i-2];
         for ( int i = nv; --i >= 0; )
         {
            v[i] = 0;
            for ( int j = nk; --j >= 0; )
               v[i] += vc[j] * *--k;
            v[i] *= vf;
         }
      }
   }
};
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif  // __PCL_JPLEphemeris_h
 
// ----------------------------------------------------------------------------
// EOF pcl/JPLEphemeris.h - Released 2024-06-18T15:48:54Z