From 1a6113968cad7f40ccb2657414c498c9a8e68e7f Mon Sep 17 00:00:00 2001 From: Enne Hekma Date: Tue, 8 Mar 2016 17:47:26 +0100 Subject: [PATCH 1/7] Add departure epoch loop to lambertScanner --- src/lambertScanner.cpp | 289 +++++++++++++++++++++-------------------- 1 file changed, 148 insertions(+), 141 deletions(-) diff --git a/src/lambertScanner.cpp b/src/lambertScanner.cpp index d1a9d9b..4c2dbf7 100644 --- a/src/lambertScanner.cpp +++ b/src/lambertScanner.cpp @@ -205,150 +205,157 @@ void executeLambertScanner( const rapidjson::Document& config ) SGP4 sgp4Arrival( arrivalObject ); const int arrivalObjectId = static_cast< int >( arrivalObject.NoradNumber( ) ); - // Loop over time-of-flight grid. - for ( int k = 0; k < input.timeOfFlightSteps; k++ ) + // Loop over departure epochs + for (int l = 0; l < 5; ++l) { - const double timeOfFlight - = input.timeOfFlightMinimum + k * input.timeOfFlightStepSize; - - const DateTime arrivalEpoch = departureEpoch.AddSeconds( timeOfFlight ); - const Eci tleArrivalState = sgp4Arrival.FindPosition( arrivalEpoch ); - const Vector6 arrivalState = getStateVector( tleArrivalState ); - - Vector3 arrivalPosition; - std::copy( arrivalState.begin( ), - arrivalState.begin( ) + 3, - arrivalPosition.begin( ) ); - - Vector3 arrivalVelocity; - std::copy( arrivalState.begin( ) + 3, - arrivalState.end( ), - arrivalVelocity.begin( ) ); - const Vector6 arrivalStateKepler - = astro::convertCartesianToKeplerianElements( arrivalState, - earthGravitationalParameter ); - - kep_toolbox::lambert_problem targeter( departurePosition, - arrivalPosition, - timeOfFlight, - earthGravitationalParameter, - !input.isPrograde, - input.revolutionsMaximum ); - - const int numberOfSolutions = targeter.get_v1( ).size( ); - - // Compute Delta-Vs for transfer and determine index of lowest. - typedef std::vector< Vector3 > VelocityList; - VelocityList departureDeltaVs( numberOfSolutions ); - VelocityList arrivalDeltaVs( numberOfSolutions ); - - typedef std::vector< double > TransferDeltaVList; - TransferDeltaVList transferDeltaVs( numberOfSolutions ); - - for ( int i = 0; i < numberOfSolutions; i++ ) + DateTime departureEpoch = input.departureEpoch; + departureEpoch = departureEpoch.AddSeconds( 60*l ); // variing departure time + + // Loop over time-of-flight grid. + for ( int k = 0; k < input.timeOfFlightSteps; k++ ) { - // Compute Delta-V for transfer. - const Vector3 transferDepartureVelocity = targeter.get_v1( )[ i ]; - const Vector3 transferArrivalVelocity = targeter.get_v2( )[ i ]; - - departureDeltaVs[ i ] = sml::add( transferDepartureVelocity, - sml::multiply( departureVelocity, -1.0 ) ); - arrivalDeltaVs[ i ] = sml::add( transferArrivalVelocity, - sml::multiply( arrivalVelocity, -1.0 ) ); - - transferDeltaVs[ i ] - = sml::norm< double >( departureDeltaVs[ i ] ) - + sml::norm< double >( arrivalDeltaVs[ i ] ); + const double timeOfFlight + = input.timeOfFlightMinimum + k * input.timeOfFlightStepSize; + + const DateTime arrivalEpoch = departureEpoch.AddSeconds( timeOfFlight ); + const Eci tleArrivalState = sgp4Arrival.FindPosition( arrivalEpoch ); + const Vector6 arrivalState = getStateVector( tleArrivalState ); + + Vector3 arrivalPosition; + std::copy( arrivalState.begin( ), + arrivalState.begin( ) + 3, + arrivalPosition.begin( ) ); + + Vector3 arrivalVelocity; + std::copy( arrivalState.begin( ) + 3, + arrivalState.end( ), + arrivalVelocity.begin( ) ); + const Vector6 arrivalStateKepler + = astro::convertCartesianToKeplerianElements( arrivalState, + earthGravitationalParameter ); + + kep_toolbox::lambert_problem targeter( departurePosition, + arrivalPosition, + timeOfFlight, + earthGravitationalParameter, + !input.isPrograde, + input.revolutionsMaximum ); + + const int numberOfSolutions = targeter.get_v1( ).size( ); + + // Compute Delta-Vs for transfer and determine index of lowest. + typedef std::vector< Vector3 > VelocityList; + VelocityList departureDeltaVs( numberOfSolutions ); + VelocityList arrivalDeltaVs( numberOfSolutions ); + + typedef std::vector< double > TransferDeltaVList; + TransferDeltaVList transferDeltaVs( numberOfSolutions ); + + for ( int i = 0; i < numberOfSolutions; i++ ) + { + // Compute Delta-V for transfer. + const Vector3 transferDepartureVelocity = targeter.get_v1( )[ i ]; + const Vector3 transferArrivalVelocity = targeter.get_v2( )[ i ]; + + departureDeltaVs[ i ] = sml::add( transferDepartureVelocity, + sml::multiply( departureVelocity, -1.0 ) ); + arrivalDeltaVs[ i ] = sml::add( transferArrivalVelocity, + sml::multiply( arrivalVelocity, -1.0 ) ); + + transferDeltaVs[ i ] + = sml::norm< double >( departureDeltaVs[ i ] ) + + sml::norm< double >( arrivalDeltaVs[ i ] ); + } + + const TransferDeltaVList::iterator minimumDeltaVIterator + = std::min_element( transferDeltaVs.begin( ), transferDeltaVs.end( ) ); + const int minimumDeltaVIndex + = std::distance( transferDeltaVs.begin( ), minimumDeltaVIterator ); + + const int revolutions = std::floor( ( minimumDeltaVIndex + 1 ) / 2 ); + + Vector6 transferState; + std::copy( departurePosition.begin( ), + departurePosition.begin( ) + 3, + transferState.begin( ) ); + std::copy( targeter.get_v1( )[ minimumDeltaVIndex ].begin( ), + targeter.get_v1( )[ minimumDeltaVIndex ].begin( ) + 3, + transferState.begin( ) + 3 ); + + const Vector6 transferStateKepler + = astro::convertCartesianToKeplerianElements( transferState, + earthGravitationalParameter ); + + // Bind values to SQL insert query. + query.bind( ":departure_object_id", departureObjectId ); + query.bind( ":arrival_object_id", arrivalObjectId ); + query.bind( ":departure_epoch", departureEpoch.ToJulian( ) ); + query.bind( ":time_of_flight", timeOfFlight ); + query.bind( ":revolutions", revolutions ); + query.bind( ":prograde", input.isPrograde ); + query.bind( ":departure_position_x", departureState[ astro::xPositionIndex ] ); + query.bind( ":departure_position_y", departureState[ astro::yPositionIndex ] ); + query.bind( ":departure_position_z", departureState[ astro::zPositionIndex ] ); + query.bind( ":departure_velocity_x", departureState[ astro::xVelocityIndex ] ); + query.bind( ":departure_velocity_y", departureState[ astro::yVelocityIndex ] ); + query.bind( ":departure_velocity_z", departureState[ astro::zVelocityIndex ] ); + query.bind( ":departure_semi_major_axis", + departureStateKepler[ astro::semiMajorAxisIndex ] ); + query.bind( ":departure_eccentricity", + departureStateKepler[ astro::eccentricityIndex ] ); + query.bind( ":departure_inclination", + departureStateKepler[ astro::inclinationIndex ] ); + query.bind( ":departure_argument_of_periapsis", + departureStateKepler[ astro::argumentOfPeriapsisIndex ] ); + query.bind( ":departure_longitude_of_ascending_node", + departureStateKepler[ astro::longitudeOfAscendingNodeIndex ] ); + query.bind( ":departure_true_anomaly", + departureStateKepler[ astro::trueAnomalyIndex ] ); + query.bind( ":arrival_position_x", arrivalState[ astro::xPositionIndex ] ); + query.bind( ":arrival_position_y", arrivalState[ astro::yPositionIndex ] ); + query.bind( ":arrival_position_z", arrivalState[ astro::zPositionIndex ] ); + query.bind( ":arrival_velocity_x", arrivalState[ astro::xVelocityIndex ] ); + query.bind( ":arrival_velocity_y", arrivalState[ astro::yVelocityIndex ] ); + query.bind( ":arrival_velocity_z", arrivalState[ astro::zVelocityIndex ] ); + query.bind( ":arrival_semi_major_axis", + arrivalStateKepler[ astro::semiMajorAxisIndex ] ); + query.bind( ":arrival_eccentricity", + arrivalStateKepler[ astro::eccentricityIndex ] ); + query.bind( ":arrival_inclination", + arrivalStateKepler[ astro::inclinationIndex ] ); + query.bind( ":arrival_argument_of_periapsis", + arrivalStateKepler[ astro::argumentOfPeriapsisIndex ] ); + query.bind( ":arrival_longitude_of_ascending_node", + arrivalStateKepler[ astro::longitudeOfAscendingNodeIndex ] ); + query.bind( ":arrival_true_anomaly", + arrivalStateKepler[ astro::trueAnomalyIndex ] ); + query.bind( ":transfer_semi_major_axis", + transferStateKepler[ astro::semiMajorAxisIndex ] ); + query.bind( ":transfer_eccentricity", + transferStateKepler[ astro::eccentricityIndex ] ); + query.bind( ":transfer_inclination", + transferStateKepler[ astro::inclinationIndex ] ); + query.bind( ":transfer_argument_of_periapsis", + transferStateKepler[ astro::argumentOfPeriapsisIndex ] ); + query.bind( ":transfer_longitude_of_ascending_node", + transferStateKepler[ astro::longitudeOfAscendingNodeIndex ] ); + query.bind( ":transfer_true_anomaly", + transferStateKepler[ astro::trueAnomalyIndex ] ); + query.bind( ":departure_delta_v_x", departureDeltaVs[ minimumDeltaVIndex ][ 0 ] ); + query.bind( ":departure_delta_v_y", departureDeltaVs[ minimumDeltaVIndex ][ 1 ] ); + query.bind( ":departure_delta_v_z", departureDeltaVs[ minimumDeltaVIndex ][ 2 ] ); + query.bind( ":arrival_delta_v_x", arrivalDeltaVs[ minimumDeltaVIndex ][ 0 ] ); + query.bind( ":arrival_delta_v_y", arrivalDeltaVs[ minimumDeltaVIndex ][ 1 ] ); + query.bind( ":arrival_delta_v_z", arrivalDeltaVs[ minimumDeltaVIndex ][ 2 ] ); + query.bind( ":transfer_delta_v", *minimumDeltaVIterator ); + + // Execute insert query. + query.executeStep( ); + + // Reset SQL insert query. + query.reset( ); } - - const TransferDeltaVList::iterator minimumDeltaVIterator - = std::min_element( transferDeltaVs.begin( ), transferDeltaVs.end( ) ); - const int minimumDeltaVIndex - = std::distance( transferDeltaVs.begin( ), minimumDeltaVIterator ); - - const int revolutions = std::floor( ( minimumDeltaVIndex + 1 ) / 2 ); - - Vector6 transferState; - std::copy( departurePosition.begin( ), - departurePosition.begin( ) + 3, - transferState.begin( ) ); - std::copy( targeter.get_v1( )[ minimumDeltaVIndex ].begin( ), - targeter.get_v1( )[ minimumDeltaVIndex ].begin( ) + 3, - transferState.begin( ) + 3 ); - - const Vector6 transferStateKepler - = astro::convertCartesianToKeplerianElements( transferState, - earthGravitationalParameter ); - - // Bind values to SQL insert query. - query.bind( ":departure_object_id", departureObjectId ); - query.bind( ":arrival_object_id", arrivalObjectId ); - query.bind( ":departure_epoch", departureEpoch.ToJulian( ) ); - query.bind( ":time_of_flight", timeOfFlight ); - query.bind( ":revolutions", revolutions ); - query.bind( ":prograde", input.isPrograde ); - query.bind( ":departure_position_x", departureState[ astro::xPositionIndex ] ); - query.bind( ":departure_position_y", departureState[ astro::yPositionIndex ] ); - query.bind( ":departure_position_z", departureState[ astro::zPositionIndex ] ); - query.bind( ":departure_velocity_x", departureState[ astro::xVelocityIndex ] ); - query.bind( ":departure_velocity_y", departureState[ astro::yVelocityIndex ] ); - query.bind( ":departure_velocity_z", departureState[ astro::zVelocityIndex ] ); - query.bind( ":departure_semi_major_axis", - departureStateKepler[ astro::semiMajorAxisIndex ] ); - query.bind( ":departure_eccentricity", - departureStateKepler[ astro::eccentricityIndex ] ); - query.bind( ":departure_inclination", - departureStateKepler[ astro::inclinationIndex ] ); - query.bind( ":departure_argument_of_periapsis", - departureStateKepler[ astro::argumentOfPeriapsisIndex ] ); - query.bind( ":departure_longitude_of_ascending_node", - departureStateKepler[ astro::longitudeOfAscendingNodeIndex ] ); - query.bind( ":departure_true_anomaly", - departureStateKepler[ astro::trueAnomalyIndex ] ); - query.bind( ":arrival_position_x", arrivalState[ astro::xPositionIndex ] ); - query.bind( ":arrival_position_y", arrivalState[ astro::yPositionIndex ] ); - query.bind( ":arrival_position_z", arrivalState[ astro::zPositionIndex ] ); - query.bind( ":arrival_velocity_x", arrivalState[ astro::xVelocityIndex ] ); - query.bind( ":arrival_velocity_y", arrivalState[ astro::yVelocityIndex ] ); - query.bind( ":arrival_velocity_z", arrivalState[ astro::zVelocityIndex ] ); - query.bind( ":arrival_semi_major_axis", - arrivalStateKepler[ astro::semiMajorAxisIndex ] ); - query.bind( ":arrival_eccentricity", - arrivalStateKepler[ astro::eccentricityIndex ] ); - query.bind( ":arrival_inclination", - arrivalStateKepler[ astro::inclinationIndex ] ); - query.bind( ":arrival_argument_of_periapsis", - arrivalStateKepler[ astro::argumentOfPeriapsisIndex ] ); - query.bind( ":arrival_longitude_of_ascending_node", - arrivalStateKepler[ astro::longitudeOfAscendingNodeIndex ] ); - query.bind( ":arrival_true_anomaly", - arrivalStateKepler[ astro::trueAnomalyIndex ] ); - query.bind( ":transfer_semi_major_axis", - transferStateKepler[ astro::semiMajorAxisIndex ] ); - query.bind( ":transfer_eccentricity", - transferStateKepler[ astro::eccentricityIndex ] ); - query.bind( ":transfer_inclination", - transferStateKepler[ astro::inclinationIndex ] ); - query.bind( ":transfer_argument_of_periapsis", - transferStateKepler[ astro::argumentOfPeriapsisIndex ] ); - query.bind( ":transfer_longitude_of_ascending_node", - transferStateKepler[ astro::longitudeOfAscendingNodeIndex ] ); - query.bind( ":transfer_true_anomaly", - transferStateKepler[ astro::trueAnomalyIndex ] ); - query.bind( ":departure_delta_v_x", departureDeltaVs[ minimumDeltaVIndex ][ 0 ] ); - query.bind( ":departure_delta_v_y", departureDeltaVs[ minimumDeltaVIndex ][ 1 ] ); - query.bind( ":departure_delta_v_z", departureDeltaVs[ minimumDeltaVIndex ][ 2 ] ); - query.bind( ":arrival_delta_v_x", arrivalDeltaVs[ minimumDeltaVIndex ][ 0 ] ); - query.bind( ":arrival_delta_v_y", arrivalDeltaVs[ minimumDeltaVIndex ][ 1 ] ); - query.bind( ":arrival_delta_v_z", arrivalDeltaVs[ minimumDeltaVIndex ][ 2 ] ); - query.bind( ":transfer_delta_v", *minimumDeltaVIterator ); - - // Execute insert query. - query.executeStep( ); - - // Reset SQL insert query. - query.reset( ); - } + } } ++showProgress; From 1798b7f6ff41168764e84b9ccbd022d406e39961 Mon Sep 17 00:00:00 2001 From: Enne Hekma Date: Wed, 9 Mar 2016 15:05:04 +0100 Subject: [PATCH 2/7] Add departure epoch loop to lambert_scanner mode. Associated test has to be updated. --- config/lambert_scanner.json.empty | 4 +++ include/D2D/lambertScanner.hpp | 42 ++++++++++++++-------- src/lambertScanner.cpp | 59 ++++++++++++++++++------------- 3 files changed, 65 insertions(+), 40 deletions(-) diff --git a/config/lambert_scanner.json.empty b/config/lambert_scanner.json.empty index 8c0e88d..18fc8d8 100644 --- a/config/lambert_scanner.json.empty +++ b/config/lambert_scanner.json.empty @@ -13,6 +13,7 @@ // WARNING: if the database file already exists, it will be overwritten! "database" : "../data/test_lambert_scanner_catalog.db", + // Set departure epoch for transfers (common to all transfers computed). // Format: [year,month,day,hours,minutes,seconds] (all elements in the array must be integers)/ // The elements year, month and day are required. The others are optional but must be included @@ -20,6 +21,9 @@ // If departure_epoch array is left empty, the TLE epoch for the departure object is assumed. "departure_epoch" : [], + // Set departure epoch grid: [range (s), # of steps] + "departure_epoch_grid" : [,], + // Set time-of-flight grid: [min (s), max (s), # of steps]. "time_of_flight_grid" : [,,], diff --git a/include/D2D/lambertScanner.hpp b/include/D2D/lambertScanner.hpp index 7321ffe..bd2dd6f 100644 --- a/include/D2D/lambertScanner.hpp +++ b/include/D2D/lambertScanner.hpp @@ -52,22 +52,26 @@ struct LambertScannerInput * Constructs data struct based on verified input parameters. * * @sa checkLambertScannerInput, executeLambertScanner - * @param[in] aCatalogPath Path to TLE catalog - * @param[in] aDatabasePath Path to SQLite database - * @param[in] aDepartureEpoch Departure epoch for all transfers - * @param[in] aTimeOfFlightMinimum Minimum time-of-flight [s] - * @param[in] aTimeOfFlightMaximum Maximum time-of-flight [s] - * @param[in] someTimeOfFlightSteps Number of steps to take in time-of-flight grid - * @param[in] aTimeOfFlightStepSize Time-of-flight step size (derived parameter) [s] - * @param[in] progradeFlag Flag indicating if prograde transfer should be computed - * (false = retrograde) - * @param[in] aRevolutionsMaximum Maximum number of revolutions - * @param[in] aShortlistLength Number of transfers to include in shortlist - * @param[in] aShortlistPath Path to shortlist file + * @param[in] aCatalogPath Path to TLE catalog + * @param[in] aDatabasePath Path to SQLite database + * @param[in] aDepartureEpochInitial Departure epoch grid initial epoch + * @param[in] aDepartureEpochSteps Number of steps to take in departure epoch grid + * @param[in] aDepartureEpochStepSize Departure epoch grid step size (derived parameter) [s] + * @param[in] aTimeOfFlightMinimum Minimum time-of-flight [s] + * @param[in] aTimeOfFlightMaximum Maximum time-of-flight [s] + * @param[in] someTimeOfFlightSteps Number of steps to take in time-of-flight grid + * @param[in] aTimeOfFlightStepSize Time-of-flight step size (derived parameter) [s] + * @param[in] progradeFlag Flag indicating if prograde transfer should be computed + * (false = retrograde) + * @param[in] aRevolutionsMaximum Maximum number of revolutions + * @param[in] aShortlistLength Number of transfers to include in shortlist + * @param[in] aShortlistPath Path to shortlist file */ LambertScannerInput( const std::string& aCatalogPath, const std::string& aDatabasePath, - const DateTime& aDepartureEpoch, + const DateTime& aDepartureEpochInitial, + const double aDepartureEpochSteps, + const double aDepartureEpochStepSize, const double aTimeOfFlightMinimum, const double aTimeOfFlightMaximum, const double someTimeOfFlightSteps, @@ -78,7 +82,9 @@ struct LambertScannerInput const std::string& aShortlistPath ) : catalogPath( aCatalogPath ), databasePath( aDatabasePath ), - departureEpoch( aDepartureEpoch ), + departureEpochInitial( aDepartureEpochInitial ), + departureEpochSteps( aDepartureEpochSteps ), + departureEpochStepSize( aDepartureEpochStepSize ), timeOfFlightMinimum( aTimeOfFlightMinimum ), timeOfFlightMaximum( aTimeOfFlightMaximum ), timeOfFlightSteps( someTimeOfFlightSteps ), @@ -96,7 +102,13 @@ struct LambertScannerInput const std::string databasePath; //! Departure epoch. - const DateTime departureEpoch; + const DateTime departureEpochInitial; + + //! Number of departure epoch steps. + const double departureEpochSteps; + + //! Departure epoch grid step size. + const double departureEpochStepSize; //! Minimum time-of-flight [s]. const double timeOfFlightMinimum; diff --git a/src/lambertScanner.cpp b/src/lambertScanner.cpp index 4c2dbf7..161e6fc 100644 --- a/src/lambertScanner.cpp +++ b/src/lambertScanner.cpp @@ -33,7 +33,6 @@ namespace d2d //! Execute lambert_scanner. void executeLambertScanner( const rapidjson::Document& config ) { - // Verify config parameters. Exception is thrown if any of the parameters are missing. const LambertScannerInput input = checkLambertScannerInput( config ); @@ -41,7 +40,7 @@ void executeLambertScanner( const rapidjson::Document& config ) const double earthGravitationalParameter = kMU; std::cout << "Earth gravitational parameter " << earthGravitationalParameter << " kg m^3 s^-2" << std::endl; - + std::cout << std::endl; std::cout << "******************************************************************" << std::endl; std::cout << " Simulation & Output " << std::endl; @@ -168,29 +167,12 @@ void executeLambertScanner( const rapidjson::Document& config ) Tle departureObject = tleObjects[ i ]; SGP4 sgp4Departure( departureObject ); - DateTime departureEpoch = input.departureEpoch; - if ( input.departureEpoch == DateTime( ) ) + DateTime departureEpoch = input.departureEpochInitial; + if ( input.departureEpochInitial == DateTime( ) ) { departureEpoch = departureObject.Epoch( ); } - const Eci tleDepartureState = sgp4Departure.FindPosition( departureEpoch ); - const Vector6 departureState = getStateVector( tleDepartureState ); - - Vector3 departurePosition; - std::copy( departureState.begin( ), - departureState.begin( ) + 3, - departurePosition.begin( ) ); - - Vector3 departureVelocity; - std::copy( departureState.begin( ) + 3, - departureState.end( ), - departureVelocity.begin( ) ); - - const Vector6 departureStateKepler - = astro::convertCartesianToKeplerianElements( departureState, - earthGravitationalParameter ); - const int departureObjectId = static_cast< int >( departureObject.NoradNumber( ) ); // Loop over arrival objects. for ( unsigned int j = 0; j < tleObjects.size( ); j++ ) @@ -205,11 +187,29 @@ void executeLambertScanner( const rapidjson::Document& config ) SGP4 sgp4Arrival( arrivalObject ); const int arrivalObjectId = static_cast< int >( arrivalObject.NoradNumber( ) ); - // Loop over departure epochs - for (int l = 0; l < 5; ++l) + // Loop over departure epoch grid. + for (int l = 0; l < input.departureEpochSteps; ++l) { - DateTime departureEpoch = input.departureEpoch; - departureEpoch = departureEpoch.AddSeconds( 60*l ); // variing departure time + DateTime departureEpoch = input.departureEpochInitial; + departureEpoch = departureEpoch.AddSeconds( input.departureEpochStepSize * l ); + + const Eci tleDepartureState = sgp4Departure.FindPosition( departureEpoch ); + const Vector6 departureState = getStateVector( tleDepartureState ); + + Vector3 departurePosition; + std::copy( departureState.begin( ), + departureState.begin( ) + 3, + departurePosition.begin( ) ); + + Vector3 departureVelocity; + std::copy( departureState.begin( ) + 3, + departureState.end( ), + departureVelocity.begin( ) ); + + const Vector6 departureStateKepler + = astro::convertCartesianToKeplerianElements( departureState, + earthGravitationalParameter ); + const int departureObjectId = static_cast< int >( departureObject.NoradNumber( ) ); // Loop over time-of-flight grid. for ( int k = 0; k < input.timeOfFlightSteps; k++ ) @@ -439,6 +439,13 @@ LambertScannerInput checkLambertScannerInput( const rapidjson::Document& config std::cout << "Departure epoch " << departureEpoch << std::endl; } + const double departureEpochRange + = find( config, "departure_epoch_grid" )->value[ 0 ].GetDouble( ); + std::cout << "Departure epoch grid range " << departureEpochRange << std::endl; + const double departureGridSteps + = find( config, "departure_epoch_grid" )->value[ 1 ].GetDouble( ); + std::cout << "Departure epoch grid steps " << departureGridSteps << std::endl; + const double timeOfFlightMinimum = find( config, "time_of_flight_grid" )->value[ 0 ].GetDouble( ); std::cout << "Minimum Time-of-Flight " << timeOfFlightMinimum << std::endl; @@ -481,6 +488,8 @@ LambertScannerInput checkLambertScannerInput( const rapidjson::Document& config return LambertScannerInput( catalogPath, databasePath, departureEpoch, + departureGridSteps, + departureEpochRange/departureGridSteps, timeOfFlightMinimum, timeOfFlightMaximum, timeOfFlightSteps, From cda34a79596ee2844ed62ad630bc847fe4193a85 Mon Sep 17 00:00:00 2001 From: Enne Hekma Date: Wed, 9 Mar 2016 15:25:35 +0100 Subject: [PATCH 3/7] Correct blank lines --- src/lambertScanner.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/lambertScanner.cpp b/src/lambertScanner.cpp index 161e6fc..0fd4158 100644 --- a/src/lambertScanner.cpp +++ b/src/lambertScanner.cpp @@ -33,6 +33,7 @@ namespace d2d //! Execute lambert_scanner. void executeLambertScanner( const rapidjson::Document& config ) { + // Verify config parameters. Exception is thrown if any of the parameters are missing. const LambertScannerInput input = checkLambertScannerInput( config ); @@ -40,7 +41,7 @@ void executeLambertScanner( const rapidjson::Document& config ) const double earthGravitationalParameter = kMU; std::cout << "Earth gravitational parameter " << earthGravitationalParameter << " kg m^3 s^-2" << std::endl; - + std::cout << std::endl; std::cout << "******************************************************************" << std::endl; std::cout << " Simulation & Output " << std::endl; From c703fe84e9eb16eebe51644ddd32380fa8e771cf Mon Sep 17 00:00:00 2001 From: Enne Hekma Date: Wed, 9 Mar 2016 15:31:37 +0100 Subject: [PATCH 4/7] Correct blank lines --- config/lambert_scanner.json.empty | 1 - 1 file changed, 1 deletion(-) diff --git a/config/lambert_scanner.json.empty b/config/lambert_scanner.json.empty index 18fc8d8..224029f 100644 --- a/config/lambert_scanner.json.empty +++ b/config/lambert_scanner.json.empty @@ -13,7 +13,6 @@ // WARNING: if the database file already exists, it will be overwritten! "database" : "../data/test_lambert_scanner_catalog.db", - // Set departure epoch for transfers (common to all transfers computed). // Format: [year,month,day,hours,minutes,seconds] (all elements in the array must be integers)/ // The elements year, month and day are required. The others are optional but must be included From d11d5ee5077658a49272aba51b85c08d3ea1026a Mon Sep 17 00:00:00 2001 From: Enne Hekma Date: Wed, 9 Mar 2016 16:33:24 +0100 Subject: [PATCH 5/7] Change variable naming to be consistent with the rest --- include/D2D/lambertScanner.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/D2D/lambertScanner.hpp b/include/D2D/lambertScanner.hpp index bd2dd6f..6f1218e 100644 --- a/include/D2D/lambertScanner.hpp +++ b/include/D2D/lambertScanner.hpp @@ -55,7 +55,7 @@ struct LambertScannerInput * @param[in] aCatalogPath Path to TLE catalog * @param[in] aDatabasePath Path to SQLite database * @param[in] aDepartureEpochInitial Departure epoch grid initial epoch - * @param[in] aDepartureEpochSteps Number of steps to take in departure epoch grid + * @param[in] someDepartureEpochSteps Number of steps to take in departure epoch grid * @param[in] aDepartureEpochStepSize Departure epoch grid step size (derived parameter) [s] * @param[in] aTimeOfFlightMinimum Minimum time-of-flight [s] * @param[in] aTimeOfFlightMaximum Maximum time-of-flight [s] @@ -70,7 +70,7 @@ struct LambertScannerInput LambertScannerInput( const std::string& aCatalogPath, const std::string& aDatabasePath, const DateTime& aDepartureEpochInitial, - const double aDepartureEpochSteps, + const double someDepartureEpochSteps, const double aDepartureEpochStepSize, const double aTimeOfFlightMinimum, const double aTimeOfFlightMaximum, @@ -83,7 +83,7 @@ struct LambertScannerInput : catalogPath( aCatalogPath ), databasePath( aDatabasePath ), departureEpochInitial( aDepartureEpochInitial ), - departureEpochSteps( aDepartureEpochSteps ), + departureEpochSteps( someDepartureEpochSteps ), departureEpochStepSize( aDepartureEpochStepSize ), timeOfFlightMinimum( aTimeOfFlightMinimum ), timeOfFlightMaximum( aTimeOfFlightMaximum ), @@ -101,7 +101,7 @@ struct LambertScannerInput //! Path to SQLite database to store output. const std::string databasePath; - //! Departure epoch. + //! Initial departure epoch. const DateTime departureEpochInitial; //! Number of departure epoch steps. From e76ef32dfe6417e9af20cb3c87fbb336c8b7c0cf Mon Sep 17 00:00:00 2001 From: Enne Hekma Date: Wed, 9 Mar 2016 20:08:03 +0100 Subject: [PATCH 6/7] Add loop to scan map plot script to generate scanmap for each departure epoch --- python/plot_lambert_scan_maps.py | 135 ++++++++++++++++--------------- 1 file changed, 70 insertions(+), 65 deletions(-) diff --git a/python/plot_lambert_scan_maps.py b/python/plot_lambert_scan_maps.py index 9f46125..be87374 100644 --- a/python/plot_lambert_scan_maps.py +++ b/python/plot_lambert_scan_maps.py @@ -33,7 +33,8 @@ print "------------------------------------------------------------------" print " D2D " print " 0.0.2 " -print " Copyright (c) 2015, K. Kumar (me@kartikkumar.com) " +print " Copyright (c) 2015-2016, K. Kumar (me@kartikkumar.com) " +print " Copyright (c) 2016, E.J. Hekma (ennehekma@gmail.com) " print "------------------------------------------------------------------" print "" @@ -73,75 +74,79 @@ print "Error %s:" % e.args[0] sys.exit(1) -# Fetch scan data. -map_order = "departure_" + config['map_order'] -scan_data = pd.read_sql("SELECT departure_object_id, arrival_object_id, min(transfer_delta_v) \ - FROM lambert_scanner_results \ - GROUP BY departure_object_id, arrival_object_id;", \ - database) -scan_data.columns = ['departure_object_id','arrival_object_id','transfer_delta_v'] -scan_map = scan_data.pivot(index='departure_object_id', \ - columns='arrival_object_id', \ - values='transfer_delta_v') -scan_order = pd.read_sql("SELECT DISTINCT departure_object_id, " + map_order + " \ - FROM lambert_scanner_results \ - ORDER BY " + map_order + " ASC", \ + +departure_epochs = pd.read_sql("SELECT DISTINCT departure_epoch \ + FROM lambert_scanner_results;", \ database) -scan_map = scan_map.reindex(index=scan_order['departure_object_id'], \ - columns=scan_order['departure_object_id']) +for i in xrange(0,departure_epochs.size): + c = departure_epochs['departure_epoch'][i] + print "Plotting scanmap with departure epoch: ",c,"Julian Date" + + # Fetch scan data. + map_order = "departure_" + config['map_order'] + scan_data = pd.read_sql("SELECT departure_object_id, arrival_object_id, min(transfer_delta_v), "+ map_order + " \ + FROM lambert_scanner_results \ + WHERE departure_epoch BETWEEN " + str(c-0.00001) +" AND "+str(c+0.00001) +" \ + GROUP BY departure_object_id, arrival_object_id;", \ + database) + scan_data.columns = ['departure_object_id','arrival_object_id','transfer_delta_v',str(map_order)] + scan_order = scan_data.sort_values(str(map_order)).drop_duplicates('departure_object_id')[['departure_object_id',str(map_order)]] + scan_map = scan_data.pivot(index='departure_object_id', \ + columns='arrival_object_id', \ + values='transfer_delta_v') + scan_map = scan_map.reindex(index=scan_order['departure_object_id'], \ + columns=scan_order['departure_object_id']) + + # Set up color map. + bins = np.linspace(scan_data['transfer_delta_v'].min(), scan_data['transfer_delta_v'].max(), 10) + groups = scan_data['transfer_delta_v'].groupby(np.digitize(scan_data['transfer_delta_v'], bins)) + levels = groups.mean().values + cmap_lin = plt.get_cmap(config['colormap']) + cmap = nlcmap(cmap_lin, levels) + + + # Plot heat map. + ax1 = plt.subplot2grid((15,15), (2, 0),rowspan=13,colspan=14) + heatmap = ax1.pcolormesh(scan_map.values, cmap=cmap, \ + vmin=scan_data['transfer_delta_v'].min(), \ + vmax=scan_data['transfer_delta_v'].max()) + ax1.set_xticks(np.arange(scan_map.shape[1] + 1)+0.5) + ax1.set_xticklabels(scan_map.columns, rotation=90) + ax1.set_yticks([]) + ax1.tick_params(axis='both', which='major', labelsize=config['tick_label_size']) + ax1.set_xlim(0, scan_map.shape[1]) + ax1.set_ylim(0, scan_map.shape[0]) + ax1.set_xlabel('Departure object',fontsize=config['axis_label_size']) + ax1.set_ylabel('Arrival object',fontsize=config['axis_label_size']) + + # Plot axis ordering. + ax2 = plt.subplot2grid((15,15), (0, 0),rowspan=2,colspan=14,sharex=ax1) + ax2.step(np.arange(0.5,scan_map.shape[1]+.5),scan_order[str(map_order)],'k',linewidth=2.0) + ax2.get_yaxis().set_major_formatter(plt.FormatStrFormatter('%.2e')) + ax2.tick_params(axis='both', which='major', labelsize=config['tick_label_size']) + plt.setp(ax2.get_xticklabels(), visible=False) + ax2.set_ylabel(config['map_order_axis_label'],fontsize=config['axis_label_size']) + + # Plot color bar. + ax3 = plt.subplot2grid((15,15), (0, 14), rowspan=15) + color_bar = plt.colorbar(heatmap,cax=ax3) + color_bar.ax.get_yaxis().labelpad = 20 + color_bar.ax.set_ylabel(r'Total transfer $\Delta V$ [km/s]', rotation=270) + + plt.tight_layout() + + # Save figure to file. + plt.savefig(config["output_directory"] + "/" + config["scan_figure"] + "_"+str(i+1) + ".png", \ + dpi=config["figure_dpi"]) + plt.close() + # database.close() + print "Figure ",i+1," generated successfully...." +print "Figure generated successfully!" +print "" # Close SQLite database if it's still open. if database: database.close() - -print "Plotting heat map ..." - -# Plot heat map. - -# Set up color map. -bins = np.linspace(scan_data['transfer_delta_v'].min(), scan_data['transfer_delta_v'].max(), 10) -groups = scan_data['transfer_delta_v'].groupby(np.digitize(scan_data['transfer_delta_v'], bins)) -levels = groups.mean().values -cmap_lin = plt.get_cmap(config['colormap']) -cmap = nlcmap(cmap_lin, levels) - -# Plot heat map. -ax1 = plt.subplot2grid((15,15), (2, 0),rowspan=13,colspan=14) -heatmap = ax1.pcolormesh(scan_map.values, cmap=cmap, \ - vmin=scan_data['transfer_delta_v'].min(), \ - vmax=scan_data['transfer_delta_v'].max()) -ax1.set_xticks(np.arange(scan_map.shape[1] + 1)+0.5) -ax1.set_xticklabels(scan_map.columns, rotation=90) -ax1.set_yticks([]) -ax1.tick_params(axis='both', which='major', labelsize=config['tick_label_size']) -ax1.set_xlim(0, scan_map.shape[1]) -ax1.set_ylim(0, scan_map.shape[0]) -ax1.set_xlabel('Departure object',fontsize=config['axis_label_size']) -ax1.set_ylabel('Arrival object',fontsize=config['axis_label_size']) - -# Plot axis ordering. -ax2 = plt.subplot2grid((15,15), (0, 0),rowspan=2,colspan=14,sharex=ax1) -ax2.plot(scan_order[str(map_order)],'k',linewidth=2.0) -ax2.get_yaxis().set_major_formatter(plt.FormatStrFormatter('%.2e')) -ax2.tick_params(axis='both', which='major', labelsize=config['tick_label_size']) -plt.setp(ax2.get_xticklabels(), visible=False) -ax2.set_ylabel(config['map_order_axis_label'],fontsize=config['axis_label_size']) - -# Plot color bar. -ax3 = plt.subplot2grid((15,15), (0, 14), rowspan=15) -color_bar = plt.colorbar(heatmap,cax=ax3) -color_bar.ax.get_yaxis().labelpad = 20 -color_bar.ax.set_ylabel(r'Total transfer $\Delta V$ [km/s]', rotation=270) - -plt.tight_layout() - -# Save figure to file. -plt.savefig(config["output_directory"] + "/" + config["scan_figure"] + ".pdf", \ - dpi=config["figure_dpi"]) - -print "Figure generated successfully!" -print "" - # Stop timer end_time = time.time( ) From 0a9739105af45a3a946e2869e9ffe3fd366fbf64 Mon Sep 17 00:00:00 2001 From: Enne Hekma Date: Thu, 10 Mar 2016 12:26:56 +0100 Subject: [PATCH 7/7] Fix small formatting issues and update Panda version in requirements. --- config/lambert_scanner.json.empty | 2 +- python/plot_lambert_scan_maps.py | 52 ++++++++++++++++++------------- python/requirements.txt | 2 +- src/lambertScanner.cpp | 6 ++-- 4 files changed, 34 insertions(+), 28 deletions(-) diff --git a/config/lambert_scanner.json.empty b/config/lambert_scanner.json.empty index 224029f..a4c1542 100644 --- a/config/lambert_scanner.json.empty +++ b/config/lambert_scanner.json.empty @@ -20,7 +20,7 @@ // If departure_epoch array is left empty, the TLE epoch for the departure object is assumed. "departure_epoch" : [], - // Set departure epoch grid: [range (s), # of steps] + // Set departure epoch grid: [range (s), # of steps]. "departure_epoch_grid" : [,], // Set time-of-flight grid: [min (s), max (s), # of steps]. diff --git a/python/plot_lambert_scan_maps.py b/python/plot_lambert_scan_maps.py index be87374..5657627 100644 --- a/python/plot_lambert_scan_maps.py +++ b/python/plot_lambert_scan_maps.py @@ -75,40 +75,47 @@ sys.exit(1) -departure_epochs = pd.read_sql("SELECT DISTINCT departure_epoch \ - FROM lambert_scanner_results;", \ - database) +departure_epochs = pd.read_sql("SELECT DISTINCT departure_epoch \ + FROM lambert_scanner_results;", \ + database) for i in xrange(0,departure_epochs.size): c = departure_epochs['departure_epoch'][i] - print "Plotting scanmap with departure epoch: ",c,"Julian Date" + print "Plotting scan map with departure epoch: ",c,"Julian Date" # Fetch scan data. map_order = "departure_" + config['map_order'] - scan_data = pd.read_sql("SELECT departure_object_id, arrival_object_id, min(transfer_delta_v), "+ map_order + " \ - FROM lambert_scanner_results \ - WHERE departure_epoch BETWEEN " + str(c-0.00001) +" AND "+str(c+0.00001) +" \ - GROUP BY departure_object_id, arrival_object_id;", \ + scan_data = pd.read_sql("SELECT departure_object_id, arrival_object_id, \ + min(transfer_delta_v), "+ map_order + " \ + FROM lambert_scanner_results \ + WHERE departure_epoch BETWEEN " + str(c-0.00001) + " \ + AND "+str(c+0.00001) +" \ + GROUP BY departure_object_id, arrival_object_id;", \ database) - scan_data.columns = ['departure_object_id','arrival_object_id','transfer_delta_v',str(map_order)] - scan_order = scan_data.sort_values(str(map_order)).drop_duplicates('departure_object_id')[['departure_object_id',str(map_order)]] - scan_map = scan_data.pivot(index='departure_object_id', \ - columns='arrival_object_id', \ + scan_data.columns = ['departure_object_id','arrival_object_id', \ + 'transfer_delta_v',str(map_order)] + scan_order = scan_data.sort_values(str(map_order)) \ + .drop_duplicates('departure_object_id')[ \ + ['departure_object_id',str(map_order)]] + + scan_map = scan_data.pivot(index='departure_object_id', \ + columns='arrival_object_id', values='transfer_delta_v') - scan_map = scan_map.reindex(index=scan_order['departure_object_id'], \ + scan_map = scan_map.reindex(index=scan_order['departure_object_id'], \ columns=scan_order['departure_object_id']) # Set up color map. - bins = np.linspace(scan_data['transfer_delta_v'].min(), scan_data['transfer_delta_v'].max(), 10) - groups = scan_data['transfer_delta_v'].groupby(np.digitize(scan_data['transfer_delta_v'], bins)) + bins = np.linspace(scan_data['transfer_delta_v'].min(), \ + scan_data['transfer_delta_v'].max(), 10) + groups = scan_data['transfer_delta_v'].groupby( \ + np.digitize(scan_data['transfer_delta_v'], bins)) levels = groups.mean().values cmap_lin = plt.get_cmap(config['colormap']) cmap = nlcmap(cmap_lin, levels) - # Plot heat map. ax1 = plt.subplot2grid((15,15), (2, 0),rowspan=13,colspan=14) - heatmap = ax1.pcolormesh(scan_map.values, cmap=cmap, \ - vmin=scan_data['transfer_delta_v'].min(), \ + heatmap = ax1.pcolormesh(scan_map.values, cmap=cmap, \ + vmin=scan_data['transfer_delta_v'].min(), \ vmax=scan_data['transfer_delta_v'].max()) ax1.set_xticks(np.arange(scan_map.shape[1] + 1)+0.5) ax1.set_xticklabels(scan_map.columns, rotation=90) @@ -136,17 +143,18 @@ plt.tight_layout() # Save figure to file. - plt.savefig(config["output_directory"] + "/" + config["scan_figure"] + "_"+str(i+1) + ".png", \ - dpi=config["figure_dpi"]) + plt.savefig(config["output_directory"] + "/" + config["scan_figure"] + "_"+str(i+1) + \ + ".png", dpi=config["figure_dpi"]) plt.close() - # database.close() print "Figure ",i+1," generated successfully...." print "Figure generated successfully!" print "" + # Close SQLite database if it's still open. if database: database.close() + # Stop timer end_time = time.time( ) @@ -157,4 +165,4 @@ print "------------------------------------------------------------------" print " Exited successfully! " print "------------------------------------------------------------------" -print "" \ No newline at end of file +print "" diff --git a/python/requirements.txt b/python/requirements.txt index 9764c44..fc3d1f9 100644 --- a/python/requirements.txt +++ b/python/requirements.txt @@ -1,4 +1,4 @@ commentjson==0.4 matplotlib==1.4.2 -pandas==0.15.2 +pandas==0.17.1 scipy==0.14.1 diff --git a/src/lambertScanner.cpp b/src/lambertScanner.cpp index 0fd4158..a61f963 100644 --- a/src/lambertScanner.cpp +++ b/src/lambertScanner.cpp @@ -33,7 +33,6 @@ namespace d2d //! Execute lambert_scanner. void executeLambertScanner( const rapidjson::Document& config ) { - // Verify config parameters. Exception is thrown if any of the parameters are missing. const LambertScannerInput input = checkLambertScannerInput( config ); @@ -174,7 +173,6 @@ void executeLambertScanner( const rapidjson::Document& config ) departureEpoch = departureObject.Epoch( ); } - // Loop over arrival objects. for ( unsigned int j = 0; j < tleObjects.size( ); j++ ) { @@ -189,10 +187,10 @@ void executeLambertScanner( const rapidjson::Document& config ) const int arrivalObjectId = static_cast< int >( arrivalObject.NoradNumber( ) ); // Loop over departure epoch grid. - for (int l = 0; l < input.departureEpochSteps; ++l) + for ( int m = 0; m < input.departureEpochSteps; ++m ) { DateTime departureEpoch = input.departureEpochInitial; - departureEpoch = departureEpoch.AddSeconds( input.departureEpochStepSize * l ); + departureEpoch = departureEpoch.AddSeconds( input.departureEpochStepSize * m ); const Eci tleDepartureState = sgp4Departure.FindPosition( departureEpoch ); const Vector6 departureState = getStateVector( tleDepartureState );