testEnvironment.cpp

An EnvFunction is used to store a function pointer and the list of parameters that are used in evaluating an environmental disturbance function. The function can either be a force or torque disturbance. Examples include:

(force) Gravity, aerodynamic drag, solar radiation pressure

(torque) Gravity gradient, magnetic, aerodynamic, fuel slosh <p>

The parameters are the constants (such as Ballistic coefficient, magnetic dipole, cross-sectional area), in the appropriate function order, and can be set or removed individually. Furthermore, because they are actually pointers to data, once the pointer is set, the data can still change (for example the mass of the spacecraft) without having to alter the EnvFunction parameter list.

Usage:

To use the EnvFunction class, the user first creates a disturbance function that matches the prototype of pt2EnvFunc, which takes the current time, orbit state, attitude state, and an STL vector of pointers to the parameters. They then must create an instance of an EnvFunction with a pointer to the function:

                EnvFunction DragForce(&myDisturbanceFunction);

and then fill out the parameter list:

                DragForce.AddParameter(reinterpret_cast<void*>(*CrossArea), 1);

Notice it is important to only give a pointer to a value that will stay "in scope" as long as needed for the simulation. Also, the reinterpret_cast<void*>() is needed to make the value 'look like' a void pointer. Don't worry too much about this, but read Stroustrop if you'd like to know more.

Example disturbance function:

// Force function that will be called each timestep
Vector DragForceFunction(const ssfTime* _currentTime, const OrbitState* _currentOrbitState, const AttitudeState* _currentAttitudeState, EnvFuncParamaterType _parameterList)
{
    Vector Forces(3);
    double BC = *(reinterpret_cast<double*>(_parameterList[1]));
    double Density = *(reinterpret_cast<double*>(_parameterList[2])); // kg/m^3
    Vector Vrel(3); Vrel = _currentOrbitState->GetState()(_(VectorIndexBase+3,VectorIndexBase+5));
    double Vrel_mag = norm2(Vrel);
    Forces = -1/2 * rho / BC * pow(Vrel_mag,2) * Vrel / Vrel_mag; 
    return Forces;
}

Example using an Environment object for integration.

testPropagation.cpp Example propagation of a combined satellite orbit & attitude.

/*  
*/
#include "orbitmodels/TwoBodyDynamics.h"
#include "Matrix.h"
#include "RungeKuttaIntegrator.h"
#include "PositionVelocity.h"
#include "Plot.h"
#include "Environment.h"
#include "OrbitStateRepresentation.h"
#include "EarthCentralBody.h"
#include "OrbitState.h"
#include "orbitframes/OrbitFrameIJK.h"

// Force function that will be called each timestep
Vector GravityForceFunction(ssfTime* _currentTime, OrbitState* _currentOrbitState, AttitudeState* _currentAttitudeState, EnvFuncParamaterType _parameterList)
{
    Vector state = _currentOrbitState->GetState();
    Vector Forces(3);
    Vector Position(3); Position(_) = state(_(VectorIndexBase,VectorIndexBase+2));
    Forces = -398600.441 / pow(norm2(Position),3) * Position;
    return Forces;
}

Vector DragForceFunction(ssfTime* _currentTime, OrbitState* _currentOrbitState, AttitudeState* _currentAttitudeState, EnvFuncParamaterType _parameterList)
{
    Vector Forces(3);
    double BC = *(reinterpret_cast<double*>(_parameterList[1]));
    Vector Vrel(3); Vrel = _currentOrbitState->GetState()(_(VectorIndexBase+3,VectorIndexBase+5));
    double Vrel_mag = norm2(Vrel);
    Forces = -1/2 *1 / BC * pow(Vrel_mag,2) * Vrel / Vrel_mag; 
    return Forces;
}

int main()
{
    // Setup an integrator and any special parameters
    RungeKuttaIntegrator myIntegrator; 
        int numOrbits = 10;
        int numSteps = 100;
  
    cout << "Number of Orbits: " << flush;
    cin >> numOrbits;
    cout << "Number of Integration Steps: "<< flush;
    cin >> numSteps;

    myIntegrator.SetNumSteps(numSteps);
    
    // Create & initialize an orbit type
    OrbitState myOrbit;
    myOrbit.SetStateRepresentation(new PositionVelocity);
    myOrbit.SetOrbitFrame(new OrbitFrameIJK);
    Vector initPV(6);
        initPV(VectorIndexBase+0) = -5776.6; // km/s
        initPV(VectorIndexBase+1) = -157; // km/s        
        initPV(VectorIndexBase+2) = 3496.9; // km/s    
        initPV(VectorIndexBase+3) = -2.595;  // km/s
        initPV(VectorIndexBase+4) = -5.651;  // km/s        
        initPV(VectorIndexBase+5) = -4.513; // km/s
    myOrbit.SetState(initPV);
    
    AttitudeState myAttitude;
    myAttitude.SetRotation(Rotation(Quaternion(0,0,0,1)));
    
    // Create the matrix of parameters needed for the RHS
    Matrix Parameters(1,1);
    Parameters(MatrixIndexBase,MatrixIndexBase) = 398600.441; //km / s^2

    
    // ENVIRONMENT TESTING
    Environment* pEarthEnv = new Environment;
    EarthCentralBody *pCBEarth = new EarthCentralBody;
    pEarthEnv->SetCentralBody(pCBEarth);
    cout << "Filling Parameters" << endl;
    EnvFunction TwoBodyGravity(&GravityForceFunction);
    pEarthEnv->AddForceFunction(TwoBodyGravity);
    EnvFunction DragForce(&DragForceFunction);
    double BC = 200;
    DragForce.AddParameter(reinterpret_cast<void*>(&BC), 1);
    pEarthEnv->AddForceFunction(DragForce);
    cout << "Finished Setting up Earth Environment" << endl;


    vector<ssfTime> integrationTimes;
    ssfTime begin(Now());
    ssfTime end(begin + 92*60*numOrbits);
    integrationTimes.push_back(begin);
    integrationTimes.push_back(end);
    cout << "PropTime = " << begin.GetSeconds() << " s -> " << end.GetSeconds() << " s" << endl;
    cout << "Orbit State: " << ~myOrbit.GetState() << endl;

    ObjectFunctor<Environment> OrbitForcesFunctor(pEarthEnv, &Environment::GetForces);
    // Integrate over the desired time interval
    tick();
    Matrix history = myIntegrator.Integrate(
                            integrationTimes,           // seconds
                            &TwoBodyDynamics, 
                            myOrbit.GetState(),
                            NULL,
                            NULL,
                            Parameters,
                            OrbitForcesFunctor
                        );

    cout << "finished propagating in " << tock() << " seconds." << endl;                                 
    // Output the flight history
    //cout << history;
    Matrix plotting = history(_,_(MatrixIndexBase+1,MatrixIndexBase+3));

    Plot3D(plotting);
    delete pEarthEnv;
    
    return 0;
                            
}

// Do not change the comments below - they will be added automatically by CVS
/*****************************************************************************
*       $Log: testEnvironment.cpp,v $
*       Revision 1.1  2003/04/08 23:06:26  nilspace
*       Initial Submission.
*       
*
******************************************************************************/

                        

---