Position.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.7.0
// ----------------------------------------------------------------------------
// pcl/Position.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_Position_h
#define __PCL_Position_h
 
 
#include <pcl/Defs.h>
 
#include <pcl/AutoPointer.h>
#include <pcl/EphemerisFile.h>
#include <pcl/Matrix.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
struct PCL_CLASS StarPosition
{
   double    alpha = 0;    
   double    delta = 0;    
   double    muAlpha = 0;  
   double    muDelta = 0;  
   double    p = 0;        
   double    v = 0;        
   TimePoint t0 = TimePoint::J2000(); 
 
   StarPosition() = default;
 
   StarPosition( const StarPosition& x )
      : alpha( x.alpha )
      , delta( x.delta )
      , muAlpha( x.muAlpha )
      , muDelta( x.muDelta )
      , p( x.p )
      , v( x.v )
      , t0( x.t0 )
   {
   }
 
   StarPosition& operator =( const StarPosition& x )
   {
      alpha = x.alpha;
      delta = x.delta;
      muAlpha = x.muAlpha;
      muDelta = x.muDelta;
      p = x.p;
      v = x.v;
      t0 = x.t0;
      return *this;
   }
 
   StarPosition( double ra, double dec,
                 double properMotionRA = 0, double properMotionDec = 0,
                 double parallax = 0,
                 double radialVelocity = 0,
                 TimePoint epoch = TimePoint::J2000() )
      : alpha( ((ra = Mod( ra, 360.0 )) < 0) ? ra + 360 : ra )
      , delta( Range( dec, -90.0, +90.0 ) )
      , muAlpha( properMotionRA )
      , muDelta( properMotionDec )
      , p( parallax )
      , v( radialVelocity )
      , t0( epoch )
   {
   }
 
private:
 
   uint64 m_uniqueId = UniqueId();
 
   static uint64 UniqueId();
 
   friend class PCL_CLASS Position;
};
 
// ----------------------------------------------------------------------------
 
struct PCL_CLASS ObserverPosition
{
   double lambda = 0;             
   double phi = 0;                
   double h = 0;                  
   double a = ea_eq_radius_IAU2009;   
   double f = ea_flattening_IERS2010; 
   Vector r0 = Vector( 0, 0, 0 ); 
   bool   cioBased = false;       
 
   constexpr static double ea_eq_radius_IAU2009 = 6378136.6;
 
   constexpr static double ea_flattening_IERS2010 = 1/298.25642;
 
   ObserverPosition() = default;
 
   ObserverPosition( const ObserverPosition& ) = default;
 
   ObserverPosition( ObserverPosition&& ) = default;
 
   ObserverPosition& operator =( const ObserverPosition& ) = default;
 
   ObserverPosition& operator =( ObserverPosition&& ) = default;
 
   ObserverPosition( double longitude, double latitude,
                     double height = 0,
                     double equatorialRadius = ea_eq_radius_IAU2009,
                     double flattening = ea_flattening_IERS2010,
                     const Vector& regionalCenter = Vector( 0, 0, 0 ),
                     bool useCIO = false )
      : lambda( ((longitude = Mod( longitude, 360.0 )) < 0) ? longitude + 360 : longitude )
      , phi( Range( latitude, -90.0, +90.0 ) )
      , h( Max( 0.0, height ) )
      , a( Max( 0.0, equatorialRadius ) )
      , f( Max( 0.0, flattening ) )
      , r0( regionalCenter )
      , cioBased( useCIO )
   {
   }
};
 
// ----------------------------------------------------------------------------
 
class PCL_CLASS Position
{
public:
 
   Position( TimePoint t, const IsoString& timescale = "TT" );
 
   Position( const Position& ) = default;
 
   Position( Position&& ) = default;
 
   Position& operator =( const Position& ) = default;
 
   Position& operator =( Position&& ) = default;
 
   Vector True( EphemerisFile::Handle& H )
   {
      (void)Geometric( H );
      return m_U0;
   }
 
   Vector True( const StarPosition& S )
   {
      (void)Geometric( S );
      return m_U0;
   }
 
   double TrueDistance( EphemerisFile::Handle& H )
   {
      return True( H ).L2Norm();
   }
 
