SpiceRotation.h

#ifndef SpiceRotation_h
#define SpiceRotation_h
 
/* SPDX-License-Identifier: CC0-1.0 */
#include <string>
#include <vector>
 
#include <nlohmann/json.hpp>
#include <ale/Orientations.h>
 
#include "Angle.h"
#include "Table.h"
#include "PolynomialUnivariate.h"
#include "Quaternion.h"
 
 
 
#define J2000Code    1
 
namespace Isis {
  class SpiceRotation {
    public:
      // Constructors
      SpiceRotation(int frameCode);
      /*      SpiceRotation( int NaifCode );
      We would like to call refchg instead to avoid the strings.  Currently Naif does
      not have refchg_c, but only the f2c'd refchg.c.*/
      SpiceRotation(int frameCode, int targetCode);
      SpiceRotation(const SpiceRotation &rotToCopy);
 
      // Destructor
      virtual ~SpiceRotation();
 
      // Change the frame (has no effect if cached)
      void SetFrame(int frameCode);
      int Frame();
 
      void SetTimeBias(double timeBias);
      double TimeBias() const;
 
      enum Source {
        Spice,                   
        Nadir,                   
        Memcache,                
        PolyFunction,            
        PolyFunctionOverSpice ,  
        PckPolyFunction          
      };
 
      enum PartialType {
        WRT_RightAscension, 
        WRT_Declination,    
        WRT_Twist           
      };
 
      enum DownsizeStatus {
        Yes,  
        Done, 
        No    
      };
 
      enum FrameType {
        UNKNOWN = 0,      
        INERTL = 1,       
        PCK  = 2,         
        CK = 3,           
        TK = 4,           
        DYN = 5,          
        BPC = 6,          
        NOTJ2000PCK = 7   
      };
 
      void SetEphemerisTime(double et);
      double EphemerisTime() const;
 
      std::vector<double> GetCenterAngles();
 
      std::vector<double> Matrix();
      std::vector<double> AngularVelocity();
 
      // TC
      std::vector<double> ConstantRotation();
      std::vector<double> &ConstantMatrix();
      void SetConstantMatrix(std::vector<double> constantMatrix);
 
      // CJ
      std::vector<double> TimeBasedRotation();
      std::vector<double> &TimeBasedMatrix();
      void SetTimeBasedMatrix(std::vector<double> timeBasedMatrix);
 
      std::vector<double> J2000Vector(const std::vector<double> &rVec);
 
      std::vector<Angle> poleRaCoefs();
 
      std::vector<Angle> poleDecCoefs();
 
      std::vector<Angle> pmCoefs();
 
      std::vector<double> poleRaNutPrecCoefs();
 
      std::vector<double> poleDecNutPrecCoefs();
 
      std::vector<double> pmNutPrecCoefs();
 
      std::vector<Angle> sysNutPrecConstants();
 
      std::vector<Angle> sysNutPrecCoefs();
 
      std::vector<double> ReferenceVector(const std::vector<double> &jVec);
 
      std::vector<double> EvaluatePolyFunction();
 
      void loadPCFromSpice(int CenterBodyCode);
      void loadPCFromTable(const PvlObject &Label);
 
      void MinimizeCache(DownsizeStatus status);
 
      void LoadCache(double startTime, double endTime, int size);
 
      void LoadCache(double time);
 
      void LoadCache(Table &table);
 
      void LoadCache(nlohmann::json &isd);
 
      Table LineCache(const QString &tableName);
 
      void ReloadCache();
 
      Table Cache(const QString &tableName);
      void CacheLabel(Table &table);
 
      void LoadTimeCache();
 
      std::vector<double> Angles(int axis3, int axis2, int axis1);
      void SetAngles(std::vector<double> angles, int axis3, int axis2, int axis1);
 
      bool IsCached() const;
 
      void SetPolynomial(const Source type=PolyFunction);
 
      void SetPolynomial(const std::vector<double> &abcAng1,
                         const std::vector<double> &abcAng2,
                         const std::vector<double> &abcAng3,
                         const Source type = PolyFunction);
 
      void usePckPolynomial();
      void setPckPolynomial(const std::vector<Angle> &raCoeff,
                            const std::vector<Angle> &decCoeff,
                            const std::vector<Angle> &pmCoeff);
 
      void GetPolynomial(std::vector<double> &abcAng1,
                         std::vector<double> &abcAng2,
                         std::vector<double> &abcAng3);
 
      void getPckPolynomial(std::vector<Angle> &raCoeff,
                            std::vector<Angle> &decCoeff,
                            std::vector<Angle> &pmCoeff);
 
      // Set the polynomial degree
      void SetPolynomialDegree(int degree);
      Source GetSource();
      void SetSource(Source source);
      void ComputeBaseTime();
      FrameType getFrameType();
      double GetBaseTime();
      double GetTimeScale();
 
      void SetOverrideBaseTime(double baseTime, double timeScale);
      void SetCacheTime(std::vector<double> cacheTime);
 
