Skip to content

Commit

Permalink
Added first iteration of one-way Doppler model
Browse files Browse the repository at this point in the history
  • Loading branch information
DominicDirkx committed Feb 26, 2017
1 parent cfac0a8 commit d99f13c
Show file tree
Hide file tree
Showing 4 changed files with 353 additions and 7 deletions.
4 changes: 4 additions & 0 deletions Tudat/Astrodynamics/ObservationModels/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,10 @@ if(USE_CSPICE)
setup_custom_test_program(test_OneWayRangeModel "${SRCROOT}${OBSERVATIONMODELSDIR}")
target_link_libraries(test_OneWayRangeModel ${TUDAT_ESTIMATION_LIBRARIES} ${Boost_LIBRARIES})

add_executable(test_OneWayDopplerModel "${SRCROOT}${OBSERVATIONMODELSDIR}/UnitTests/unitTestOneWayDopplerModel.cpp")
setup_custom_test_program(test_OneWayDopplerModel "${SRCROOT}${OBSERVATIONMODELSDIR}")
target_link_libraries(test_OneWayDopplerModel ${TUDAT_ESTIMATION_LIBRARIES} ${Boost_LIBRARIES})

add_executable(test_PositionObservationModel "${SRCROOT}${OBSERVATIONMODELSDIR}/UnitTests/unitTestPositionObservationModel.cpp")
setup_custom_test_program(test_PositionObservationModel "${SRCROOT}${OBSERVATIONMODELSDIR}")
target_link_libraries(test_PositionObservationModel ${TUDAT_ESTIMATION_LIBRARIES} ${Boost_LIBRARIES})
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,151 @@
/* Copyright (c) 2010-2017, Delft University of Technology
* All rigths reserved
*
* This file is part of the Tudat. Redistribution and use in source and
* binary forms, with or without modification, are permitted exclusively
* under the terms of the Modified BSD license. You should have received
* a copy of the license with this file. If not, please or visit:
* http://tudat.tudelft.nl/LICENSE.
*/

#define BOOST_TEST_MAIN

#include <limits>
#include <string>

#include <boost/format.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/make_shared.hpp>

#include "Tudat/Basics/testMacros.h"

#include "Tudat/InputOutput/basicInputOutput.h"

#include "Tudat/SimulationSetup/EnvironmentSetup/body.h"
#include "Tudat/Astrodynamics/ObservationModels/oneWayRangeObservationModel.h"
#include "Tudat/SimulationSetup/EstimationSetup/createObservationModel.h"
#include "Tudat/SimulationSetup/EnvironmentSetup/defaultBodies.h"
#include "Tudat/SimulationSetup/EnvironmentSetup/createBodies.h"