   double TrueDistance( const StarPosition& S )
   {
      return True( S ).L2Norm();
   }
 
   double LightTravelTime( EphemerisFile::Handle& H )
   {
      (void)Geometric( H );
      return m_tau;
   }
 
   Vector Geometric( EphemerisFile::Handle& H );
 
   Vector Geometric( const StarPosition& S );
 
   Vector Astrometric( EphemerisFile::Handle& H );
 
   Vector Astrometric( const StarPosition& S );
 
   Vector Proper( EphemerisFile::Handle& H );
 
   Vector Proper( const StarPosition& S );
 
   Vector Apparent( EphemerisFile::Handle& H );
 
   Vector Apparent( const StarPosition& S );
 
   Vector Intermediate( EphemerisFile::Handle& H );
 
   Vector Intermediate( const StarPosition& S );
 
   void SetObserver( const ObserverPosition& observer );
 
   void SetGeocentric();
 
   ObserverPosition Observer() const
   {
      if ( m_observer )
         return *m_observer;
      return ObserverPosition();
   }
 
   bool IsTopocentric() const
   {
      return m_observer;
   }
 
   bool IsPolarMotionEnabled() const
   {
      return m_usePolarMotion;
   }
 
   void EnablePolarMotion( bool enable = true )
   {
      m_usePolarMotion = enable;
   }
 
   void DisablePolarMotion( bool disable = true )
   {
      EnablePolarMotion( !disable );
   }
 
   TimePoint TDB() const
   {
      return m_t;
   }
 
   TimePoint Teph() const
   {
      return m_t;
   }
 
   TimePoint TT() const
   {
      return m_tt;
   }
 
   TimePoint UTC() const
   {
      return m_utc;
   }
 
   TimePoint UT1() const
   {
      return m_ut1;
   }
 
   Vector BarycentricPositionOfEarth() const
   {
      return m_Eb;
   }
 
   Vector BarycentricVelocityOfEarth() const
   {
      return m_Edb;
   }
 
   Vector BarycentricPositionOfSun() const
   {
      return m_Sb;
   }
 
   Vector HeliocentricPositionOfEarth() const
   {
      return m_Eh;
   }
 
   void InitEquinoxBasedParameters();
 
   void InitCIOBasedParameters();
 
   Matrix EquinoxBiasPrecessionNutationMatrix()
   {
      InitEquinoxBasedParameters();
      return m_M;
   }
 
   Matrix EquinoxBiasPrecessionNutationInverseMatrix()
   {
      InitEquinoxBasedParameters();
      return m_Minv;
   }
 
   Matrix CIOBiasPrecessionNutationMatrix()
   {
      InitCIOBasedParameters();
      return m_C;
   }
 
   Matrix CIOBiasPrecessionNutationInverseMatrix()
   {
      InitCIOBasedParameters();
      return m_Cinv;
   }
 
   Vector CIP()
   {
      InitEquinoxBasedParameters();
      return Vector( m_X, m_Y );
   }
 
   Vector CIP_ITRS() const;
 
   double CIO()
   {
      InitEquinoxBasedParameters();
      return m_s;
   }
 
   double EO()
   {
      InitEquinoxBasedParameters();
      return m_EO;
   }
 
   double ERA()
   {
      InitEquinoxBasedParameters();
      return m_ERA;
   }
 
   double GAST()
   {
      InitEquinoxBasedParameters();
      return m_GAST;
   }
 
   double EpsA() const
   {
      // Mean obliquity of the ecliptic, IAU 2006 precession model.
      if ( m_M.IsEmpty() )
         return AsRad( Poly( m_TT, { 84381.406, -46.836769, -0.0001831, 0.00200340, -0.000000576, -0.0000000434 } ) );
      return m_epsa;
   }
 
   DPoint NutationAngles()
   {
      InitEquinoxBasedParameters();
      return DPoint( m_dpsi, m_deps );
   }
 
   double PhaseAngle( EphemerisFile::Handle& H );
 
   double PhaseAngle( const StarPosition& S );
 
   bool CanComputeApparentVisualMagnitude( const EphemerisFile::Handle& H ) const;
 
