ForceModel.hpp

Go to the documentation of this file.

#pragma ident "$Id: ForceModel.hpp 2535 2011-03-25 15:58:06Z ccutlip $"

#ifndef GPSTK_FORCE_MODEL_HPP
#define GPSTK_FORCE_MODEL_HPP

//============================================================================
//
//  This file is part of GPSTk, the GPS Toolkit.
//
//  The GPSTk is free software; you can redistribute it and/or modify
//  it under the terms of the GNU Lesser General Public License as published
//  by the Free Software Foundation; either version 2.1 of the License, or
//  any later version.
//
//  The GPSTk is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with GPSTk; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
//  Wei Yan - Chinese Academy of Sciences . 2009, 2010
//
//============================================================================


#include "Vector.hpp"
#include "Matrix.hpp"
#include "Spacecraft.hpp"
#include "EarthBody.hpp"

namespace gpstk
{

   class ForceModel
   {
   public:

      enum ForceModelIndex
      {
         FMI_BASE       = 1000,  

         FMI_GEOEARTH,           
         FMI_GEOSUN,             
         FMI_GEOMOON,            
         FMI_DRAG,               
         FMI_SRP,                
         FMI_RELATIVE,           
         FMT_EMPIRICAL,          
         
         //... add more here

         FMI_LIST   = 2000      
      };

      enum ForceModelType
      {
         Cd,               // Coefficient of drag
         Cr               // Coefficient of Reflectivity
      };

      ForceModel()
      {         
         a.resize(3,0.0);
         da_dr.resize(3,3,0.0);
         da_dv.resize(3,3,0.0);
         da_dp.resize(3,0,0.0);      // default np = 0;

         da_dcd.resize(3,1,0.0);
         da_dcr.resize(3,1,0.0);
      }

      virtual ~ForceModel(){}


      virtual void doCompute(UTCTime t, EarthBody& bRef, Spacecraft& sc)
      {
         a.resize(3,0.0);
         da_dr.resize(3,3,0.0);
         da_dv.resize(3,3,0.0);
         da_dp.resize(3,0,0.0);      // default np = 0;

         da_dcd.resize(3,1,0.0);
         da_dcr.resize(3,1,0.0);

      }
         
      virtual std::string modelName() const
      { return "ForceModel"; }


      virtual int forceIndex() const
      { return FMI_BASE; }


      virtual Vector<double> getAccel() const
      { return a; }

      virtual Matrix<double> partialR() const
      { return da_dr; }

      virtual Matrix<double> partialV() const
      { return da_dv; }

      virtual Matrix<double> partialP() const
      { return da_dp; }

      virtual Matrix<double> partialCd() const
      { return da_dcd; } 

      virtual Matrix<double> partialCr() const
      { return da_dcr; }

      int getNP() const
      { return da_dp.cols(); }

      Matrix<double> getAMatrix() const
      {
            /* A Matrix
            |                        |
            | 0         I      0      |
            |                        |
         A =| da_dr      da_dv   da_dp  |
            |                        |
            | 0         0      0      |
            |                        |
            */

         const int np = da_dp.cols();

         gpstk::Matrix<double> A(6+np,6+np,0.0);

         A(0,3) = 1.0;
         A(1,4) = 1.0;
         A(2,5) = 1.0;

         // da_dr
         A(3,0) = da_dr(0,0);
         A(3,1) = da_dr(0,1);
         A(3,2) = da_dr(0,2);
         A(4,0) = da_dr(1,0);
         A(4,1) = da_dr(1,1);
         A(4,2) = da_dr(1,2);
         A(5,0) = da_dr(2,0);
         A(5,1) = da_dr(2,1);
         A(5,2) = da_dr(2,2);

         // da_dv
         A(3,3) = da_dv(0,0);
         A(3,4) = da_dv(0,1);
         A(3,5) = da_dv(0,2);
         A(4,3) = da_dv(1,0);
         A(4,4) = da_dv(1,1);
         A(4,5) = da_dv(1,2);
         A(5,3) = da_dv(2,0);
         A(5,4) = da_dv(2,1);
         A(5,5) = da_dv(2,2);

         // da_dp
         for(int i=0;i<np;i++)
         {
            A(3,6+i) = da_dp(0,i);
            A(4,6+i) = da_dp(1,i);
            A(5,6+i) = da_dp(2,i);
         }

         return A;

      }  // End of method 'getAMatrix()'

      void test()
      {
         /*
         cout<<"test Force Model"<<endl;

         a.resize(3,2.0);
         da_dr.resize(3,3,3.0);
         da_dv.resize(3,3,4.0);
         da_dp.resize(3,2,5.0);
         writeToFile("default.fm");

         // it work well
         */
      }

   protected:

      Vector<double> a;         // 3
      
      Matrix<double> da_dr;      // 3*3
      
      Matrix<double> da_dv;      // 3*3
      
      Matrix<double> da_dp;      // 3*np
         
      Matrix<double> da_dcd;      // 3*1
         
      Matrix<double> da_dcr;      // 3*1

      

   }; // End of 'class ForceModel'

   inline std::ostream& operator<<( std::ostream& s,
                                    const gpstk::ForceModel& fm )
   {
      Vector<double> a = fm.getAccel();
      Matrix<double> da_dr = fm.partialR();
      Matrix<double> da_dv = fm.partialV();
      Matrix<double> da_dp = fm.partialP();

      s<<"a ["<<a.size()<<"]\n{\n"
         <<a<<endl<<"}\n\n";

      s<<"da/dr ["<<da_dr.rows()<<","<<da_dr.cols()<<"]\n{\n"
         <<da_dr<<endl<<"}\n\n";

      s<<"da/dv ["<<da_dv.rows()<<","<<da_dv.cols()<<"]\n{\n"
         <<da_dv<<endl<<"}\n\n";

      s<<"da/dp ["<<da_dp.rows()<<","<<da_dp.cols()<<"]\n{\n"
         <<da_dp<<endl<<"}\n\n";

      Matrix<double> A = fm.getAMatrix();

      s<<"A = ["<<A.rows()<<","<<A.cols()<<"]\n{\n"
         <<A<<endl<<"}\n\n";

      return s;

   }

      // @}

}  // End of namespace 'gpstk'


#endif  // GPSTK_FORCE_MODEL_HPP

---