diff --git a/config/lambert_scanner.json.empty b/config/lambert_scanner.json.empty index a4c1542..53c2daf 100644 --- a/config/lambert_scanner.json.empty +++ b/config/lambert_scanner.json.empty @@ -7,18 +7,24 @@ "mode" : "lambert_scanner", // Set path to TLE catalog file. - "catalog" : "../data/catalog/test_catalog.txt", + "catalog_path" : "../data/test_catalog.txt", // Set path to output database (SQLite). // WARNING: if the database file already exists, it will be overwritten! - "database" : "../data/test_lambert_scanner_catalog.db", + "database_path" : "../data/test_catalog_lambert_scanner.db", + + // Set length of sequence. + // The scanner will generate a sequence of grid searches to map the transfer costs for a + // succession of transfer legs. + // N.B.: The number of legs is: sequence_length - 1. + "sequence_length" : , // 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 left to right. // If departure_epoch array is left empty, the TLE epoch for the departure object is assumed. - "departure_epoch" : [], + "departure_epoch" : [,,], // Set departure epoch grid: [range (s), # of steps]. "departure_epoch_grid" : [,], @@ -26,6 +32,11 @@ // Set time-of-flight grid: [min (s), max (s), # of steps]. "time_of_flight_grid" : [,,], + // Set fixed stay time [s]. + // The stay time is added to the end of each leg, i.e., the departure epoch for the next leg is + // delayed by stay_time. + "stay_time" : , + // Set flag indicating if transfers are prograde. If set to false, retrograde transfers are // computed. // N.B.: Since SQLite doesn't have a boolean type, the flag is stored as an integer @@ -35,8 +46,7 @@ // Set maximum number of transfer revolutions (N). "revolutions_maximum" : , - // Set number of transfers to include in shortlist and absolute path to output file [N, file]. - // The shortlist is based on the N transfers specified with the lowest transfers Delta-V. - // If N is set to 0 no output will be written to file. - "shortlist" : [0,""] + // Set path to sequences file with summary of best multi-leg transfers for each sequence. + // The sequences are sorted from lowest to highest Delta-V. + "sequences_path" : "" } \ No newline at end of file diff --git a/include/D2D/lambertScanner.hpp b/include/D2D/lambertScanner.hpp index 6f1218e..4e5dedd 100644 --- a/include/D2D/lambertScanner.hpp +++ b/include/D2D/lambertScanner.hpp @@ -12,14 +12,26 @@ #include #include +#include #include #include +#include "D2D/tools.hpp" + namespace d2d { +//! Forward declaration. +struct LambertScannerInput; +struct LambertPorkChopPlotGridPoint; + +//! Lambert pork-chop plot, consisting of list of grid points. +typedef std::vector< LambertPorkChopPlotGridPoint > LambertPorkChopPlot; +//! Collection of all Lambert pork-chop plots with leg, departure object & arrival object IDs. +typedef std::map< PorkChopPlotId, LambertPorkChopPlot > AllLambertPorkChopPlots; + //! Execute lambert_scanner. /*! * Executes lambert_scanner application mode that performs a grid search to compute \f$\Delta V\f$ @@ -35,6 +47,141 @@ namespace d2d */ void executeLambertScanner( const rapidjson::Document& config ); +//! Check lambert_scanner input parameters. +/*! + * Checks that all inputs for the lambert_scanner application mode are valid. If not, an error is + * thrown with a short description of the problem. If all inputs are valid, a data struct + * containing all the inputs is returned, which is subsequently used to execute lambert_scanner + * and related functions. + * + * @sa executeLambertScanner, LambertScannerInput + * @param[in] config User-defined configuration options (extracted from JSON input file) + * @return Struct containing all valid input to execute lambert_scanner + */ +LambertScannerInput checkLambertScannerInput( const rapidjson::Document& config ); + +//! Create lambert_scanner table. +/*! + * Creates lambert_scanner and lambert_sequences tables in SQLite database used to store results + * obtained from running the lambert_scanner application mode. + * + * @sa executeLambertScanner + * @param[in] database SQLite database handle + * @param[in] sequenceLength Number of targets in sequence + */ +void createLambertScannerTables( SQLite::Database& database, const int sequenceLength ); + +//! Write best multi-leg Lambert transfers for each sequence to file. +/*! + * Writes list of best multi-leg Lambert transfers for each sequence to file. The list is based on + * retrieving the multi-leg transfers with the lowest \f$\Delta V\f$ for each sequence from the + * "sequences" tables in the database. The list is sorted from lowest to highest \f$\Delta V\f$. + * + * @sa executeLambertScanner, createLambertScannerTables + * @param[in] database SQLite database handle + * @param[in] sequencesPath Path to sequences file + * @param[in] sequenceLength Length of each sequence + */ +void writeSequencesToFile( SQLite::Database& database, + const std::string& sequencesPath, + const int sequenceLength ); + +//! Recurse through sequences leg-by-leg and compute pork-chop plots. +/*! + * Recurses through pool of TLE objects to compute sequences and generates pork-chop plots based + * on specified departure and arrival objects for a given leg, whilst traversing the tree of + * sequences. The departure and arrival epochs have to be generated apriori for each leg and + * are provided by the user. The transfers are computed using the Lambert solver (Izzo, 2014). + * + * All of the pork-chop plots are stored in the container provided by the user (AllPorkChopPlots). + * + * @sa executeLambertScanner + * @param[in] currentSequencePosition Current position along sequence + * @param[in] tleObjects Pool of TLE objects to select from + * @param[in] allEpochs Pre-computed departure and arrival epochs for each + * leg + * @param[in] isPrograde Flag indicating if prograde transfer should be + * computed (false = retrograde) + * @param[in] revolutionsMaximum Maximum number of revolutions + * @param[out] sequence Sequence of TLE objects + * @param[in] transferId Counter for unique ID assigned to transfers + * @param[out] allPorkChopPlots Set of all Lambert pork-chop plots + */ +void recurseLambertTransfers( const int currentSequencePosition, + const TleObjects& tleObjects, + const AllEpochs& allEpochs, + const bool isPrograde, + const int revolutionsMaximum, + Sequence& sequence, + int& transferId, + AllLambertPorkChopPlots& allPorkChopPlots ); + +//! Compute Lambert pork-chop plot. +/*! + * Computes a Lambert pork-chop plot, given TLE departure and arrival objects and a list of + * departure and arrival epoch pairs. The transfers are computed using the Lambert solver + * implemented in PyKep (Izzo, 2014). The sense (prograde/retrograde) of the transfer and the + * maximum number of revolutions permitted are specified by the user in the configuration file. + * The list of departure and arrival epoch pairs are pre-computed. The 2D pork-chop plot grid is + * flattened into a 1D list with the grid point coordinates specified as a departure-arrival epoch + * pair. + * + * The result is stored as a list of LambertPorkChopPlotGridPoint points, which contain the + * departure epoch, time-of-flight (departure-arrival epoch) and Lambert transfer data. + * + * @sa LambertPorkChopPlotGridPoint, computeAllPorkChopPlotEpochs + * @param[in] departureObject Departure TLE object + * @param[in] arrivalObject Arrival TLE object + * @param[in] listOfEpochs Flattened list of departure-arrival epoch pairs + * @param[in] isPrograde Flag indicating if prograde Lambert transfer should be computed + * (false = retrograde) + * @param[in] revolutionsMaximum Maximum number of revolutions + * @param[in] transferId Counter for unique ID assigned to transfers + * @return List of LambertPorkChopPlotGridPoint objects + */ +LambertPorkChopPlot computeLambertPorkChopPlot( const Tle& departureObject, + const Tle& arrivalObject, + const ListOfEpochs& listOfEpochs, + const bool isPrograde, + const int revolutionsMaximum, + int& transferId ); + +//! Recurse through sequence leg-by-leg and compute multi-leg transfers. +/*! + * Recurses leg-by-leg through a given sequence and computes multi-leg transfers, taking into + * account matching between the arrival epoch from the previous leg, stay time at target and + * departure epoch for next leg. All possible multi-leg transfers are computed by traversing the + * tree of pork-chop plots for each leg. + * + * All of the multi-leg transfers for the specified transfer are stored in a list. + * + * The multiLegTransferData, launchEpoch and lastArrivalEpoch input parameters are storage variables + * that are used during the recursion to keep track of values whilst traversing the sequence tree. + * The multiLegTransferData parameter should be passed to the function without any values. The + * launchEpoch and lastArrivalEpoch parameters are set to default values and should not be passed + * by the user. + * + * @sa executeLambertScanner, ListOfSequences, AllPorkChopPlots, ListOfMultiLegTransfers + * @param[in] currentSequencePosition Current position along sequence + * @param[in] sequence Sequence of TLE objects + * @param[in] allPorkChopPlots Collection of all pork-chop plots for all sequences + * @param[in] stayTime Fixed stay time at each target + * @param[out] listOfMultiLegTransfers List of all multi-leg transfers for given sequence + * @param[out] multiLegTransferData Storage container for multi-leg transfer data (for + * internal use) + * @param[in] launchEpoch Storage variable for launch epoch (default: DateTime( )) + * @param[in] lastArrivalEpoch Storage variable for arrival epoch from previous leg + * (default: DateTime( )) + */ +void recurseMuiltiLegLambertTransfers( const int currentSequencePosition, + const Sequence& sequence, + AllLambertPorkChopPlots& allPorkChopPlots, + const double stayTime, + ListOfMultiLegTransfers& listOfMultiLegTransfers, + MultiLegTransferData& multiLegTransferData, + DateTime launchEpoch = DateTime( ), + DateTime lastArrivalEpoch = DateTime( ) ); + //! Input for lambert_scanner application mode. /*! * Data struct containing all valid lambert_scanner input parameters. This struct is populated by @@ -54,6 +201,7 @@ struct LambertScannerInput * @sa checkLambertScannerInput, executeLambertScanner * @param[in] aCatalogPath Path to TLE catalog * @param[in] aDatabasePath Path to SQLite database + * @param[in] aSequenceLength Sequence length * @param[in] aDepartureEpochInitial Departure epoch grid initial epoch * @param[in] someDepartureEpochSteps Number of steps to take in departure epoch grid * @param[in] aDepartureEpochStepSize Departure epoch grid step size (derived parameter) [s] @@ -61,14 +209,15 @@ struct LambertScannerInput * @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] aStayTime Fixed stay time at arrival object [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] aSequencesPath Path to sequences file */ LambertScannerInput( const std::string& aCatalogPath, const std::string& aDatabasePath, + const int aSequenceLength, const DateTime& aDepartureEpochInitial, const double someDepartureEpochSteps, const double aDepartureEpochStepSize, @@ -76,12 +225,13 @@ struct LambertScannerInput const double aTimeOfFlightMaximum, const double someTimeOfFlightSteps, const double aTimeOfFlightStepSize, + const double aStayTime, const bool progradeFlag, const int aRevolutionsMaximum, - const int aShortlistLength, - const std::string& aShortlistPath ) + const std::string& aSequencesPath ) : catalogPath( aCatalogPath ), databasePath( aDatabasePath ), + sequenceLength( aSequenceLength ), departureEpochInitial( aDepartureEpochInitial ), departureEpochSteps( someDepartureEpochSteps ), departureEpochStepSize( aDepartureEpochStepSize ), @@ -89,10 +239,10 @@ struct LambertScannerInput timeOfFlightMaximum( aTimeOfFlightMaximum ), timeOfFlightSteps( someTimeOfFlightSteps ), timeOfFlightStepSize( aTimeOfFlightStepSize ), + stayTime( aStayTime ), isPrograde( progradeFlag ), revolutionsMaximum( aRevolutionsMaximum ), - shortlistLength( aShortlistLength ), - shortlistPath( aShortlistPath ) + sequencesPath( aSequencesPath ) { } //! Path to TLE catalog. @@ -101,6 +251,9 @@ struct LambertScannerInput //! Path to SQLite database to store output. const std::string databasePath; + //! Length of sequence (numberOfLegs = sequenceLength - 1). + const int sequenceLength; + //! Initial departure epoch. const DateTime departureEpochInitial; @@ -122,61 +275,149 @@ struct LambertScannerInput //! Time-of-flight step size [s]. const double timeOfFlightStepSize; + //! Fixed stay time at arrival object, i.e., departure for next leg is delayed by stay time. + const double stayTime; + //! Flag indicating if transfers are prograde. False indicates retrograde. const bool isPrograde; //! Maximum number of revolutions (N) for transfer. Number of revolutions is 2*N+1. const int revolutionsMaximum; - //! Number of entries (lowest transfer \f$\Delta V\f$) to include in shortlist. - const int shortlistLength; - - //! Path to shortlist file. - const std::string shortlistPath; + //! Path to sequences file. + const std::string sequencesPath; protected: private: }; -//! Check lambert_scanner input parameters. +//! Grid point in pork-chop plot. /*! - * Checks that all inputs for the lambert_scanner application mode are valid. If not, an error is - * thrown with a short description of the problem. If all inputs are valid, a data struct - * containing all the inputs is returned, which is subsequently used to execute lambert_scanner - * and related functions. + * Data struct containing all data pertaining to a grid point in a pork-chop plot, i.e., departure + * epoch, time-of-flight and all data related to the computed transfer. * - * @sa executeLambertScanner, LambertScannerInput - * @param[in] config User-defined configuration options (extracted from JSON input file) - * @return Struct containing all valid input to execute lambert_scanner + * @sa recurseTransfers */ -LambertScannerInput checkLambertScannerInput( const rapidjson::Document& config ); +struct LambertPorkChopPlotGridPoint +{ +public: -//! Create lambert_scanner table. -/*! - * Creates lambert_scanner table in SQLite database used to store results obtaned from running - * the lambert_scanner application mode. - * - * @sa executeLambertScanner - * @param[in] database SQLite database handle - */ -void createLambertScannerTable( SQLite::Database& database ); + //! Construct data struct. + /*! + * Constructs data struct based on departure epoch, time-of-flight and transfer data for grid + * point in pork-chop plot. + * + * @param[in] aTransferId A unique transfer ID + * @param[in] aDepartureEpoch A departure epoch corresponding to a grid point + * @param[in] anArrivalEpoch An arrival epoch corresponding to a grid point + * @param[in] aTimeOfFlight A time-of-flight (arrival-departure epoch) for a grid point + * @param[in] someRevolutions Number of revolutions for transfer + * @param[in] progradeFlag Flag indicating if transfer is prograde (false = retrograde) + * @param[in] aDepartureState Departure state corresponding to a grid point + * @param[in] aDepartureStateKepler Departure state in Keplerian elements corresponding to a + * grid point + * @param[in] anArrivalState Arrival state corresponding to a grid point + * @param[in] anArrivalStateKepler Arrival state in Keplerian elements corresponding to a grid + * point + * @param[in] aTransferStateKepler Transfer state in Keplerian elements corresponding to a grid + * point + * @param[in] aDepartureDeltaV Computed departure \f$\Delta V\f$ + * @param[in] anArrivalDeltaV Computed arrival \f$\Delta V\f$ + * @param[in] aTransferDeltaV Total computed transfer \f$\Delta V\f$ + */ + LambertPorkChopPlotGridPoint( const int aTransferId, + const DateTime& aDepartureEpoch, + const DateTime& anArrivalEpoch, + const double aTimeOfFlight, + const int someRevolutions, + const bool progradeFlag, + const Vector6& aDepartureState, + const Vector6& aDepartureStateKepler, + const Vector6& anArrivalState, + const Vector6& anArrivalStateKepler, + const Vector6& aTransferStateKepler, + const Vector3& aDepartureDeltaV, + const Vector3& anArrivalDeltaV, + const double aTransferDeltaV ) + : transferId( aTransferId ), + departureEpoch( aDepartureEpoch ), + arrivalEpoch( anArrivalEpoch ), + timeOfFlight( aTimeOfFlight ), + revolutions( someRevolutions ), + isPrograde( progradeFlag ), + departureState( aDepartureState ), + departureStateKepler( aDepartureStateKepler ), + arrivalState( anArrivalState ), + arrivalStateKepler( anArrivalStateKepler ), + transferStateKepler( aTransferStateKepler ), + departureDeltaV( aDepartureDeltaV ), + arrivalDeltaV( anArrivalDeltaV ), + transferDeltaV( aTransferDeltaV ) + { } -//! Write transfer shortlist to file. -/*! - * Writes shortlist of debris-to-debris Lambert transfers to file. The shortlist is based on the - * requested number of transfers with the lowest transfer \f$\Delta V\f$, retrieved by sorting the - * transfers in the SQLite database. - * - * @sa executeLambertScanner, createLambertScannerTable - * @param[in] database SQLite database handle - * @param[in] shortlistNumber Number of entries to include in shortlist (if it exceeds number of - * entries in database table, the whole table is written to file) - * @param[in] shortlistPath Path to shortlist file - */ -void writeTransferShortlist( SQLite::Database& database, - const int shortlistNumber, - const std::string& shortlistPath ); + //! Overload operator-=. + /*! + * Overloads operator-= to assign current object to object provided as input. + * + * WARNING: this is a dummy overload to get by the problem of adding a + * LambertPorkChopPlotGridPoint object to a STL container! It does not correctly assign + * the current object to the dummy grid point provided! + * + * @sa executeLambertScanner + * @param[in] dummyGridPoint Dummy grid point that is ignored + * @return The current object + */ + LambertPorkChopPlotGridPoint& operator=( const LambertPorkChopPlotGridPoint& dummyGridPoint ) + { + return *this; + } + + //! Unique transfer ID. + const int transferId; + + //! Departure epoch. + const DateTime departureEpoch; + + //! Arrival epoch. + const DateTime arrivalEpoch; + + //! Time of flight [s]. + const double timeOfFlight; + + //! Number of revolutions (N) for transfer. + const int revolutions; + + //! Flag indicating if transfers are prograde. False indicates retrograde. + const bool isPrograde; + + //! Departure state [km; km/s]. + const Vector6 departureState; + + //! Departure state in Keplerian elements [km; -; rad; rad; rad; rad]. + const Vector6 departureStateKepler; + + //! Arrival state [km; km/s]. + const Vector6 arrivalState; + + //! Arrival state in Keplerian elements [km; -; rad; rad; rad; rad]. + const Vector6 arrivalStateKepler; + + //! Transfer state in Keplerian elements [km; -; rad; rad; rad; rad]. + const Vector6 transferStateKepler; + + //! Departure \f$\Delta V\f$ [km/s]. + const Vector3 departureDeltaV; + + //! Arrival \f$\Delta V\f$ [km/s]. + const Vector3 arrivalDeltaV; + + //! Total transfer \f$\Delta V\f$ [km/s]. + const double transferDeltaV; + +protected: +private: +}; } // namespace d2d diff --git a/include/D2D/tools.hpp b/include/D2D/tools.hpp index 035c7b0..8bd5d2d 100644 --- a/include/D2D/tools.hpp +++ b/include/D2D/tools.hpp @@ -13,11 +13,14 @@ #include #include #include +#include +#include #include #include +#include #include #include @@ -30,6 +33,31 @@ namespace d2d { +//! Forward declarations +struct MultiLegTransfer; +struct TransferData; + +//! List of TLE objects generated from TLE strings. +typedef std::vector< Tle > TleObjects; +//! Sequence of TLE Objects. +typedef std::vector< Tle > Sequence; +//! List of TLE object ID sequences (key=sequence ID). +typedef std::map< int, Sequence > ListOfSequences; + +//! Departure-Arrival epoch pair. +typedef std::pair< DateTime, DateTime > Epochs; +//! List of departure-arrival epoch pairs. +typedef std::vector< Epochs > ListOfEpochs; +//! Collection of lists of departure-arrival epoch pairs (key=leg ID). +typedef std::map< int, ListOfEpochs > AllEpochs; + +//! List of multi-leg transfer data. +typedef std::vector< TransferData > MultiLegTransferData; +//! List of multi-leg transfers. +typedef std::vector< MultiLegTransfer > ListOfMultiLegTransfers; +//! Collection of all multi-leg transfers for all sequences (key=sequence ID). +typedef std::map< int, ListOfMultiLegTransfers > AllMultiLegTransfers; + //! Create custom CATCH Approx object with tolerance for comparing doubles. /*! * Creates a custom CATCH Approx object that can be used for comparing doubles in unit tests. The @@ -84,7 +112,7 @@ StateHistory sampleKeplerOrbit( const Vector6& initialState, * Samples a SGP4 orbit and generates a state-history stored in a STL map (key=epoch). The * SGP4 orbit is sampled by using FindPosition() provided with libsgp4. * - * @param[in] tle Two-line element data of the object to be propagated + * @param[in] tle Two-line element data of the object to be propagated * @param[in] initialEpochJulian Starting epoch for the SGP4 propagator * (default = 0.0) [Julian date] * @param[in] propagationTime Total propagation time [s] @@ -102,17 +130,17 @@ StateHistory sampleSGP4Orbit( const Tle& tle, * convertCartesianStateToTwoLineElements function in the Atom library. * * @sa atom::convertCartesianStateToTwoLineElements - * @param[in] propagatedCartesianState The state obtained after propagating virtual TLE + * @param[in] propagatedCartesianState The state obtained after propagating virtual TLE * using SGP4 with time-of-flight = 0.0 - * @param[in] trueCartesianState The true cartesian state corresponding to zero + * @param[in] trueCartesianState The true cartesian state corresponding to zero * time-of-flight - * @param[in] relativeTolerance Relative difference between the propagated and true - * Cartesian state is checked against the relative + * @param[in] relativeTolerance Relative difference between the propagated and true + * Cartesian state is checked against the relative * tolerance - * @param[in] absoluteTolerance Absolute difference between the propagated and true - * Cartesian state is checked against the absolute + * @param[in] absoluteTolerance Absolute difference between the propagated and true + * Cartesian state is checked against the absolute * tolerance - * @return Returns boolean 'true' if the test passed, + * @return Returns boolean 'true' if the test passed, * 'false' otherwise */ bool executeVirtualTleConvergenceTest( const Vector6& propagatedCartesianState, @@ -228,6 +256,260 @@ void removeNewline( std::string& string ); */ int getTleCatalogType( const std::string& catalogFirstLine ); +//! Recurse leg-by-leg to generate list of TLE sequences. +/*! + * Recurses through pool of TLE objects to generate list of sequences containing TLE IDs (NORAD + * number). The sequences are all of a specified length. The final list of sequences generated + * contains a full enumeration of all possible sequences using the TLE object pool. + * + * @sa executeLambertScanner + * @param[in] currentSequencePosition Current position in sequence + * @param[in] tleObjects Pool of TLE objects to select from + * @param[in] sequence Sequence of TLE object IDs + * @param[in] sequenceId Counter for unique ID assigned to sequences + * @param[out] listOfSequences List of sequences generated + */ +void recurseSequences( const int currentSequencePosition, + const TleObjects& tleObjects, + Sequence& sequence, + int& sequenceId, + ListOfSequences& listOfSequences ); + +//! Compute departure-arrival epoch pairs for all pork-chop plots. +/*! + * Computes all departure and arrival epochs for each pork-chop plot at each leg. Since multiple + * combinations of departure epoch and time-of-flight can lead the same arrival epoch, the unique + * list of arrival epochs are extracted from each leg. A fixed stay time is added to each of these + * arrival epochs, yields the list of departure epochs for the subsequence leg. The function returns + * a map containing the departure-arrival epoch pairs belonging to the pork-chop plots for each leg. + * + * @param[in] sequenceLength Fixed length of multi-leg sequence + * @param[in] stayTime Fixed stay time at arrival object [s] + * @param[in] departureEpochInitial Initial departure epoch at start of grid + * @param[in] departureEpochSteps Number of steps within departure epoch grid + * @param[in] departureEpochStepSize Step size within departure epoch grid [s] + * @param[in] timeOfFlightMinimum Minimum time-of-flight [s] + * @param[in] timeOfFlightSteps Number of steps within time-of-flight grid + * @param[in] timeOfFlightStepSize Step size within time-of-flight grid [s] + */ +AllEpochs computeAllPorkChopPlotEpochs( const int sequenceLength, + const double stayTime, + const DateTime& departureEpochInitial, + const int departureEpochSteps, + const double departureEpochStepSize, + const double timeOfFlightMinimum, + const int timeOfFlightSteps, + const double timeOfFlightStepSize ); + +//! Pork-chop plot leg, departure object and arrival object IDs. +/*! + * Data struct containing leg, departure object and arrival object IDs that define the location of a + * given pork-chop plot along the sequence tree. These ID numbers serve to uniquely identify a + * pork-chop plot. The data struct is used as the key in a map to associate pork-chop plots along + * the sequence tree. This ensures that the pork-chop plots are only computed once. + * + * @sa recurseSequences, recurseTransfers, PorkChopPlotGridPoint + */ +struct PorkChopPlotId +{ +public: + + //! Construct data struct. + /*! + * Constructs data struct based on leg, departure object and arrival object IDs that define the + * location of a pork-chop plot within the sequence tree. + */ + PorkChopPlotId( const int aLegID, + const int aDepartureObjectId, + const int anArrivalObjectId ) + : legId( aLegID ), + departureObjectId( aDepartureObjectId ), + arrivalObjectId( anArrivalObjectId ) + { } + + //! Leg ID. + const int legId; + + //! Departure object ID. + const int departureObjectId; + + //! Arrival object ID + const int arrivalObjectId; + +protected: +private: +}; + +//! Overload ==-operator to compare PorkChopPlotId objects. +/*! + * Overloads ==-operator to compare two PorkChopPlotId objects. The comparison is based on + * sequentially checking that the leg, departure object and arrival object IDs are equal in both + * objects. + * + * @sa PorkChopPlotId + * @param[in] id1 First PorkChopPlotId object + * @param[in] id2 Second PorkChopPlotId object + * @return True if PorkChopPlotId objects are equal, false otherwise + */ +bool operator==( const PorkChopPlotId& id1, const PorkChopPlotId& id2 ); + +//! Overload !=-operator to compare PorkChopPlotId objects. +/*! + * Overloads !=-operator to compare two PorkChopPlotId objects. The comparison is based on + * negating the result obtained from calling the ==-operator on both objects. + * + * @sa PorkChopPlotId, + * @param[in] id1 First PorkChopPlotId object + * @param[in] id2 Second PorkChopPlotId object + * @return False if PorkChopPlotId objects are equal, true otherwise + */ +bool operator!=( const PorkChopPlotId& id1, const PorkChopPlotId& id2 ); + +//! Overload <-operator to compare PorkChopPlotId objects. +/*! + * Overloads <-operator to compare two PorkChopPlotId objects. The comparison is based on + * sequentially checking that the leg, departure object and arrival object IDs of the first + * PorkChopPlotId object (first argument) are less than the corresponding values stored in the + * second PorkChopPlotId object (second argument). + * + * @sa PorkChopPlotId + * @param[in] id1 First PorkChopPlotId object + * @param[in] id2 Second PorkChopPlotId object + * @return True if id1 is less than id2, false otherwise + */ +bool operator<( const PorkChopPlotId& id1, const PorkChopPlotId& id2 ); + +//! Overload >=-operator to compare PorkChopPlotId objects. +/*! + * Overloads <=-operator to compare two PorkChopPlotId objects. The comparison is based on + * checking that the results obtained from calling the ==-operator or the <-operator are true. + * + * @sa PorkChopPlotId + * @param[in] id1 First PorkChopPlotId object + * @param[in] id2 Second PorkChopPlotId object + * @return True if id1 is less than or equal to id2, false otherwise + */ +bool operator<=( const PorkChopPlotId& id1, const PorkChopPlotId& id2 ); + +//! Overload >-operator to compare PorkChopPlotId objects. +/*! + * Overloads >-operator to compare two PorkChopPlotId objects. The comparison is based on + * checking negating the results obtained from calling the <=-operator on both objects. + * + * @sa PorkChopPlotId + * @param[in] id1 First PorkChopPlotId object + * @param[in] id2 Second PorkChopPlotId object + * @return True if id1 is greater than id2, false otherwise + */ +bool operator>( const PorkChopPlotId& id1, const PorkChopPlotId& id2 ); + +//! Overload >=-operator to compare PorkChopPlotId objects. +/*! + * Overloads >=-operator to compare two PorkChopPlotId objects. The comparison is based on + * checking that the results obtained from calling the ==-operator or the >-operator are true. + * + * @sa PorkChopPlotId + * @param[in] id1 First PorkChopPlotId object + * @param[in] id2 Second PorkChopPlotId object + * @return True if id1 is greater than or equal to id2, false otherwise + */ +bool operator>=( const PorkChopPlotId& id1, const PorkChopPlotId& id2 ); + + +//! Data for single transfer. +/*! + * Data struct based on a single transfer leg. The transfer ID is set using the value stored in a + * pork-chop plot. The time-of-flight and transfer \f$\Delta V\f$ are extracting from the pork-chop + * plot. + */ +struct TransferData +{ +public: + + //! Construct data struct. + /*! + * Constructs data based on transfer ID for specified transfer in pork-chop plot and + * corresponding time-of-flight and \f$\Delta V\f$. + * + * @sa MultiLegTransfer + * @param[in] aTransferId A transfer ID extracted from a pork-chop plot + * @param[in] aTimeOfFlight Time-of-flight associated with a given transfer ID + * @param[in] aTransferDeltaV \f$\Delta V\f$ associated with a given transfer ID + */ + TransferData( const int aTransferId, + const double aTimeOfFlight, + const double aTransferDeltaV ) + : transferId( aTransferId ), + timeOfFlight( aTimeOfFlight ), + transferDeltaV( aTransferDeltaV ) + { } + + //! Overload operator-=. + /*! + * Overloads operator-= to assign current object to object provided as input. + * + * WARNING: this is a dummy overload to get by the problem of adding a TransferData object to a + * STL container! It does not correctly assign the current object to the dummy transfer + * data object provided! + * + * @sa MultiLegTransferData, MultiLegTransfer + * @param[in] dummyTransferData Dummy transfer data object that is ignored + * @return The current object + */ + TransferData& operator=( const TransferData& dummyTransferData ) + { + return *this; + } + + //! Transfer ID. + const int transferId; + + //! Time-of-flight [s]. + const double timeOfFlight; + + //! Transfer \f$\Delta V\f$ [km/s] + const double transferDeltaV; + +protected: +private: +}; + +//! Data for multi-leg transfer. +/*! + * Data struct based on a multi-leg transfer. The launch epoch is set as the departure epoch for the + * first leg. The list of transfer data is based on the transfer IDs for the individual legs and the + * associated time-of-flight and transfer \f$\Delta V\f$ values. + */ +struct MultiLegTransfer +{ +public: + + //! Construct data struct. + /*! + * Constructs data struct based on launch epoch from departure epoch window for the first leg, + * and a list of transfer data based on the transfer IDs for individual legs and associated + * time-of-flight and \f$\Delta V\f$ values. + * + * @param[in] aLaunchEpoch A launch epoch for a multi-leg transfer + * @param[in] someMultiLegTransferData A list of multiple transfers IDs and associated + * time-of-flight and \f$\Delta V\f$ values + */ + MultiLegTransfer( const DateTime& aLaunchEpoch, + const MultiLegTransferData& someMultiLegTransferData ) + : launchEpoch( aLaunchEpoch ), + multiLegTransferData( someMultiLegTransferData ) + { } + + //! Launch epoch (departure epoch for first leg). + DateTime launchEpoch; + + //! List of data per leg. + MultiLegTransferData multiLegTransferData; + +protected: +private: +}; + } // namespace d2d #endif // D2D_TOOLS_HPP diff --git a/src/atomScanner.cpp b/src/atomScanner.cpp index 5984f15..b10156b 100644 --- a/src/atomScanner.cpp +++ b/src/atomScanner.cpp @@ -23,6 +23,8 @@ #include #include +#include + #include #include diff --git a/src/j2Analysis.cpp b/src/j2Analysis.cpp index 8ed768c..795748c 100644 --- a/src/j2Analysis.cpp +++ b/src/j2Analysis.cpp @@ -17,10 +17,12 @@ #include #include -#include - #include +#include + +#include + #include "D2D/j2Analysis.hpp" #include "D2D/tools.hpp" #include "D2D/typedefs.hpp" diff --git a/src/lambertFetch.cpp b/src/lambertFetch.cpp index d7e153c..6066b4b 100644 --- a/src/lambertFetch.cpp +++ b/src/lambertFetch.cpp @@ -12,6 +12,8 @@ #include +#include + #include #include "D2D/lambertFetch.hpp" diff --git a/src/lambertScanner.cpp b/src/lambertScanner.cpp index 073e8e0..b8bc6dd 100644 --- a/src/lambertScanner.cpp +++ b/src/lambertScanner.cpp @@ -14,18 +14,16 @@ #include #include -#include - #include #include #include -#include + +#include #include #include #include "D2D/lambertScanner.hpp" -#include "D2D/tools.hpp" namespace d2d { @@ -60,10 +58,11 @@ void executeLambertScanner( const rapidjson::Document& config ) // Reset file stream to start of file. catalogFile.seekg( 0, std::ios::beg ); - typedef std::vector< std::string > TleStrings; - typedef std::vector< Tle > TleObjects; TleObjects tleObjects; + //! List of TLE strings extracted from catalog. + typedef std::vector< std::string > TleStrings; + while ( std::getline( catalogFile, catalogLine ) ) { TleStrings tleStrings; @@ -94,19 +93,124 @@ void executeLambertScanner( const rapidjson::Document& config ) SQLite::Database database( input.databasePath.c_str( ), SQLITE_OPEN_READWRITE|SQLITE_OPEN_CREATE ); - // Create table for Lambert scanner results in SQLite database. - std::cout << "Creating SQLite database table if needed ... " << std::endl; - createLambertScannerTable( database ); + // Create tables for Lambert scanner results in SQLite database. + std::cout << "Creating SQLite database tables if needed ... " << std::endl; + createLambertScannerTables( database, input.sequenceLength ); std::cout << "SQLite database set up successfully!" << std::endl; - // Start SQL transaction. - SQLite::Transaction transaction( database ); + std::cout << "Enumerating sequences and populating database ... " << std::endl; - // Setup insert query. - std::ostringstream lambertScannerTableInsert; - lambertScannerTableInsert - << "INSERT INTO lambert_scanner_results VALUES (" - << "NULL," + // Start SQL transaction to populate table of sequences. + SQLite::Transaction sequencesTransaction( database ); + + // Set up insert query to add sequences to table in database. + std::ostringstream lambertScannerSequencesTableInsert; + lambertScannerSequencesTableInsert + << "INSERT INTO sequences VALUES (" + << ":sequence_id,"; + for ( int i = 0; i < input.sequenceLength - 1; ++i ) + { + lambertScannerSequencesTableInsert + << ":target_" << i << ","; + } + lambertScannerSequencesTableInsert + << ":target_" << input.sequenceLength - 1 << ","; + for ( int i = 0; i < input.sequenceLength - 2; ++i ) + { + lambertScannerSequencesTableInsert + << "0,"; + } + lambertScannerSequencesTableInsert + << "0," + << "0.0," + << "0.0," + << "0.0," + << "0.0" + << ");"; + + SQLite::Statement sequencesQuery( database, lambertScannerSequencesTableInsert.str( ) ); + + // Enumerate all sequences for the given pool of TLE objects using a recursive function that + // generates IDs leg-by-leg. + ListOfSequences listOfSequences; + Sequence sequence( input.sequenceLength ); + int currentSequencePosition = 0; + int sequenceId = 1; + + recurseSequences( currentSequencePosition, tleObjects, sequence, sequenceId, listOfSequences ); + + // Write list of sequences to table in database. + for ( ListOfSequences::iterator iteratorSequences = listOfSequences.begin( ); + iteratorSequences != listOfSequences.end( ); + ++iteratorSequences ) + { + sequencesQuery.bind( ":sequence_id", iteratorSequences->first ); + + for ( unsigned int j = 0; j < sequence.size( ); ++j ) + { + std::ostringstream sequencesQueryTargetNumber; + sequencesQueryTargetNumber << ":target_" << j; + sequencesQuery.bind( + sequencesQueryTargetNumber.str( ), + static_cast< int >( iteratorSequences->second[ j ].NoradNumber( ) ) ); + } + + // Execute insert query. + sequencesQuery.executeStep( ); + + // Reset SQL insert query. + sequencesQuery.reset( ); + } + + // Commit sequences transaction. + sequencesTransaction.commit( ); + + std::cout << "Sequences successfully enumerated!" << std::endl; + + std::cout << "Pre-computing all departure-arrival epochs for each leg ... " << std::endl; + + // Pre-compute all departure and arrival epochs for each leg. + const AllEpochs allEpochs = computeAllPorkChopPlotEpochs( input.sequenceLength, + input.stayTime, + input.departureEpochInitial, + input.departureEpochSteps, + input.departureEpochStepSize, + input.timeOfFlightMinimum, + input.timeOfFlightSteps, + input.timeOfFlightStepSize ); + + std::cout << "All departure-arrival epochs successfully computed!" << std::endl; + + std::cout << "Computing all pork-chop plot transfers for each leg ... " << std::endl; + + // Compute pork-chop sets for 1-to-1 transfers on a leg-by-leg basis, by stepping through a + // given sequence using recursion. + currentSequencePosition = 0; + + // Compute total set of all Lambert pork-chop data. + AllLambertPorkChopPlots allPorkChopPlots; + int transferId = 1; + recurseLambertTransfers( currentSequencePosition, + tleObjects, + allEpochs, + input.isPrograde, + input.revolutionsMaximum, + sequence, + transferId, + allPorkChopPlots ); + + std::cout << "All pork-chop plot transfers successfully computed! " << std::endl; + + std::cout << "Populating database with pork-chop plot transfers ... " << std::endl; + + // Start SQL transaction to populate database with computed transfers. + SQLite::Transaction transfersTransaction( database ); + + // Set up insert query to add transfer data to table in database. + std::ostringstream lambertScannerTransfersTableInsert; + lambertScannerTransfersTableInsert + << "INSERT INTO lambert_scanner_transfers VALUES (" + << ":transfer_id," << ":departure_object_id," << ":arrival_object_id," << ":departure_epoch," @@ -149,241 +253,274 @@ void executeLambertScanner( const rapidjson::Document& config ) << ":arrival_delta_v_x," << ":arrival_delta_v_y," << ":arrival_delta_v_z," - << ":transfer_delta_v" + << ":transfer_delta_v," + << ":leg_id" << ");"; - SQLite::Statement query( database, lambertScannerTableInsert.str( ) ); - - std::cout << "Computing Lambert transfers and populating database ... " << std::endl; + SQLite::Statement transfersQuery( database, lambertScannerTransfersTableInsert.str( ) ); - // Loop over TLE objects and compute transfers based on Lambert targeter across time-of-flight - // grid. - boost::progress_display showProgress( tleObjects.size( ) ); - - // Loop over all departure objects. - for ( unsigned int i = 0; i < tleObjects.size( ); i++ ) + for ( AllLambertPorkChopPlots::iterator it = allPorkChopPlots.begin( ); + it != allPorkChopPlots.end( ); + it++ ) { - // Compute departure state. - Tle departureObject = tleObjects[ i ]; - SGP4 sgp4Departure( departureObject ); + const PorkChopPlotId porkChopId = it->first; + const LambertPorkChopPlot porkChopPlot = it->second; - DateTime departureEpoch = input.departureEpochInitial; - if ( input.departureEpochInitial == DateTime( ) ) + for ( unsigned int i = 0; i < porkChopPlot.size( ); ++i ) { - departureEpoch = departureObject.Epoch( ); + // Bind values to SQL insert query. + transfersQuery.bind( ":transfer_id", porkChopPlot[ i ].transferId ); + transfersQuery.bind( ":departure_object_id", porkChopId.departureObjectId ); + transfersQuery.bind( ":arrival_object_id", porkChopId.arrivalObjectId ); + transfersQuery.bind( ":departure_epoch", + porkChopPlot[ i ].departureEpoch.ToJulian( ) ); + transfersQuery.bind( ":time_of_flight", porkChopPlot[ i ].timeOfFlight ); + transfersQuery.bind( ":revolutions", porkChopPlot[ i ].revolutions ); + transfersQuery.bind( ":prograde", porkChopPlot[ i ].isPrograde ); + transfersQuery.bind( ":departure_position_x", + porkChopPlot[ i ].departureState[ astro::xPositionIndex ] ); + transfersQuery.bind( ":departure_position_y", + porkChopPlot[ i ].departureState[ astro::yPositionIndex ] ); + transfersQuery.bind( ":departure_position_z", + porkChopPlot[ i ].departureState[ astro::zPositionIndex ] ); + transfersQuery.bind( ":departure_velocity_x", + porkChopPlot[ i ].departureState[ astro::xVelocityIndex ] ); + transfersQuery.bind( ":departure_velocity_y", + porkChopPlot[ i ].departureState[ astro::yVelocityIndex ] ); + transfersQuery.bind( ":departure_velocity_z", + porkChopPlot[ i ].departureState[ astro::zVelocityIndex ] ); + transfersQuery.bind( ":departure_semi_major_axis", + porkChopPlot[ i ].departureStateKepler[ astro::semiMajorAxisIndex ] ); + transfersQuery.bind( ":departure_eccentricity", + porkChopPlot[ i ].departureStateKepler[ astro::eccentricityIndex ] ); + transfersQuery.bind( ":departure_inclination", + porkChopPlot[ i ].departureStateKepler[ astro::inclinationIndex ] ); + transfersQuery.bind( ":departure_argument_of_periapsis", + porkChopPlot[ i ].departureStateKepler[ astro::argumentOfPeriapsisIndex ] ); + transfersQuery.bind( ":departure_longitude_of_ascending_node", + porkChopPlot[ i ].departureStateKepler[ astro::longitudeOfAscendingNodeIndex ] ); + transfersQuery.bind( ":departure_true_anomaly", + porkChopPlot[ i ].departureStateKepler[ astro::trueAnomalyIndex ] ); + transfersQuery.bind( ":arrival_position_x", + porkChopPlot[ i ].arrivalState[ astro::xPositionIndex ] ); + transfersQuery.bind( ":arrival_position_y", + porkChopPlot[ i ].arrivalState[ astro::yPositionIndex ] ); + transfersQuery.bind( ":arrival_position_z", + porkChopPlot[ i ].arrivalState[ astro::zPositionIndex ] ); + transfersQuery.bind( ":arrival_velocity_x", + porkChopPlot[ i ].arrivalState[ astro::xVelocityIndex ] ); + transfersQuery.bind( ":arrival_velocity_y", + porkChopPlot[ i ].arrivalState[ astro::yVelocityIndex ] ); + transfersQuery.bind( ":arrival_velocity_z", + porkChopPlot[ i ].arrivalState[ astro::zVelocityIndex ] ); + transfersQuery.bind( ":arrival_semi_major_axis", + porkChopPlot[ i ].arrivalStateKepler[ astro::semiMajorAxisIndex ] ); + transfersQuery.bind( ":arrival_eccentricity", + porkChopPlot[ i ].arrivalStateKepler[ astro::eccentricityIndex ] ); + transfersQuery.bind( ":arrival_inclination", + porkChopPlot[ i ].arrivalStateKepler[ astro::inclinationIndex ] ); + transfersQuery.bind( ":arrival_argument_of_periapsis", + porkChopPlot[ i ].arrivalStateKepler[ astro::argumentOfPeriapsisIndex ] ); + transfersQuery.bind( ":arrival_longitude_of_ascending_node", + porkChopPlot[ i ].arrivalStateKepler[ astro::longitudeOfAscendingNodeIndex ] ); + transfersQuery.bind( ":arrival_true_anomaly", + porkChopPlot[ i ].arrivalStateKepler[ astro::trueAnomalyIndex ] ); + transfersQuery.bind( ":transfer_semi_major_axis", + porkChopPlot[ i ].transferStateKepler[ astro::semiMajorAxisIndex ] ); + transfersQuery.bind( ":transfer_eccentricity", + porkChopPlot[ i ].transferStateKepler[ astro::eccentricityIndex ] ); + transfersQuery.bind( ":transfer_inclination", + porkChopPlot[ i ].transferStateKepler[ astro::inclinationIndex ] ); + transfersQuery.bind( ":transfer_argument_of_periapsis", + porkChopPlot[ i ].transferStateKepler[ astro::argumentOfPeriapsisIndex ] ); + transfersQuery.bind( ":transfer_longitude_of_ascending_node", + porkChopPlot[ i ].transferStateKepler[ astro::longitudeOfAscendingNodeIndex ] ); + transfersQuery.bind( ":transfer_true_anomaly", + porkChopPlot[ i ].transferStateKepler[ astro::trueAnomalyIndex ] ); + transfersQuery.bind( ":departure_delta_v_x", porkChopPlot[ i ].departureDeltaV[ 0 ] ); + transfersQuery.bind( ":departure_delta_v_y", porkChopPlot[ i ].departureDeltaV[ 1 ] ); + transfersQuery.bind( ":departure_delta_v_z", porkChopPlot[ i ].departureDeltaV[ 2 ] ); + transfersQuery.bind( ":arrival_delta_v_x", porkChopPlot[ i ].arrivalDeltaV[ 0 ] ); + transfersQuery.bind( ":arrival_delta_v_y", porkChopPlot[ i ].arrivalDeltaV[ 1 ] ); + transfersQuery.bind( ":arrival_delta_v_z", porkChopPlot[ i ].arrivalDeltaV[ 2 ] ); + transfersQuery.bind( ":transfer_delta_v", porkChopPlot[ i ].transferDeltaV ); + transfersQuery.bind( ":leg_id", porkChopId.legId ); + + // Execute insert query. + transfersQuery.executeStep( ); + + // Reset SQL insert query. + transfersQuery.reset( ); } + } + + // Commit transfers transaction. + transfersTransaction.commit( ); + + std::cout << "Database successfully populated with pork-chop plot transfers! " << std::endl; - // Loop over arrival objects. - for ( unsigned int j = 0; j < tleObjects.size( ); j++ ) + std::cout << "Computing all multi-leg transfers for each sequence ..." << std::endl; + + AllMultiLegTransfers allMultiLegTransfers; + ListOfMultiLegTransfers listOfMultiLegTransfers; + MultiLegTransferData multiLegTransferData; + + // Loop over all sequences and call recursive function to compute multi-leg transfers. + for ( ListOfSequences::iterator iteratorSequences = listOfSequences.begin( ); + iteratorSequences != listOfSequences.end( ); + iteratorSequences++ ) + { + recurseMuiltiLegLambertTransfers( 0, + iteratorSequences->second, + allPorkChopPlots, + input.stayTime, + listOfMultiLegTransfers, + multiLegTransferData ); + + allMultiLegTransfers[ iteratorSequences->first ] = listOfMultiLegTransfers; + + listOfMultiLegTransfers.clear( ); + multiLegTransferData.clear( ); + } + + std::cout << "All multi-leg transfers successfully computed!" << std::endl; + + std::cout << "Populating the database with all multi-leg transfers ..." << std::endl; + + // Start SQL transaction to populate database with computed multi-leg transfers. + SQLite::Transaction multiLegTransfersTransaction( database ); + + // Set up insert query to add multi-leg transfer data to table in database. + std::ostringstream lambertScannerMultiLegTransfersTableInsert; + lambertScannerMultiLegTransfersTableInsert + << "INSERT INTO lambert_scanner_multi_leg_transfers VALUES (" + << "NULL," + << ":sequence_id,"; + for ( int i = 0; i < input.sequenceLength - 1; ++i ) + { + + lambertScannerMultiLegTransfersTableInsert + << ":transfer_id_" << i + 1 << ","; + } + lambertScannerMultiLegTransfersTableInsert + << ":launch_epoch," + << ":mission_duration," + << ":total_transfer_delta_v" + << ");"; + + SQLite::Statement multiLegTransferQuery( database, + lambertScannerMultiLegTransfersTableInsert.str( ) ); + + for ( AllMultiLegTransfers::iterator iteratorMultiLegTransfers = allMultiLegTransfers.begin( ); + iteratorMultiLegTransfers != allMultiLegTransfers.end( ); + iteratorMultiLegTransfers++ ) + { + const int sequenceId = iteratorMultiLegTransfers->first; + const ListOfMultiLegTransfers listOfMultiLegTransfers = iteratorMultiLegTransfers->second; + + for ( unsigned int i = 0; i < listOfMultiLegTransfers.size( ); ++i ) { - // Skip the case of the departure and arrival objects being the same. - if ( i == j ) - { - continue; - } - Tle arrivalObject = tleObjects[ j ]; - SGP4 sgp4Arrival( arrivalObject ); - const int arrivalObjectId = static_cast< int >( arrivalObject.NoradNumber( ) ); + // Bind values to SQL insert query. + multiLegTransferQuery.bind( ":sequence_id", sequenceId ); + multiLegTransferQuery.bind( + ":launch_epoch", listOfMultiLegTransfers[ i ].launchEpoch.ToJulian( ) ); - // Loop over departure epoch grid. - for ( int m = 0; m < input.departureEpochSteps; ++m ) + double missionDuration = 0.0; + double totalTransferDeltaV = 0.0; + + for ( int j = 0; j < input.sequenceLength - 1; ++j ) { - DateTime departureEpoch = input.departureEpochInitial; - departureEpoch = departureEpoch.AddSeconds( input.departureEpochStepSize * m ); - - 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++ ) - { - 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 ); + std::ostringstream multiLegTransfersQueryTransferNumber; + multiLegTransfersQueryTransferNumber << ":transfer_id_" << j + 1; + multiLegTransferQuery.bind( + multiLegTransfersQueryTransferNumber.str( ), + listOfMultiLegTransfers[ i ].multiLegTransferData[ j ].transferId ); + + missionDuration + += listOfMultiLegTransfers[ i ].multiLegTransferData[ j ].timeOfFlight; + totalTransferDeltaV + += listOfMultiLegTransfers[ i ].multiLegTransferData[ j ].transferDeltaV; + } - 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( arrivalVelocity, - sml::multiply( transferArrivalVelocity, -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 ); + multiLegTransferQuery.bind( ":mission_duration", missionDuration ); + multiLegTransferQuery.bind( ":total_transfer_delta_v", totalTransferDeltaV ); - // 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( ); - } - } - } + // Execute insert query. + multiLegTransferQuery.executeStep( ); - ++showProgress; + // Reset SQL insert query. + multiLegTransferQuery.reset( ); + } } - // Commit transaction. - transaction.commit( ); + // Commit multi-leg transfers transaction. + multiLegTransfersTransaction.commit( ); - std::cout << std::endl; - std::cout << "Database populated successfully!" << std::endl; - std::cout << std::endl; + std::cout << "Database successfully populated with all multi-leg transfers!" << std::endl; + + std::cout << "Update sequences table with best multi-leg transfers ... " << std::endl; - // Check if shortlist file should be created; call function to write output. - if ( input.shortlistLength > 0 ) + std::ostringstream bestMultiLegTransfersReplace; + bestMultiLegTransfersReplace + << "REPLACE INTO sequences " + << "SELECT sequence_id, "; + for ( int i = 0; i < input.sequenceLength; ++i ) + { + bestMultiLegTransfersReplace + << " target_" << i << ", "; + } + for ( int i = 0; i < input.sequenceLength - 1; ++i ) + { + bestMultiLegTransfersReplace + << " transfer_id_" << i + 1 << "_, "; + } + bestMultiLegTransfersReplace + << " launch_epoch_, " + << " lambert_delta_v_, " + << " mission_duration_, " + << " atom_transfer_delta_v " + << "FROM (SELECT * " + << " FROM sequences " + << " AS SEQ " + << " JOIN (SELECT sequence_id AS \"sequence_id_match\", "; + for ( int i = 0; i < input.sequenceLength - 1; ++i ) { - std::cout << "Writing shortlist to file ... " << std::endl; - writeTransferShortlist( database, input.shortlistLength, input.shortlistPath ); - std::cout << "Shortlist file created successfully!" << std::endl; + bestMultiLegTransfersReplace + << " transfer_id_" << i + 1 << " AS transfer_id_" + << i + 1 << "_, "; } + bestMultiLegTransfersReplace + << " launch_epoch AS \"launch_epoch_\", " + << " min(total_transfer_delta_v) AS \"lambert_delta_v_\", " + << " mission_duration AS \"mission_duration_\" " + << " FROM lambert_scanner_multi_leg_transfers " + << " GROUP BY lambert_scanner_multi_leg_transfers.sequence_id) " + << " AS MULTI " + << " ON SEQ.sequence_id = MULTI.sequence_id_match);"; + + + database.exec( bestMultiLegTransfersReplace.str( ).c_str( ) ); + + std::cout << "Sequences table successfully updated with best multi-leg transfers!" << std::endl; + + // Write sequences file. + std::cout << "Writing best multi-leg transfers for each sequence to file ... " << std::endl; + writeSequencesToFile( database, input.sequencesPath, input.sequenceLength ); + std::cout << "Sequences file created successfully!" << std::endl; } //! Check lambert_scanner input parameters. LambertScannerInput checkLambertScannerInput( const rapidjson::Document& config ) { - const std::string catalogPath = find( config, "catalog" )->value.GetString( ); - std::cout << "Catalog " << catalogPath << std::endl; + const std::string catalogPath = find( config, "catalog_path" )->value.GetString( ); + std::cout << "Catalog path " << catalogPath << std::endl; + + const std::string databasePath = find( config, "database_path" )->value.GetString( ); + std::cout << "Database path " << databasePath << std::endl; - const std::string databasePath = find( config, "database" )->value.GetString( ); - std::cout << "Database " << databasePath << std::endl; + const int sequenceLength = find( config, "sequence_length" )->value.GetInt( ); + std::cout << "Sequence length " << sequenceLength << std::endl; const ConfigIterator departureEpochIterator = find( config, "departure_epoch" ); std::map< std::string, int > departureEpochElements; @@ -447,10 +584,10 @@ LambertScannerInput checkLambertScannerInput( const rapidjson::Document& config const double timeOfFlightMinimum = find( config, "time_of_flight_grid" )->value[ 0 ].GetDouble( ); - std::cout << "Minimum Time-of-Flight " << timeOfFlightMinimum << std::endl; + std::cout << "Minimum time-of-flight " << timeOfFlightMinimum << std::endl; const double timeOfFlightMaximum = find( config, "time_of_flight_grid" )->value[ 1 ].GetDouble( ); - std::cout << "Maximum Time-of-Flight " << timeOfFlightMaximum << std::endl; + std::cout << "Maximum time-of-flight " << timeOfFlightMaximum << std::endl; if ( timeOfFlightMinimum > timeOfFlightMaximum ) { @@ -459,7 +596,10 @@ LambertScannerInput checkLambertScannerInput( const rapidjson::Document& config const double timeOfFlightSteps = find( config, "time_of_flight_grid" )->value[ 2 ].GetDouble( ); - std::cout << "# Time-of-Flight steps " << timeOfFlightSteps << std::endl; + std::cout << "# Time-of-flight steps " << timeOfFlightSteps << std::endl; + + const double stayTime = find( config, "stay_time" )->value.GetDouble( ); + std::cout << "Stay time " << stayTime << std::endl; const bool isPrograde = find( config, "is_prograde" )->value.GetBool( ); if ( isPrograde ) @@ -474,18 +614,12 @@ LambertScannerInput checkLambertScannerInput( const rapidjson::Document& config const int revolutionsMaximum = find( config, "revolutions_maximum" )->value.GetInt( ); std::cout << "Maximum revolutions " << revolutionsMaximum << std::endl; - const int shortlistLength = find( config, "shortlist" )->value[ 0 ].GetInt( ); - std::cout << "# of shortlist transfers " << shortlistLength << std::endl; - - std::string shortlistPath = ""; - if ( shortlistLength > 0 ) - { - shortlistPath = find( config, "shortlist" )->value[ 1 ].GetString( ); - std::cout << "Shortlist " << shortlistPath << std::endl; - } + const std::string sequencesPath = find( config, "sequences_path" )->value.GetString( ); + std::cout << "Sequences path " << sequencesPath << std::endl; return LambertScannerInput( catalogPath, databasePath, + sequenceLength, departureEpoch, departureGridSteps, departureEpochRange/departureGridSteps, @@ -493,25 +627,56 @@ LambertScannerInput checkLambertScannerInput( const rapidjson::Document& config timeOfFlightMaximum, timeOfFlightSteps, ( timeOfFlightMaximum - timeOfFlightMinimum ) / timeOfFlightSteps, + stayTime, isPrograde, revolutionsMaximum, - shortlistLength, - shortlistPath ); + sequencesPath ); } -//! Create lambert_scanner table. -void createLambertScannerTable( SQLite::Database& database ) +//! Create lambert_scanner tables. +void createLambertScannerTables( SQLite::Database& database, const int sequenceLength ) { // Drop table from database if it exists. - database.exec( "DROP TABLE IF EXISTS lambert_scanner_results;" ); - - // Set up SQL command to create table to store lambert_scanner results. - std::ostringstream lambertScannerTableCreate; - lambertScannerTableCreate - << "CREATE TABLE lambert_scanner_results (" - << "\"transfer_id\" INTEGER PRIMARY KEY AUTOINCREMENT," - << "\"departure_object_id\" TEXT," - << "\"arrival_object_id\" TEXT," + database.exec( "DROP TABLE IF EXISTS sequences;" ); + database.exec( "DROP TABLE IF EXISTS lambert_scanner_transfers;" ); + database.exec( "DROP TABLE IF EXISTS lambert_scanner_multi_leg_transfers;" ); + + // Set up SQL command to create table to store lambert_sequences. + std::ostringstream lambertScannerSequencesTableCreate; + lambertScannerSequencesTableCreate + << "CREATE TABLE sequences (" + << "\"sequence_id\" INTEGER PRIMARY KEY ,"; + for ( int i = 0; i < sequenceLength; ++i ) + { + lambertScannerSequencesTableCreate + << "\"target_" << i << "\" INTEGER ,"; + } + for ( int i = 0; i < sequenceLength - 1; ++i ) + { + lambertScannerSequencesTableCreate + << "\"transfer_id_" << i + 1 << "\" INTEGER ,"; + } + lambertScannerSequencesTableCreate + << "\"launch_epoch\" REAL ," + << "\"lambert_transfer_delta_v\" REAL ," + << "\"mission_duration\" REAL ," + << "\"atom_transfer_delta_v\" REAL );"; + + // // Execute command to create table. + database.exec( lambertScannerSequencesTableCreate.str( ).c_str( ) ); + + if ( !database.tableExists( "sequences" ) ) + { + throw std::runtime_error( "ERROR: Creating table 'sequences' failed!" ); + } + + // Set up SQL command to create table to store lambert_scanner transfers. + std::ostringstream lambertScannerTransfersTableCreate; + lambertScannerTransfersTableCreate + << "CREATE TABLE lambert_scanner_transfers (" + << "\"transfer_id\" INTEGER PRIMARY KEY," + << "\"departure_object_id\" INTEGER," + << "\"arrival_object_id\" INTEGER," << "\"departure_epoch\" REAL," << "\"time_of_flight\" REAL," << "\"revolutions\" INTEGER," @@ -553,182 +718,396 @@ void createLambertScannerTable( SQLite::Database& database ) << "\"arrival_delta_v_x\" REAL," << "\"arrival_delta_v_y\" REAL," << "\"arrival_delta_v_z\" REAL," - << "\"transfer_delta_v\" REAL" + << "\"transfer_delta_v\" REAL," + << "\"leg_id\" INTEGER" << ");"; // Execute command to create table. - database.exec( lambertScannerTableCreate.str( ).c_str( ) ); + database.exec( lambertScannerTransfersTableCreate.str( ).c_str( ) ); // Execute command to create index on transfer Delta-V column. std::ostringstream transferDeltaVIndexCreate; transferDeltaVIndexCreate << "CREATE INDEX IF NOT EXISTS \"transfer_delta_v\" on " - << "lambert_scanner_results (transfer_delta_v ASC);"; + << "lambert_scanner_transfers (transfer_delta_v ASC);"; database.exec( transferDeltaVIndexCreate.str( ).c_str( ) ); - if ( !database.tableExists( "lambert_scanner_results" ) ) + if ( !database.tableExists( "lambert_scanner_transfers" ) ) + { + throw std::runtime_error( "ERROR: Creating table 'lambert_scanner_transfers' failed!" ); + } + + // Set up SQL command to create table to store lambert_scanner multi-leg transfers. + std::ostringstream lambertScannerMultiLegTransfersTableCreate; + lambertScannerMultiLegTransfersTableCreate + << "CREATE TABLE lambert_scanner_multi_leg_transfers (" + << "\"multi_leg_transfer_id\" INTEGER PRIMARY KEY AUTOINCREMENT," + << "\"sequence_id\" INTEGER,"; + for ( int i = 0; i < sequenceLength - 1; ++i ) + { + lambertScannerMultiLegTransfersTableCreate + << "\"transfer_id_" << i + 1 << "\" INTEGER,"; + } + lambertScannerMultiLegTransfersTableCreate + << "\"launch_epoch\" REAL," + << "\"mission_duration\" REAL," + << "\"total_transfer_delta_v\" REAL" + << ");"; + + // Execute command to create table. + database.exec( lambertScannerMultiLegTransfersTableCreate.str( ).c_str( ) ); + + // Execute command to create index on total transfer Delta-V column. + std::ostringstream totalTransferDeltaVIndexCreate; + totalTransferDeltaVIndexCreate << "CREATE INDEX IF NOT EXISTS \"total_transfer_delta_v\" on " + << "lambert_scanner_multi_leg_transfers " + << "(total_transfer_delta_v ASC);"; + database.exec( totalTransferDeltaVIndexCreate.str( ).c_str( ) ); + + if ( !database.tableExists( "lambert_scanner_multi_leg_transfers" ) ) { - throw std::runtime_error( "ERROR: Creating table 'lambert_scanner_results' failed!" ); + throw std::runtime_error( + "ERROR: Creating table 'lambert_scanner_multi_leg_transfers' failed!" ); } } -//! Write transfer shortlist to file. -void writeTransferShortlist( SQLite::Database& database, - const int shortlistNumber, - const std::string& shortlistPath ) +//! Write best multi-leg Lambert transfers for each sequence to file. +void writeSequencesToFile( SQLite::Database& database, + const std::string& sequencesPath, + const int sequenceLength ) { - // Fetch transfers to include in shortlist. - // Database table is sorted by transfer_delta_v, in ascending order. - std::ostringstream shortlistSelect; - shortlistSelect << "SELECT * FROM lambert_scanner_results ORDER BY transfer_delta_v ASC LIMIT " - << shortlistNumber << ";"; - SQLite::Statement query( database, shortlistSelect.str( ) ); + // Fetch sequences tables from database and sort from lowest to highest Delta-V. + std::ostringstream sequencesSelect; + sequencesSelect << "SELECT * FROM sequences ORDER BY lambert_transfer_delta_v ASC;"; + SQLite::Statement query( database, sequencesSelect.str( ) ); - // Write fetch data to file. - std::ofstream shortlistFile( shortlistPath.c_str( ) ); + // Write sequences to file. + std::ofstream sequencesFile( sequencesPath.c_str( ) ); // Print file header. - shortlistFile << "transfer_id," - << "departure_object_id," - << "arrival_object_id," - << "departure_epoch," - << "time_of_flight," - << "revolutions," - << "prograde," - << "departure_position_x," - << "departure_position_y," - << "departure_position_z," - << "departure_velocity_x," - << "departure_velocity_y," - << "departure_velocity_z," - << "departure_semi_major_axis," - << "departure_eccentricity," - << "departure_inclination," - << "departure_argument_of_periapsis," - << "departure_longitude_of_ascending_node," - << "departure_true_anomaly," - << "arrival_position_x," - << "arrival_position_y," - << "arrival_position_z," - << "arrival_velocity_x," - << "arrival_velocity_y," - << "arrival_velocity_z," - << "arrival_semi_major_axis," - << "arrival_eccentricity," - << "arrival_inclination," - << "arrival_argument_of_periapsis," - << "arrival_longitude_of_ascending_node," - << "arrival_true_anomaly," - << "transfer_semi_major_axis," - << "transfer_eccentricity," - << "transfer_inclination," - << "transfer_argument_of_periapsis," - << "transfer_longitude_of_ascending_node," - << "transfer_true_anomaly," - << "departure_delta_v_x," - << "departure_delta_v_y," - << "departure_delta_v_z," - << "arrival_delta_v_x," - << "arrival_delta_v_y," - << "arrival_delta_v_z," - << "transfer_delta_v" + sequencesFile << "sequence_id,"; + for ( int i = 0; i < sequenceLength; ++i ) + { + sequencesFile << "target_" << i << ","; + } + for ( int i = 0; i < sequenceLength - 1; ++i ) + { + sequencesFile << "transfer_id_" << i + 1 << ","; + } + sequencesFile << "launch_epoch," + << "lambert_transfer_delta_v," + << "mission_duration" << std::endl; // Loop through data retrieved from database and write to file. while( query.executeStep( ) ) { - const int transferId = query.getColumn( 0 ); - const int departureObjectId = query.getColumn( 1 ); - const int arrivalObjectId = query.getColumn( 2 ); - const double departureEpoch = query.getColumn( 3 ); - const double timeOfFlight = query.getColumn( 4 ); - const int revolutions = query.getColumn( 5 ); - const int prograde = query.getColumn( 6 ); - const double departurePositionX = query.getColumn( 7 ); - const double departurePositionY = query.getColumn( 8 ); - const double departurePositionZ = query.getColumn( 9 ); - const double departureVelocityX = query.getColumn( 10 ); - const double departureVelocityY = query.getColumn( 11 ); - const double departureVelocityZ = query.getColumn( 12 ); - const double departureSemiMajorAxis = query.getColumn( 13 ); - const double departureEccentricity = query.getColumn( 14 ); - const double departureInclination = query.getColumn( 15 ); - const double departureArgumentOfPeriapsis = query.getColumn( 16 ); - const double departureLongitudeOfAscendingNode = query.getColumn( 17 ); - const double departureTrueAnomaly = query.getColumn( 18 ); - const double arrivalPositionX = query.getColumn( 19 ); - const double arrivalPositionY = query.getColumn( 20 ); - const double arrivalPositionZ = query.getColumn( 21 ); - const double arrivalVelocityX = query.getColumn( 22 ); - const double arrivalVelocityY = query.getColumn( 23 ); - const double arrivalVelocityZ = query.getColumn( 24 ); - const double arrivalSemiMajorAxis = query.getColumn( 25 ); - const double arrivalEccentricity = query.getColumn( 26 ); - const double arrivalInclination = query.getColumn( 27 ); - const double arrivalArgumentOfPeriapsis = query.getColumn( 28 ); - const double arrivalLongitudeOfAscendingNode = query.getColumn( 29 ); - const double arrivalTrueAnomaly = query.getColumn( 30 ); - const double transferSemiMajorAxis = query.getColumn( 31 ); - const double transferEccentricity = query.getColumn( 32 ); - const double transferInclination = query.getColumn( 33 ); - const double transferArgumentOfPeriapsis = query.getColumn( 34 ); - const double transferLongitudeOfAscendingNode = query.getColumn( 35 ); - const double transferTrueAnomaly = query.getColumn( 36 ); - const double departureDeltaVX = query.getColumn( 37 ); - const double departureDeltaVY = query.getColumn( 38 ); - const double departureDeltaVZ = query.getColumn( 39 ); - const double arrivalDeltaVX = query.getColumn( 40 ); - const double arrivalDeltaVY = query.getColumn( 41 ); - const double arrivalDeltaVZ = query.getColumn( 42 ); - const double transferDeltaV = query.getColumn( 43 ); - - shortlistFile << transferId << "," - << departureObjectId << "," - << arrivalObjectId << "," - << departureEpoch << "," - << timeOfFlight << "," - << revolutions << "," - << prograde << "," - << departurePositionX << "," - << departurePositionY << "," - << departurePositionZ << "," - << departureVelocityX << "," - << departureVelocityY << "," - << departureVelocityZ << "," - << departureSemiMajorAxis << "," - << departureEccentricity << "," - << departureInclination << "," - << departureArgumentOfPeriapsis << "," - << departureLongitudeOfAscendingNode << "," - << departureTrueAnomaly << "," - << arrivalPositionX << "," - << arrivalPositionY << "," - << arrivalPositionZ << "," - << arrivalVelocityX << "," - << arrivalVelocityY << "," - << arrivalVelocityZ << "," - << arrivalSemiMajorAxis << "," - << arrivalEccentricity << "," - << arrivalInclination << "," - << arrivalArgumentOfPeriapsis << "," - << arrivalLongitudeOfAscendingNode << "," - << arrivalTrueAnomaly << "," - << transferSemiMajorAxis << "," - << transferEccentricity << "," - << transferInclination << "," - << transferArgumentOfPeriapsis << "," - << transferLongitudeOfAscendingNode << "," - << transferTrueAnomaly << "," - << departureDeltaVX << "," - << departureDeltaVY << "," - << departureDeltaVZ << "," - << arrivalDeltaVX << "," - << arrivalDeltaVY << "," - << arrivalDeltaVZ << "," - << transferDeltaV + const int sequenceId = query.getColumn( 0 ); + std::vector< int > targets( sequenceLength ); + for ( unsigned int i = 0; i < targets.size( ); ++i ) + { + targets[ i ] = query.getColumn( i + 1 ); + } + std::vector< int > transferIds( sequenceLength - 1 ); + for ( unsigned int i = 0; i < transferIds.size( ); ++i ) + { + transferIds[ i ] = query.getColumn( i + sequenceLength + 1 ); + } + const double launchEpoch = query.getColumn( 2 * sequenceLength ); + const double lambertTransferDeltaV = query.getColumn( 2 * sequenceLength + 1 ); + const double missionDuration = query.getColumn( 2 * sequenceLength + 2 ); + + sequencesFile << sequenceId << ","; + for ( unsigned int i = 0; i < targets.size( ); ++i ) + { + sequencesFile << targets[ i ] << ","; + } + for ( unsigned int i = 0; i < transferIds.size( ); ++i ) + { + sequencesFile << transferIds[ i ] << ","; + } + sequencesFile << launchEpoch << "," + << lambertTransferDeltaV << "," + << missionDuration << std::endl; } - shortlistFile.close( ); + sequencesFile.close( ); +} + +//! Recurse through sequences leg-by-leg and compute pork-chop plots. +void recurseLambertTransfers( const int currentSequencePosition, + const TleObjects& tleObjects, + const AllEpochs& allEpochs, + const bool isPrograde, + const int revolutionsMaximum, + Sequence& sequence, + int& transferId, + AllLambertPorkChopPlots& allPorkChopPlots ) +{ + // Check if the end of the sequence has been reached. + if ( currentSequencePosition == static_cast< int >( sequence.size( ) ) ) + { + return; + } + + for ( unsigned int i = 0; i < tleObjects.size( ); i++ ) + { + sequence[ currentSequencePosition ] = tleObjects[ i ]; + + // Check if the current position in the sequence is beyond the first. + // (currentSequencePosition = 0 means that the first TLE object ID has to be set). + if ( currentSequencePosition > 0 ) + { + const Tle departureObject = sequence[ currentSequencePosition - 1 ]; + const Tle arrivalObject = sequence[ currentSequencePosition ]; + const int currentLeg = currentSequencePosition; + const PorkChopPlotId porkChopPlotId( currentLeg, + departureObject.NoradNumber( ), + arrivalObject.NoradNumber( ) ); + + // Check if pork-chop ID exists already, i.e., combination of current leg, departure + // object ID and arrival object ID. + // If it exists already, skip the generation of the pork-chop plot. + if ( allPorkChopPlots.find( porkChopPlotId ) == allPorkChopPlots.end( ) ) + { + allPorkChopPlots[ porkChopPlotId ] + = computeLambertPorkChopPlot( departureObject, + arrivalObject, + allEpochs.at( currentLeg ), + isPrograde, + revolutionsMaximum, + transferId ); + } + } + + // Erase the selected TLE object at position currentSequencePosition in sequence. + TleObjects tleObjectsLocal = tleObjects; + tleObjectsLocal.erase( tleObjectsLocal.begin( ) + i ); + + recurseLambertTransfers( currentSequencePosition + 1, + tleObjectsLocal, + allEpochs, + isPrograde, + revolutionsMaximum, + sequence, + transferId, + allPorkChopPlots ); + } +} + +//! Compute Lambert pork-chop plot. +LambertPorkChopPlot computeLambertPorkChopPlot( const Tle& departureObject, + const Tle& arrivalObject, + const ListOfEpochs& listOfEpochs, + const bool isPrograde, + const int revolutionsMaximum, + int& transferId ) +{ + LambertPorkChopPlot porkChopPlot; + + const double earthGravitationalParameter = kMU; + + // Loop over all departure and arrival epoch pairs and compute transfers. + for ( unsigned int i = 0; i < listOfEpochs.size( ); ++i ) + { + DateTime departureEpoch = listOfEpochs[ i ].first; + DateTime arrivalEpoch = listOfEpochs[ i ].second; + + // Compute Cartesian departure position and velocity. + SGP4 sgp4Departure( departureObject ); + 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 ); + + // Compute Cartesian departure position and velocity. + SGP4 sgp4Arrival( arrivalObject ); + 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 ); + + const double timeOfFlight + = ( arrivalEpoch.ToJulian( ) - departureEpoch.ToJulian( ) ) * 24.0 * 3600.0; + + kep_toolbox::lambert_problem targeter( departurePosition, + arrivalPosition, + timeOfFlight, + earthGravitationalParameter, + !isPrograde, + 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 j = 0; j < numberOfSolutions; ++j ) + { + // Compute Delta-V for transfer. + const Vector3 transferDepartureVelocity = targeter.get_v1( )[ j ]; + const Vector3 transferArrivalVelocity = targeter.get_v2( )[ j ]; + + departureDeltaVs[ j ] = sml::add( transferDepartureVelocity, + sml::multiply( departureVelocity, -1.0 ) ); + arrivalDeltaVs[ j ] = sml::add( arrivalVelocity, + sml::multiply( transferArrivalVelocity, -1.0 ) ); + + transferDeltaVs[ j ] + = sml::norm< double >( departureDeltaVs[ j ] ) + + sml::norm< double >( arrivalDeltaVs[ j ] ); + } + + 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 ); + + porkChopPlot.push_back( + LambertPorkChopPlotGridPoint( transferId, + departureEpoch, + arrivalEpoch, + timeOfFlight, + revolutions, + isPrograde, + departureState, + departureStateKepler, + arrivalState, + arrivalStateKepler, + transferStateKepler, + departureDeltaVs[ minimumDeltaVIndex ], + arrivalDeltaVs[ minimumDeltaVIndex ], + *minimumDeltaVIterator ) ); + transferId++; + + } + return porkChopPlot; +} + +//! Recurse through sequence leg-by-leg and compute multi-leg transfers. +void recurseMuiltiLegLambertTransfers( const int currentSequencePosition, + const Sequence& sequence, + AllLambertPorkChopPlots& allPorkChopPlots, + const double stayTime, + ListOfMultiLegTransfers& listOfMultiLegTransfers, + MultiLegTransferData& multiLegTransferData, + DateTime launchEpoch, + DateTime lastArrivalEpoch ) +{ + // Check if the end of the sequence has been reached. + if ( currentSequencePosition + 1 == static_cast< int >( sequence.size( ) ) ) + { + // Add transfer to list. + listOfMultiLegTransfers.push_back( MultiLegTransfer( launchEpoch, multiLegTransferData ) ); + + return; + } + + // Extract pork-chop plot based on current position in sequence. + const PorkChopPlotId porkChopPlotId( currentSequencePosition + 1, + sequence[ currentSequencePosition ].NoradNumber( ), + sequence[ currentSequencePosition + 1 ].NoradNumber( ) ); + LambertPorkChopPlot porkChopPlot = allPorkChopPlots[ porkChopPlotId ]; + LambertPorkChopPlot porkChopPlotMatched; + + // If we are beyond the first leg, the pork-chop plot needs to be filtered to ensure that only + // the departure epochs that match the last arrival epoch + stay time are retained. + if ( currentSequencePosition > 0 ) + { + // Declare list of indices in pork-chop plot list to keep based on matching the arrival + // epoch from the previous leg with the departure epochs for the current leg. + std::vector< int > matchIndices; + + // Loop over pork-chop plot list and save indices where the last arrival epoch + stay time + // matches the departure epochs for the current leg. + for ( unsigned int i = 0; i < porkChopPlot.size( ); ++i ) + { + DateTime matchEpoch = lastArrivalEpoch; + matchEpoch.AddSeconds( stayTime ); + if ( porkChopPlot[ i ].departureEpoch == matchEpoch ) + { + matchIndices.push_back( i ); + } + } + + // Loop over the matched indices and copy the pork-chop plot grid points to the matched + // list. + for ( unsigned int i = 0; i < matchIndices.size( ); ++i ) + { + porkChopPlotMatched.push_back( porkChopPlot[ matchIndices[ i ] ] ); + } + } + else + { + porkChopPlotMatched = porkChopPlot; + } + + // Loop over all departure-arrival epoch pairs in pork-chop plot. + for ( unsigned int i = 0; i < porkChopPlotMatched.size( ); ++i ) + { + // If the current position is the start of the first leg, set the launch epoch to the + // departure epoch for the multi-leg transfer. + if ( currentSequencePosition == 0 ) + { + launchEpoch = porkChopPlotMatched[ i ].departureEpoch; + } + + // Add the time-of-flight and Delta V of the current transfer to the data container. + multiLegTransferData.push_back( TransferData( porkChopPlotMatched[ i ].transferId, + porkChopPlotMatched[ i ].timeOfFlight, + porkChopPlotMatched[ i ].transferDeltaV ) ); + + // Update the last arrival epoch to the current arrival epoch. + lastArrivalEpoch = porkChopPlotMatched[ i ].arrivalEpoch; + + recurseMuiltiLegLambertTransfers( currentSequencePosition + 1, + sequence, + allPorkChopPlots, + stayTime, + listOfMultiLegTransfers, + multiLegTransferData, + launchEpoch, + lastArrivalEpoch ); + + // Remove last entry in multi-leg transfer; + multiLegTransferData.pop_back( ); + } } } // namespace d2d diff --git a/src/sgp4Scanner.cpp b/src/sgp4Scanner.cpp index dc53987..1738731 100644 --- a/src/sgp4Scanner.cpp +++ b/src/sgp4Scanner.cpp @@ -27,6 +27,8 @@ #include #include +#include + #include #include diff --git a/src/tools.cpp b/src/tools.cpp index 1aa83c6..e258a6a 100644 --- a/src/tools.cpp +++ b/src/tools.cpp @@ -5,25 +5,19 @@ * See accompanying file LICENSE.md or copy at http://opensource.org/licenses/MIT */ -#include #include #include #include #include -#include -#include #include #include #include -#include - -#include "D2D/tools.hpp" - -#include #include +#include "D2D/tools.hpp" + namespace d2d { @@ -97,7 +91,7 @@ bool executeVirtualTleConvergenceTest( const Vector6& propagatedCartesianState, const Vector6& trueCartesianState, const double relativeTolerance, const double absoluteTolerance ) -{ +{ // Check for NAN values. for ( int i = 0; i < 6; i++ ) { @@ -109,20 +103,20 @@ bool executeVirtualTleConvergenceTest( const Vector6& propagatedCartesianState, } } - // Check if error between target and propagated Cartesian states is within specified + // Check if error between target and propagated Cartesian states is within specified // tolerances. Vector6 absoluteDifference = propagatedCartesianState; Vector6 relativeDifference = propagatedCartesianState; for ( int i = 0; i < 6; i++ ) { - absoluteDifference[ i ] - = std::fabs( propagatedCartesianState[ i ] - trueCartesianState[ i ] ); + absoluteDifference[ i ] + = std::fabs( propagatedCartesianState[ i ] - trueCartesianState[ i ] ); relativeDifference[ i ] = absoluteDifference[ i ] / std::fabs( trueCartesianState[ i ] ); if ( relativeDifference[ i ] > relativeTolerance ) { if ( absoluteDifference[ i ] > absoluteTolerance ) { - // If the relative and absolute difference for the ith element exceeds the + // If the relative and absolute difference for the ith element exceeds the // specified tolerances, the convergence test has failed. return false; } @@ -210,4 +204,202 @@ int getTleCatalogType( const std::string& catalogFirstLine ) return tleLines; } +//! Recurse leg-by-leg to generate list of TLE sequences. +void recurseSequences( const int currentSequencePosition, + const TleObjects& tleObjects, + Sequence& sequence, + int& sequenceId, + ListOfSequences& listOfSequences ) +{ + // If current leg has reached the length of the sequence, then the sequence is complete. + if ( currentSequencePosition == static_cast< int >( sequence.size( ) ) ) + { + return; + } + + // Loop through pool of TLE objects and store IDs in sequence. + for ( unsigned int i = 0; i < tleObjects.size( ); i++ ) + { + // Store ith TLE object in sequence at position of current leg. + sequence[ currentSequencePosition ] = tleObjects[ i ]; + + // Create a local copy of the list of TLE objects and erase the object selected for the + // current leg. + TleObjects tleObjectsLocal = tleObjects; + tleObjectsLocal.erase( tleObjectsLocal.begin( ) + i ); + + // Proceed to the next leg in the sequence through recursion. + recurseSequences( currentSequencePosition + 1, + tleObjectsLocal, + sequence, + sequenceId, + listOfSequences ); + + // Write the sequence to the list. + if ( currentSequencePosition == static_cast< int >( sequence.size( ) ) - 1 ) + { + listOfSequences[ sequenceId ] = sequence; + sequenceId++; + } + } +} + +//! Compute departure-arrival epoch pairs for all pork-chop plots. +AllEpochs computeAllPorkChopPlotEpochs( const int sequenceLength, + const double stayTime, + const DateTime& departureEpochInitial, + const int departureEpochSteps, + const double departureEpochStepSize, + const double timeOfFlightMinimum, + const int timeOfFlightSteps, + const double timeOfFlightStepSize ) +{ + // Create container of unique departure epochs that are used per leg to compute the pork-chop + // data set. + std::vector< DateTime > uniqueDepartureEpochs; + uniqueDepartureEpochs.push_back( departureEpochInitial ); + for ( int i = 1; i < departureEpochSteps + 1; ++i ) + { + DateTime departureEpoch = departureEpochInitial; + departureEpoch = departureEpoch.AddSeconds( departureEpochStepSize * i ); + uniqueDepartureEpochs.push_back( departureEpoch ); + } + + AllEpochs allEpochs; + + // Loop over each leg and generate the departure-arrival epoch pairs. + for ( int i = 0; i < sequenceLength - 1; ++i ) + { + ListOfEpochs listOfEpochs; + + // Loop over unique departure epochs. + for ( unsigned int j = 0; j < uniqueDepartureEpochs.size( ); ++j ) + { + // Loop over time-of-flight grid. + for ( int k = 0; k < timeOfFlightSteps + 1; k++ ) + { + const double timeOfFlight = timeOfFlightMinimum + k * timeOfFlightStepSize; + const DateTime departureEpoch = uniqueDepartureEpochs[ j ]; + const DateTime arrivalEpoch = uniqueDepartureEpochs[ j ].AddSeconds( timeOfFlight ); + + // Store pair of departure and arrival epochs in the list. + listOfEpochs.push_back( std::make_pair< DateTime, DateTime >( departureEpoch, + arrivalEpoch ) ); + } + } + + // Store list of all departure and arrival epoch combinations for the current leg in the + // map. + allEpochs[ i + 1 ] = listOfEpochs; + + // Extract all arrival epochs from list of epochs. + std::vector< DateTime > listOfArrivalEpochs; + for ( unsigned int j = 0; j < listOfEpochs.size( ); ++j ) + { + listOfArrivalEpochs.push_back( listOfEpochs[ j ].second ); + } + + // Sort arrival epochs and only keep unique entries. The unique arrival epochs are saved as + // the departure epochs for the next leg after adding a fixed stay time. + std::sort( listOfArrivalEpochs.begin( ), listOfArrivalEpochs.end( ) ); + std::vector< DateTime >::iterator arrivalEpochIterator + = std::unique( listOfArrivalEpochs.begin( ), listOfArrivalEpochs.end( ) ); + listOfArrivalEpochs.resize( std::distance( listOfArrivalEpochs.begin( ), + arrivalEpochIterator ) ); + for ( unsigned int j = 0; j < listOfArrivalEpochs.size( ); ++j ) + { + listOfArrivalEpochs[ j ] = listOfArrivalEpochs[ j ].AddSeconds( stayTime ); + } + uniqueDepartureEpochs = listOfArrivalEpochs; + } + + return allEpochs; +} + +//! Overload ==-operator to compare PorkChopPlotId objects. +bool operator==( const PorkChopPlotId& id1, const PorkChopPlotId& id2 ) +{ + bool isEqual = false; + if ( id1.legId == id2.legId ) + { + if ( id1.departureObjectId == id2.departureObjectId ) + { + if ( id1.arrivalObjectId == id2.arrivalObjectId ) + { + isEqual = true; + } + } + } + + return isEqual; +} + +//! Overload !=-operator to compare PorkChopPlotId objects. +bool operator!=( const PorkChopPlotId& id1, const PorkChopPlotId& id2 ) +{ + return !(id1 == id2); +} + +//! Overload <-operator to compare PorkChopPlotId objects. +bool operator<( const PorkChopPlotId& id1, const PorkChopPlotId& id2 ) +{ + bool isLessThan = false; + if ( id1.legId == id2.legId ) + { + if ( id1.departureObjectId == id2.departureObjectId ) + { + if ( id1.arrivalObjectId < id2.arrivalObjectId ) + { + isLessThan = true; + } + } + else if ( id1.departureObjectId < id2.departureObjectId ) + { + isLessThan = true; + } + } + else if ( id1.legId < id2.legId ) + { + isLessThan = true; + } + + return isLessThan; +} + +//! Overload >=-operator to compare PorkChopPlotId objects. +bool operator<=( const PorkChopPlotId& id1, const PorkChopPlotId& id2 ) +{ + bool isLessThanOrEqual = false; + if ( id1 < id2 ) + { + isLessThanOrEqual = true; + } + else if ( id1 == id2 ) + { + isLessThanOrEqual = true; + } + return isLessThanOrEqual; +} + +//! Overload >-operator to compare PorkChopPlotId objects. +bool operator>( const PorkChopPlotId& id1, const PorkChopPlotId& id2 ) +{ + return !( id1 <= id2 ); +} + +//! Overload >=-operator to compare PorkChopPlotId objects. +bool operator>=( const PorkChopPlotId& id1, const PorkChopPlotId& id2 ) +{ + bool isGreaterThanOrEqual = false; + if ( id1 > id2 ) + { + isGreaterThanOrEqual = true; + } + else if ( id1 == id2 ) + { + isGreaterThanOrEqual = true; + } + return isGreaterThanOrEqual; +} + } // namespace d2d