RungeKuttaIntegrator.cpp

Go to the documentation of this file.

#include "RungeKuttaIntegrator.h"
namespace O_SESSAME {

RungeKuttaIntegrator::RungeKuttaIntegrator()
{
    m_NumSteps = 1;
}

/*
Matrix RungeKuttaIntegrator::Integrate(const Vector &_propTime, odeFunc _FuncPtr, const Vector &_initialConditions, const Matrix &_constants, vectorFuncPtr _vectorFuncPtr)
{
        double timeInitial = _propTime(VectorIndexBase);
        double t = timeInitial;
        double h = (_propTime(VectorIndexBase + 1) - t) / m_NumSteps;

        int numEqs = _initialConditions.getIndexCount();
        Vector inputs = _initialConditions;
        Vector K1(numEqs);
        Vector K2(numEqs);
        Vector K3(numEqs);
        Vector K4(numEqs);

        Matrix output(m_NumSteps + 1, numEqs + 1);

        output(MatrixIndexBase,MatrixIndexBase) = timeInitial;
        output(MatrixIndexBase,_(MatrixIndexBase+1,MatrixIndexBase+numEqs)) = ~_initialConditions;

        Vector temp(numEqs);
        for (int ii = 1; ii <= m_NumSteps; ++ii)
        {
            K1 = h * _FuncPtr(t, inputs, _constants, _vectorFuncPtr);
            temp = inputs + K1 / 2.0;
            K2 = h * _FuncPtr(t + h/2, temp, _constants, _vectorFuncPtr);
            temp = inputs + K2 / 2.0;
            K3 = h * _FuncPtr(t + h/2, temp, _constants, _vectorFuncPtr);
            temp = inputs + K3;
            K4 = h * _FuncPtr(t + h, temp, _constants, _vectorFuncPtr);
            for (int jj = MatrixIndexBase; jj < MatrixIndexBase + numEqs; ++jj)
            {
                inputs(jj) += (K1(jj)
                                 + 2.0 * K2(jj)
                                 + 2.0 * K3(jj)
                                 + K4(jj)) / 6.0;
            }
            t = timeInitial + ii * h;
            output(MatrixIndexBase + ii,MatrixIndexBase) = t;
            output(MatrixIndexBase + ii,_(MatrixIndexBase+1,MatrixIndexBase+numEqs)) = ~inputs; 
        }
        return output;
}
*/

Matrix RungeKuttaIntegrator::Integrate(const vector<ssfTime>& _propTime, odeFunctor _odeFunctorPtr, const Vector& _initialConditions, Orbit* _pOrbit, Attitude* _pAttitude, const Matrix& _constants, const Functor& _functorPtr)
{
        vector<ssfTime> propTime = _propTime; // need to create copy to prevent the stl .begin() and .end() from discarding the const-ness
        vector<ssfTime>::iterator timeIterator = propTime.begin();

        ssfTime timeInitial = (*timeIterator);
        ssfTime t = timeInitial;
        timeIterator = propTime.end() - 1;
        double h = ((*timeIterator).GetSeconds() - timeInitial.GetSeconds()) / m_NumSteps; // Calculation step size (seconds)
        int numEqs = _initialConditions.getIndexCount();
        Vector inputs = _initialConditions;
        Vector K1(numEqs);
        Vector K2(numEqs);
        Vector K3(numEqs);
        Vector K4(numEqs);

        Matrix output(m_NumSteps + 1, numEqs + 1);

        output(MatrixIndexBase,MatrixIndexBase) = timeInitial.GetSeconds();
        output(MatrixIndexBase,_(MatrixIndexBase+1,MatrixIndexBase+numEqs)) = ~_initialConditions;
        ssfTime tTemp;
        Vector temp(numEqs);
        for (int ii = 1; ii <= m_NumSteps; ++ii)
        {
            K1 = h * _odeFunctorPtr(t, inputs, _pOrbit, _pAttitude, _constants, _functorPtr);
            temp = inputs + K1 / 2.0;           tTemp = t + h/2;        
            K2 = h * _odeFunctorPtr(tTemp, temp, _pOrbit, _pAttitude, _constants, _functorPtr);
            temp = inputs + K2 / 2.0;                           
            K3 = h * _odeFunctorPtr(tTemp, temp, _pOrbit, _pAttitude, _constants, _functorPtr);
            temp = inputs + K3;         tTemp = t + h;          
            K4 = h * _odeFunctorPtr(tTemp, temp, _pOrbit, _pAttitude, _constants, _functorPtr);
            for (int jj = MatrixIndexBase; jj < MatrixIndexBase + numEqs; ++jj)
            {
                inputs(jj) += (K1(jj)
                                 + 2.0 * K2(jj)
                                 + 2.0 * K3(jj)
                                 + K4(jj)) / 6.0;
            }
            t = timeInitial + ii * h;
            output(MatrixIndexBase + ii,MatrixIndexBase) = t.GetSeconds();
            output(MatrixIndexBase + ii,_(MatrixIndexBase+1,MatrixIndexBase+numEqs)) = ~inputs; 
        }
        return output;
}
} // close namespace O_SESSAME

// Do not change the comments below - they will be added automatically by CVS
/*****************************************************************************
*       $Log: RungeKuttaIntegrator.cpp,v $
*       Revision 1.7  2003/05/22 02:59:15  nilspace
*       Updated comments. Changed to pass in pointers to Orbit & Attitude objects.
*       
*       Revision 1.6  2003/05/13 18:58:27  nilspace
*       Cleaned up comments.
*       
*       Revision 1.5  2003/04/27 22:04:34  nilspace
*       Created the namespace O_SESSAME.
*       
*       Revision 1.4  2003/04/25 14:15:06  nilspace
*       Created a temporary vectro<ssfTime> to prevent the STL .begin() and .end() functions from discarding the const-ness of _propTime in Integrate().
*       
*       Revision 1.3  2003/04/25 13:45:55  nilspace
*       const'd Get() functions.
*       
*       Revision 1.2  2003/04/23 16:30:59  nilspace
*       Various bugfixes & uploading of all changed code for new programmers.
*       
*       Revision 1.1  2003/04/08 22:29:55  nilspace
*       Initial submission. Made a subclass of Integrator.
*       
*       
*
******************************************************************************/

---