      // Derivative methods
      double DPolynomial(const int coeffIndex);
      double DPckPolynomial(PartialType partialVar, const int coeffIndex);
 
      std::vector<double> toJ2000Partial(const std::vector<double> &lookT,
                                         PartialType partialVar, int coeffIndex);
      std::vector<double> ToReferencePartial(std::vector<double> &lookJ,
                                             PartialType partialVar, int coeffIndex);
      void DCJdt(std::vector<double> &dRJ);
 
      double WrapAngle(double compareAngle, double angle);
      void SetAxes(int axis1, int axis2, int axis3);
      std::vector<double> GetFullCacheTime();
      void FrameTrace(double et);
 
      // Return the frame chain for the constant part of the rotation (ends in target)
      std::vector<int>  ConstantFrameChain();
      std::vector<int>  TimeFrameChain();
      void InitConstantRotation(double et);
      bool HasAngularVelocity();
 
      void ComputeAv();
      std::vector<double> Extrapolate(double timeEt);
 
      void checkForBinaryPck();
 
      int cacheSize() {
        if (m_orientation) {
          return m_orientation->getRotations().size();
        }
        return 0;
      }
 
    protected:
      void SetFullCacheParameters(double startTime, double endTime, int cacheSize);
      void setEphemerisTimeMemcache();
      void setEphemerisTimeNadir();
      void setEphemerisTimeSpice();
      void setEphemerisTimePolyFunction();
      void setEphemerisTimePolyFunctionOverSpice();
      void setEphemerisTimePckPolyFunction();
      std::vector<double> p_cacheTime;  
      int p_degree;                     
      int p_axis1;                      
      int p_axis2;                      
      int p_axis3;                      
      ale::Orientations *m_orientation; 
 
    private:
      // method
      void setFrameType();
      std::vector<int> p_constantFrames;  
      std::vector<int> p_timeFrames;      
      double p_timeBias;                  
 
      double p_et;                           
      Quaternion p_quaternion;            
      bool p_matrixSet;                    
      bool m_tOrientationAvailable;  
 
 
      FrameType m_frameType;  
      Source p_source;                    
      int p_axisP;                        
      int p_axisV;                        
      int p_targetCode;                   
 