   Optional<double> ApparentVisualMagnitude( EphemerisFile::Handle& H );
 
   bool CanComputeCometApparentVisualTotalMagnitude( const EphemerisFile::Handle& H ) const;
 
   Optional<double> CometApparentVisualTotalMagnitude( EphemerisFile::Handle& H );
 
   bool CanComputeCometApparentVisualNuclearMagnitude( const EphemerisFile::Handle& H ) const;
 
   Optional<double> CometApparentVisualNuclearMagnitude( EphemerisFile::Handle& H );
 
   DPoint EquatorialToHorizontal( const DPoint& q )
   {
      return EquatorialToHorizontal( q.x, q.y );
   }
 
   DPoint EquatorialToHorizontal( double ra, double dec )
   {
      if ( !IsTopocentric() )
         return DPoint( 0 );
 
      if ( m_observer->cioBased )
         InitCIOBasedParameters();
      else
         InitEquinoxBasedParameters();
 
      double sh, ch; SinCos( Norm2Pi( (m_observer->cioBased ? m_ERA : m_GAST) - ra + Rad( m_observer->lambda ) ), sh, ch );
      double sd, cd; SinCos( dec, sd, cd );
      return DPoint( Norm2Pi( ArcTan( -cd*sh, sd*m_cphi - cd*ch*m_sphi ) ),
                     ArcSin( sd*m_sphi + cd*ch*m_cphi ) );
   }
 
   static Vector EquatorialToEcliptic( const Vector& q, double se, double ce )
   {
      // Rx(eps)*q
      return Vector( q[0], q[1]*ce + q[2]*se, q[2]*ce - q[1]*se );
   }
 
   static Vector EquatorialToEcliptic( const Vector& q, double eps )
   {
      double se, ce; SinCos( eps, se, ce );
      return EquatorialToEcliptic( q, se, ce );
   }
 
   static DPoint EquatorialToEcliptic( const DPoint& q, double se, double ce )
   {
      DPoint e;
      EquatorialToEcliptic( Vector::FromSpherical( q.x, q.y ), se, ce ).ToSpherical2Pi( e.x, e.y );
      return e;
   }
 
   static DPoint EquatorialToEcliptic( const DPoint& q, double eps )
   {
      double se, ce; SinCos( eps, se, ce );
      return EquatorialToEcliptic( q, se, ce );
   }
 
   static Vector EclipticToEquatorial( const Vector& q, double se, double ce )
   {
      // Rx(eps)*q
      return Vector( q[0], q[1]*ce - q[2]*se, q[2]*ce + q[1]*se );
   }
 
   static Vector EclipticToEquatorial( const Vector& q, double eps )
   {
      double se, ce; SinCos( eps, se, ce );
      return EclipticToEquatorial( q, se, ce );
   }
 
   static DPoint EclipticToEquatorial( const DPoint& q, double se, double ce )
   {
      DPoint e;
      EclipticToEquatorial( Vector::FromSpherical( q.x, q.y ), se, ce ).ToSpherical2Pi( e.x, e.y );
      return e;
   }
 
   static DPoint EclipticToEquatorial( const DPoint& q, double eps )
   {
      double se, ce; SinCos( eps, se, ce );
      return EclipticToEquatorial( q, se, ce );
   }
 
   static Vector ICRSEquatorialToGalactic( const Vector& q )
   {
      return Matrix( +0.494055821648, -0.054657353964, -0.445679169947,
                     -0.872844082054, -0.484928636070, +0.746511167077,
                     -0.867710446378, -0.198779490637, +0.455593344276 )*q;
   }
 
   static DPoint ICRSEquatorialToGalactic( const DPoint& q )
   {
      DPoint g;
      ICRSEquatorialToGalactic( Vector::FromSpherical( q.x, q.y ) ).ToSpherical2Pi( g.x, g.y );
      return g;
   }
 
   static Vector GalacticToICRSEquatorial( const Vector& g )
   {
      return Matrix( -0.823971452454, +1.289170419979, -2.918407488771,
                     -2.840810442223, -1.835975400200, +0.229340705201,
                     -2.808784424949, +1.654265575150, -3.263314653103 )*g;
   }
 
