SpiceRotation.h

Go to the documentation of this file.

#ifndef SpiceRotation_h
#define SpiceRotation_h

#include <string>
#include <vector>

//#include <SpiceUsr.h>
//#include <SpiceZfc.h>
//#include <SpiceZmc.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);

      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);

      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);

      // 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();

    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;  
      std::vector<std::vector<double> > p_cache; 
      int p_degree;                     
      int p_axis1;                      
      int p_axis2;                      
      int p_axis3;                      

    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<std::vector<double> > p_cacheAv;
      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