      double p_baseTime;                  
      double p_timeScale;                 
      bool p_degreeApplied;               
      std::vector<double> p_coefficients[3];  
      bool p_noOverride;                  
      double p_overrideBaseTime;          
      double p_overrideTimeScale;         
      DownsizeStatus p_minimizeCache;     
      double p_fullCacheStartTime;        
      double p_fullCacheEndTime;          
      int p_fullCacheSize;                
      std::vector<double> p_TC;           
      std::vector<double> p_CJ;           
      std::vector<double> p_av;           
      bool p_hasAngularVelocity;          
      std::vector<double> StateTJ();      
      // The remaining items are only used for PCK frame types.  In this case the
      // rotation is  stored as a cache, but the coefficients are available for display
      // or comparison, and the first three coefficient sets can be solved for and
      // updated in jigsaw.   The initial coefficient values are read from a Naif PCK.
      //
      // The general equation for the right ascension of the pole is
      //
      // raPole  =  raPole[0] + raPole[1]*Time  + raPole[2]*Time**2 + raNutPrec,
      //    where
      //    raNutPrec  =  raNutPrec1[0]*sin(sysNutPrec[0][0] + sysNutPrec[0][1]*Time) +
      //                  raNutPrec1[1]*sin(sysNutPrec[1][0] + sysNutPrec[1][1]*Time) + ...
      //                  raNutPrec1[N-1]*sin(sysNutPrec[N-1][0] + sysNutPrec[N-1][1]*Time) +
      // (optional for multiples of nutation precession angles)
      //                  raNutPrec2[0]*sin(2*(sysNutPrec[0][0] + sysNutPrec[0][1]*Time)) +
      //                  raNutPrec2[1]*sin(2*(sysNutPrec[1][0] + sysNutPrec[1][1]*Time)) + ...
      //                  raNutPrec2[N-1]*sin(2*(sysNutPrec[N-1][0] + sysNutPrec[N-1][1]*Time)) +
      //                  raNutPrecM[0]*sin(M*(sysNutPrec[0][0] + sysNutPrec[0][1]*Time)) +
      //                  raNutPrecM[1]*sin(M*(sysNutPrec[1][0] + sysNutPrec[1][1]*Time)) + ...
      //                  raNutPrecM[N-1]*sin(M*(sysNutPrec[N-1][0] + sysNutPrec[N-1][1]*Time)) +
      //
      // The general equation for the declination of the pole is
      //
      // decPole  =  p_decPole[0] + p_decPole[1]*Time  + p_decPole[2]*Time**2 + decNutPrec,
      //    where
      //    decNutPrec  =  decNutPrec1[0]*cos(sysNutPrec[0][0] + sysNutPrec[0][1]*Time) +
      //                   decNutPrec1[1]*cos(sysNutPrec[1][0] + sysNutPrec[1][1]*Time) + ...
      //                   decNutPrec1[N-1]*cos(sysNutPrec[N-1][0] + sysNutPrec[N-1][1]*Time) +
      //                   decNutPrec2[0]*cos(2*(sysNutPrec[0][0] + sysNutPrec[0][1]*Time)) +
      //                   decNutPrec2[1]*cos(2*(sysNutPrec[1][0] + sysNutPrec[1][1]*Time)) + ...
      //                   decNutPrec2[N-1]*cos(2*(sysNutPrec[N-1][0] + sysNutPrec[N-1][1]*Time)) +
      // (optional for multiples of nutation precession angles)
      //                   decNutPrecM[0]*sin(M*(sysNutPrec[0][0] + sysNutPrec[0][1]*Time)) +
      //                   decNutPrecM[1]*sin(M*(sysNutPrec[1][0] + sysNutPrec[1][1]*Time)) + ...
      //                   decNutPrecM[N-1]*sin(M*(sysNutPrec[N-1][0] + sysNutPrec[N-1][1]*Time))
      //
      //     and Time is julian centuries since J2000.
      //
      // The general equation for the prime meridian rotation is
      //
      // pm  =  p_pm[0] + p_pm[1]*Dtime  + p_pm[2]*Dtime**2 + pmNutPrec,
      //    where
      //    pmNutPrec  =  pmNutPrec1[0]*sin(sysNutPrec[0][0] + sysNutPrec[0][1]*Time) +
      //                  pmNutPrec1[1]*sin(sysNutPrec[1][0] + sysNutPrec[1][1]*Time) + ...
      //                  pmNutPrec1[N-1]*sin(sysNutPrec[N-1][0] + sysNutPrec[N-1][1]*Time) +
      // (optional for multiples of nutation precession angles)
      //                  pmNutPrec2[0]*sin(2*(sysNutPrec[0][0] + sysNutPrec[0][1]*Time)) +
      //                  pmNutPrec2[1]*sin(2*(sysNutPrec[1][0] + sysNutPrec[1][1]*Time)) + ...
      //                  pmNutPrec2[N-1]*sin(2*(sysNutPrec[N-1][0] + sysNutPrec[N-1][1]*Time)) +
      //                  pmNutPrecM[0]*sin(M*(sysNutPrec[0][0] + sysNutPrec[0][1]*Time)) +
      //                  pmNutPrecM[1]*sin(M*(sysNutPrec[1][0] + sysNutPrec[1][1]*Time)) + ...
      //                  pmNutPrecM[N-1]*sin(M*(sysNutPrec[N-1][0] + sysNutPrec[N-1][1]*Time))
      //
      //     Time is interval in Julian centuries since the standard epoch,
      //     dTime is interval in days from the standard epoch (J2000),
      //
      //     N is the number of nutation/precession terms for the planetary system of the target
      //     body,  (possibly including multiple angles as unique terms,
      //             ie. 2*sysNutPrec[0][0] + sysNutPrec[][1]*Time).
      //
      //     Many of the constants in this equation are 0. for a given body.
      //
      //     M is included as an option for future improvements.  M = highest multiple (period)
      //     of any of the nutation/precession angles included in the equations.
      //
      //     ***NOTE*** Currently Naif stores multiples (amplitudes) as if they were additional
      //                nutation/precession terms (periods) in the equation.  This method works as
      //                long as jigsaw does not solve for those values.  In order to solve for
      //                those values, the multiples will need to be known so that the partial
      //                derivatives can be correctly calculated.  Some possible ways of doing this
      //                are 1) Convince Naif to change their data format indicating the relation
      //                      2) Make an Isis version of the PCK data and have Isis software to
      //                          calculate the rotation and partials.
      //                      3) Have an Isis addendum file that identifies the repeated periods
      //                          and software to apply them when calculating the rotation and partials.
      //
      //                For now this software will handle any terms with the same period and different
      //                amplitudes as unique terms in the equation (raNutPrec, decNutPrec,
      //                and pmNutPrec).
      //
      // The next three vectors will have length 3 (for a quadratic polynomial) if used.
      std::vector<Angle>m_raPole;       
      std::vector<Angle>m_decPole;      
      std::vector<Angle>m_pm ;          
      //
      // Currently multiples (terms with periods matching other terms but varying amplitudes)
      // are handled as additional terms added to the end of the vector as Naif does (see
      // comments in any of the standard Naif PCK.
      std::vector<double>m_raNutPrec;    
      std::vector<double>m_decNutPrec;  
      std::vector<double>m_pmNutPrec;   
 
      // The periods of bodies in the same system are modeled with a linear equation
      std::vector<Angle>m_sysNutPrec0; 
      std::vector<Angle>m_sysNutPrec1; 
 
      // The following scalars are used in the IAU equations to convert p_et to the appropriate time
      // units for calculating target body ra, dec, and w.  These need to be initialized in every
      // constructor.
      static const double m_centScale;
      static const double m_dayScale;
 
  };
};
 
#endif