   static DPoint GalacticToICRSEquatorial( const DPoint& g )
   {
      DPoint q;
      GalacticToICRSEquatorial( Vector::FromSpherical( g.x, g.y ) ).ToSpherical2Pi( q.x, q.y );
      return q;
   }
 
private:
 
   // TDB
   TimePoint   m_t;
   // TT
   TimePoint   m_tt;
   // UTC
   TimePoint   m_utc;
   // UT1
   TimePoint   m_ut1;
   // TT in Julian centuries since J2000.0.
   double      m_TT;
   // Barycentric ICRS position and velocity of the Earth.
   Vector      m_Eb, m_Edb;
   // Geocentric ICRS position and velocity of the observer, in km and km/day.
   Vector      m_G, m_Gd;
   // Barycentric ICRS position of the Sun.
   Vector      m_Sb;
   // Heliocentric position of the Earth.
   Vector      m_Eh;
   // Distance Earth-Sun.
   double      m_E;
   // Heliocentric position of the observer (au).
   Vector      m_Oh;
   // Distance observer-Sun (au).
   double      m_O;
   // Bias+precession angles.
   double      m_gamb, m_phib, m_psib, m_epsa;
   // Nutation angles.
   double      m_dpsi, m_deps;
   // Position of the Celestial Intermediate Pole (CIP) in GCRS.
   double      m_X, m_Y;
   // CIO locator.
   double      m_s;
   // Combined bias-precession-nutation matrix (NPB).
   Matrix      m_M, m_Minv;
   // Combined bias-precession-nutation matrix (NPB_CIO).
   Matrix      m_C, m_Cinv;
   // Equation of the origins.
   double      m_EO;
   // Earth rotation angle.
   double      m_ERA;
   // Greenwich apparent sidereal time.
   double      m_GAST;
   // Light-travel time in days.
   double      m_tau = 0;
   // Barycentric position.
   Vector      m_ub;
   // True geocentric position.
   Vector      m_U0;
   // Geocentric position.
   Vector      m_U;
   // Astrometric place.
   Vector      m_u1;
   // Proper place.
   Vector      m_u2;
   // Apparent place.
   Vector      m_u3e;
   // Intermediate place.
   Vector      m_u3i;
 
   // Handles for calculation of fundamental ephemerides and nutation angles.
   AutoPointerCloner<EphemerisFile::Handle> m_TT_TDB, m_HE, m_HS, m_HN;
 
   // Current observer for calculation of topocentric coordinates.
   AutoPointerCloner<ObserverPosition> m_observer;
   double                              m_sphi, m_cphi;
 
   // Special case flags.
   bool m_isMoon = false, m_isSun = false, m_isStar = false;
 
   // Whether to account for polar motion (CIP->ITRS) in calculation of
   // topocentric positions.
   bool m_usePolarMotion = true;
 
   // Unique identifier of the object whose positions are being calculated.
   uint64 m_uniqueObjectId = 0;
 
   template <class T>
   bool Validate( const T& obj )
   {
      if ( obj.m_uniqueId != m_uniqueObjectId )
      {
         m_U0 = m_U = m_ub = m_u1 = m_u2 = m_u3e = m_u3i = Vector();
         m_tau = 0;
         m_isMoon = m_isSun = m_isStar = false;
         m_uniqueObjectId = obj.m_uniqueId;
         return false;
      }
      return true;
   }
 
   Vector Deflection();
   Vector Aberration();
 
   static double CIOLocator( double T, double X, double Y );
 
   /*
    * Astronomical constants, IAU 2009/2012 and IERS 2003/2010.
    */
   constexpr static double au_km    = 149597870.7;          // astronomical unit (km)
   constexpr static double c_km_s   = 299792.458;           // speed of light (km/s)
   constexpr static double c_km_day = c_km_s*86400;         // speed of light (km/day)
   constexpr static double c_au_day = (c_km_s/au_km)*86400; // speed of light (au/day)
   constexpr static double earth_omega = 7.292115e-5;       // angular velocity of Earth in radians/s
};
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif  // __PCL_Position_h
 
// ----------------------------------------------------------------------------
// EOF pcl/Position.h - Released 2024-06-18T15:48:54Z