namespace tudat
{
namespace unit_tests
{

using namespace tudat::observation_models;
using namespace tudat::spice_interface;
using namespace tudat::ephemerides;
using namespace tudat::simulation_setup;


BOOST_AUTO_TEST_SUITE( test_one_way_range_model )


BOOST_AUTO_TEST_CASE( testOneWayRangeModel )
{
std::string kernelsPath = input_output::getSpiceKernelPath( );
spice_interface::loadSpiceKernelInTudat( kernelsPath + "de-403-masses.tpc");
spice_interface::loadSpiceKernelInTudat( kernelsPath + "naif0009.tls");
spice_interface::loadSpiceKernelInTudat( kernelsPath + "pck00009.tpc");
spice_interface::loadSpiceKernelInTudat( kernelsPath + "de421.bsp");

// Define bodies to use.
std::vector< std::string > bodiesToCreate;
bodiesToCreate.push_back( "Earth" );
bodiesToCreate.push_back( "Sun" );
bodiesToCreate.push_back( "Mars" );

// Specify initial time
double initialEphemerisTime = 0.0;
double finalEphemerisTime = initialEphemerisTime + 7.0 * 86400.0;
double maximumTimeStep = 3600.0;
double buffer = 10.0 * maximumTimeStep;

// Create bodies settings needed in simulation
std::map< std::string, boost::shared_ptr< BodySettings > > defaultBodySettings =
getDefaultBodySettings(
bodiesToCreate, initialEphemerisTime - buffer, finalEphemerisTime + buffer );

// Create bodies
NamedBodyMap bodyMap = createBodies( defaultBodySettings );

setGlobalFrameBodyEphemerides( bodyMap, "SSB", "ECLIPJ2000" );

// Define link ends for observations.
LinkEnds linkEnds;
linkEnds[ transmitter ] = std::make_pair( "Earth" , "" );
linkEnds[ receiver ] = std::make_pair( "Mars" , "" );

// Create observation model.
std::vector< std::string > lightTimePerturbingBodies = boost::assign::list_of( "Sun" );
std::vector< boost::shared_ptr< LightTimeCorrectionSettings > > lightTimeCorrectionSettings;


boost::shared_ptr< ObservationModel< 1, double, double> > observationModel =
ObservationModelCreator< 1, double, double>::createObservationModel(
oneWayDoppler, linkEnds, bodyMap );

boost::shared_ptr< OneWayDopplerObservationModel< double, double> > dopplerObservationModel =
boost::dynamic_pointer_cast< OneWayDopplerObservationModel< double, double> >( observationModel );


// Compute observation separately with two functions.
for( unsigned testCase = 0; testCase < 2; testCase++ )
{

double observationTime = ( finalEphemerisTime + initialEphemerisTime ) / 2.0;
std::vector< double > linkEndTimes;
std::vector< Eigen::Vector6d > linkEndStates;

LinkEndType referenceLinkEnd;
if( testCase == 0 )
{
referenceLinkEnd = transmitter;
}
else
{
referenceLinkEnd = receiver;
}

double dopplerObservable = observationModel->computeObservationsWithLinkEndData(
observationTime, referenceLinkEnd, linkEndTimes, linkEndStates )( 0 );

boost::shared_ptr< LightTimeCalculator< double, double > > lightTimeCalculator =
dopplerObservationModel->getLightTimeCalculator( );
Eigen::Vector6d transmitterState, receiverState;
double lightTime = lightTimeCalculator->calculateLightTimeWithLinkEndsStates(
receiverState, transmitterState, observationTime, testCase );

{
TUDAT_CHECK_MATRIX_CLOSE_FRACTION( receiverState, linkEndStates.at( 1 ), 1.0E-15 );
TUDAT_CHECK_MATRIX_CLOSE_FRACTION( transmitterState, linkEndStates.at( 0 ), 1.0E-15 );

if( testCase == 0 )
{
BOOST_CHECK_SMALL( std::fabs( observationTime - linkEndTimes.at( 0 ) ), 1.0E-15 );
BOOST_CHECK_SMALL( std::fabs( observationTime + lightTime - linkEndTimes.at( 1 ) ), 1.0E-15 );
}
else
{
BOOST_CHECK_SMALL( std::fabs( observationTime - linkEndTimes.at( 1 ) ), 1.0E-15 );
BOOST_CHECK_SMALL( std::fabs( observationTime - lightTime - linkEndTimes.at( 0 ) ), 1.0E-15 );
}
}

double timePerturbation = 100.0;
double upPerturbedLightTime = lightTimeCalculator->calculateLightTime( linkEndTimes.at( 0 ) + timePerturbation, false );
double downPerturbedLightTime = lightTimeCalculator->calculateLightTime( linkEndTimes.at( 0 ) - timePerturbation, false );

double lightTimeSensitivity = ( upPerturbedLightTime - downPerturbedLightTime ) / ( 2.0 * timePerturbation );

BOOST_CHECK_SMALL( std::fabs( lightTimeSensitivity - dopplerObservable ), 1.0E-14 );
}
}

BOOST_AUTO_TEST_SUITE_END( )

}

}


Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,34 @@ ObservationScalarType calculateLineOfSightVelocityAsCFraction(
return lineOfSightUnitVector.dot( velocityVector ) / physical_constants::getSpeedOfLight< ObservationScalarType >( );
}

