GravityCalculator.hpp

 
 
#pragma once
 
#include "AstGlobal.h"
#include "GravityField.hpp"
#include "AstUtil/StringView.hpp"
// #include "AstMath/Matrix.hpp"
// #include "AstMath/Vector.hpp"
#include <cmath>
#include <vector>
#include <memory>
 
AST_NAMESPACE_BEGIN
 
AST_CORE_API void aPointMassField(const Vector3d& position, double GM, Vector3d& acceleration, Matrix3d& hessian);
 
class GravityCalculator1;
class GravityCalculator2;
class GravityCalculator3;
 
#ifdef _AST_USE_GRAVITY_CALCULATOR_1
    #define _AST_ENABLE_GRAVITY_CALCULATOR_1
    using GravityCalculatorDefault = GravityCalculator1;
#else
    using GravityCalculatorDefault = GravityCalculator3;
#endif
 
class AST_CORE_API GravityCalculator
{
public:
    GravityCalculator();
    
    GravityCalculator(const GravityField &gravityField, int degree=-1, int order=-1);
 
    GravityCalculator(GravityField &&gravityField, int degree=-1, int order=-1);
 
 
    const GravityField& getGravityField() const { return gravityField_; }
 
    void updateVariations(const TimePoint& tp){gravityField_.updateVariations(tp);}
 
    // 在这里不能公开提供返回可修改引用接口，否则会导致重力场系数等参数被修改
    // 
    // 例如: 在初始化重力场计算对象后，
    // 通过重力场的可修改引用调用load函数重新加载重力场模型
    // 可能导致加载重力场阶次比计算所需的阶次小，从而导致执行崩溃问题
    // 或者导致加载重力场阶次比计算所需的阶次大，但是计算时并没有使用高阶项，容易导致与预期不符
    // 
    // 因此，在这里不能公开提供可修改引用，只能提供const引用
protected:
    GravityField& getGravityField() { return gravityField_; }
public:
    int getDegree() const { return degree_; }
 
    int getOrder() const { return order_; }
 
    void setDegreeForPartial(unsigned int degree);
 
    void setOrderForPartial(unsigned int order);
 
    int getDegreeForPartial() const { return degreeForPartial_; }
 
    int getOrderForPartial() const { return orderForPartial_; }
    
    virtual ~GravityCalculator () = default;
 
    virtual void calcPertAcceleration (const Vector3d &positionCBF, Vector3d &accelerationCBF) = 0;
 
    void calcTotalAcceleration (const Vector3d &positionCBF, Vector3d &accelerationCBF);
private:
    void initDegreeOrder(int degree, int order);
protected:
    GravityField gravityField_;                 
    unsigned int degree_{0};                    
    unsigned int order_{0};                     
    unsigned int degreeForPartial_{0};          
    unsigned int orderForPartial_{0};           
};
 
#ifdef _AST_ENABLE_GRAVITY_CALCULATOR_1
 
class AST_CORE_API GravityCalculator1 final: public GravityCalculator
{
public:
    GravityCalculator1();
    GravityCalculator1(const GravityField &gravityField, int degree, int order);
    ~GravityCalculator1();
 
    void calcPertAcceleration (const Vector3d &positionCBF, Vector3d &accelerationCBF) final;
 
protected:
    class Impl;
    std::unique_ptr<Impl> impl_;
};
 
#endif
 
 
class AST_CORE_API GravityCalculator2 final: public GravityCalculator
{
public:
    GravityCalculator2(const GravityField &gravityField, int degree, int order);
    
    void calcPertAcceleration (const Vector3d &positionCBF, Vector3d &accelerationCBF) final;
 
private:
    void initializeRecursionCoefficients();
    
    void computeGradient(const Vector3d &position, Vector3d &gradient);
    
    int computeTesseral(int m, int degree, int index, double t, double u, double tOu);
 
private:
    std::vector<double> gnmOj_{};
 
    std::vector<double> hnmOj_{};
 
    std::vector<double> enm_{};
 
    std::vector<double> sectorial_{};
private:
    std::vector<double> pnm0Plus2{};
    std::vector<double> pnm0Plus1{};
    std::vector<double> pnm0{};
    std::vector<double> pnm1{};
private:
    std::vector<double> aOrN{};
    std::vector<double> cosLambdaN{};
    std::vector<double> sinLambdaN{};
};
 
 
class AST_CORE_API GravityCalculator3 final: public GravityCalculator
{
public:
    GravityCalculator3();
    GravityCalculator3(const GravityField &gravityField, int degree, int order);
    GravityCalculator3(GravityField &&gravityField, int degree, int order);
 
    GravityCalculator3(const GravityCalculator3&) = delete;
    GravityCalculator3& operator=(const GravityCalculator3&) = delete;
 
    ~GravityCalculator3() override;
    
    void calcPertAcceleration(const Vector3d &positionCBF, Vector3d &accelerationCBF) final;
    void calcPertAcceleration(const Vector3d &positionCBF, Vector3d &accelerationCBF, Matrix3d &gradient);
    template<bool fillgradient>
    void calcPertAcceleration(const Vector3d &positionCBF, Vector3d &accelerationCBF, Matrix3d *gradient);
    using GravityCalculator::calcTotalAcceleration;
    void calcTotalAcceleration (const Vector3d &positionCBF, Vector3d &accelerationCBF, Matrix3d &gradient);
protected:
    void init();
    void deinit();
private:
    double        Factor;  // Factor = 1 (magnetic) or -mu (gravity)
    double**      C{};       // Normalized harmonic coefficients
    double**      S{};       // Normalized harmonic coefficients
    double**      A{};       // Normalized 'derived' Assoc. Legendre Poly
    double**      V{};       // Normalization factor
    double*       Re{};      // powers of projection of pos onto x_ecf (re)
    double*       Im{};      // powers of projection of pos onto y_ecf (im)
    double**      N1{};      // Temporary
    double**      N2{};      // Temporary
    double**      VR01{};    // Temporary
    double**      VR11{};    // Temporary
    double**      VR02{};    // Temporary
    double**      VR12{};    // Temporary
    double**      VR22{};    // Temporary
};
 
 
class AST_CORE_API GravityCalculator4 final: public GravityCalculator
{
public:
    using GravityCalculator::GravityCalculator;
 
    void calcPertAcceleration (const Vector3d &positionCBF, Vector3d &accelerationCBF) final;
};
 
 
class AST_CORE_API GravityCalculator5 final: public GravityCalculator
{
public:
    using GravityCalculator::GravityCalculator;
 
    void calcPertAcceleration (const Vector3d &positionCBF, Vector3d &accelerationCBF) final;
};
 
 
AST_NAMESPACE_END
 
AST_DECL_TYPE_ALIAS(GravityCalculator)
AST_DECL_TYPE_ALIAS(GravityCalculator1)
AST_DECL_TYPE_ALIAS(GravityCalculator2)
AST_DECL_TYPE_ALIAS(GravityCalculator3)
AST_DECL_TYPE_ALIAS(GravityCalculator4)
AST_DECL_TYPE_ALIAS(GravityCalculator5)