template< typename ObservationScalarType = double >
ObservationScalarType computeOneWayFirstOrderDopplerTaylorSeriesExpansion(
Eigen::Matrix< ObservationScalarType, 6, 1 >& transmitterState,
Eigen::Matrix< ObservationScalarType, 6, 1 >& receiverState,
const int taylorSeriesOrder )
{
Eigen::Matrix< ObservationScalarType, 3, 1 > relativePostion = ( ( receiverState - transmitterState ).segment( 0, 3 ) ).normalized( );

ObservationScalarType transmitterTerm =
( relativePostion.dot( transmitterState.segment( 3, 3 ) ) / physical_constants::getSpeedOfLight< ObservationScalarType >( ) );
ObservationScalarType receiverTerm =
( relativePostion.dot( receiverState.segment( 3, 3 ) ) / physical_constants::getSpeedOfLight< ObservationScalarType >( ) );
ObservationScalarType currentTaylorSeriesTerm =
mathematical_constants::getFloatingInteger< ObservationScalarType >( 1 );
ObservationScalarType currentTaylorSeries = mathematical_constants::getFloatingInteger< ObservationScalarType >( 0 );

for( int i = 0; i < taylorSeriesOrder; i++ )
{
currentTaylorSeriesTerm *= receiverTerm;
currentTaylorSeries += currentTaylorSeriesTerm;
}

return -transmitterTerm + currentTaylorSeries *
( mathematical_constants::getFloatingInteger< ObservationScalarType >( 1 ) - transmitterTerm );

}

//! Computes observable as d f_{A}/d_f_{B} - 1, with a the transmitter and B the receiver, omitting proper time variations and light time corrections
template< typename ObservationScalarType = double, typename TimeType = double >
class OneWayDopplerObservationModel: public ObservationModel< 1, ObservationScalarType, TimeType >
{
Expand All @@ -50,6 +78,7 @@ class OneWayDopplerObservationModel: public ObservationModel< 1, ObservationScal
lightTimeCalculator_( lightTimeCalculator )
{
one_ = mathematical_constants::getFloatingInteger< ObservationScalarType >( 1 );
taylorSeriesExpansionOrder_ = 3;
}

~OneWayDopplerObservationModel( ){ }
Expand Down Expand Up @@ -91,11 +120,11 @@ class OneWayDopplerObservationModel: public ObservationModel< 1, ObservationScal
lightTime = lightTimeCalculator_->calculateLightTimeWithLinkEndsStates(
receiverState, transmitterState, time, true );

PositionType relativePostion = ( ( receiverState - transmitterState ).segment( 0, 3 ) ).normalized( );

return ( Eigen::Matrix< ObservationScalarType, 1, 1 >( )
<<( one_ - calculateLineOfSightVelocityAsCFraction< ObservationScalarType >( relativePostion, receiverState.segment( 3, 3 ) ) ) /
( one_ - calculateLineOfSightVelocityAsCFraction< ObservationScalarType >( relativePostion, transmitterState.segment( 3, 3 ) ) ) ).finished( );
<< computeOneWayFirstOrderDopplerTaylorSeriesExpansion( transmitterState, receiverState, taylorSeriesExpansionOrder_ ) ).finished( );
// return ( Eigen::Matrix< ObservationScalarType, 1, 1 >( )
// <<( one_ - calculateLineOfSightVelocityAsCFraction< ObservationScalarType >( relativePostion, receiverState.segment( 3, 3 ) ) ) /
// ( one_ - calculateLineOfSightVelocityAsCFraction< ObservationScalarType >( relativePostion, transmitterState.segment( 3, 3 ) ) ) ).finished( );
}

//! Function to compute one-way Doppler observable without any corrections.
Expand Down Expand Up @@ -156,9 +185,12 @@ class OneWayDopplerObservationModel: public ObservationModel< 1, ObservationScal
PositionType relativePostion = ( ( receiverState - transmitterState ).segment( 0, 3 ) ).normalized( );
relativePostion /= relativePostion.norm( );

return ( Eigen::Matrix< ObservationScalarType, 1, 1 >( )
<<( one_ - calculateLineOfSightVelocityAsCFraction< ObservationScalarType >( relativePostion, receiverState.segment( 3, 3 ) ) ) /
( one_ - calculateLineOfSightVelocityAsCFraction< ObservationScalarType >( relativePostion, transmitterState.segment( 3, 3 ) ) ) ).finished( );
ObservationScalarType dopplerObservable = computeOneWayFirstOrderDopplerTaylorSeriesExpansion< ObservationScalarType >
( transmitterState, receiverState, taylorSeriesExpansionOrder_ );


return ( Eigen::Matrix< ObservationScalarType, 1, 1 >( ) << dopplerObservable ).finished( );


}

Expand All @@ -178,6 +210,8 @@ class OneWayDopplerObservationModel: public ObservationModel< 1, ObservationScal

ObservationScalarType one_;

int taylorSeriesExpansionOrder_;

};

}
Expand Down
Loading

0 comments on commit d99f13c

Please sign in to comment.