1 /***************************************************************************** 2 * Copyright (C) 2004-2018 The pykep development team, * 3 * Advanced Concepts Team (ACT), European Space Agency (ESA) * 4 * * 5 * https://gitter.im/esa/pykep * 6 * https://github.com/esa/pykep * 7 * * 8 * act@esa.int * 9 * * 10 * This program is free software; you can redistribute it and/or modify * 11 * it under the terms of the GNU General Public License as published by * 12 * the Free Software Foundation; either version 2 of the License, or * 13 * (at your option) any later version. * 14 * * 15 * This program is distributed in the hope that it will be useful, * 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 18 * GNU General Public License for more details. * 19 * * 20 * You should have received a copy of the GNU General Public License * 21 * along with this program; if not, write to the * 22 * Free Software Foundation, Inc., * 23 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * 24 *****************************************************************************/ 25 26 #ifndef KEP_TOOLBOX_LEG_H 27 #define KEP_TOOLBOX_LEG_H 28 29 #include <boost/type_traits/is_same.hpp> 30 #include <boost/utility.hpp> 31 #include <iterator> 32 #include <vector> 33 34 #include <keplerian_toolbox/core_functions/array3D_operations.hpp> 35 #include <keplerian_toolbox/core_functions/propagate_lagrangian.hpp> 36 #include <keplerian_toolbox/core_functions/propagate_taylor.hpp> 37 #include <keplerian_toolbox/detail/visibility.hpp> 38 #include <keplerian_toolbox/epoch.hpp> 39 #include <keplerian_toolbox/exceptions.hpp> 40 #include <keplerian_toolbox/serialization.hpp> 41 #include <keplerian_toolbox/sims_flanagan/sc_state.hpp> 42 #include <keplerian_toolbox/sims_flanagan/spacecraft.hpp> 43 #include <keplerian_toolbox/sims_flanagan/throttle.hpp> 44 45 namespace kep_toolbox 46 { 47 /// Sims-Flanagan transcription of low-thrust trajectories 48 /** 49 * This namespace contains the routines that allow building and evaluating low-thrust trajectories using the 50 * Sims-Flanagan transcription method. 51 */ 52 namespace sims_flanagan 53 { 54 55 /// Single low-thrust leg (phase) 56 /** 57 * This class represents, generically, a low-thrust leg (phase) as a sequence of successive 58 * impulses of magnitude compatible with the low-thrust propulsion system of a spacecraft. 59 * The leg achieves to transfer a given spacecraft from an initial to a final state in the 60 * time given (and can be considered as feasible) whenever the method evaluate_mismatch 61 * returns all zeros and the method get_throttles_con returns all values less than zero. 62 * Th sequence of different impulses is represented by the class throttles. These represent 63 * the cartesian components \f$ \mathbf x = (x_1,y_1,z_1) \f$ of a normalized \f$ \Delta V \f$ and are thus 64 * numbers that need to satisfy the constraint \f$|\mathbf x| \le 1\f$ 65 * 66 * \image html sims_flanagan_leg.png "Visualization of a feasible leg (Earth-Mars)" 67 * \image latex sims_flanagan_leg.png "Visualization of a feasible leg (Earth-Mars)" width=5cm 68 * 69 * @author Dario Izzo (dario.izzo _AT_ googlemail.com) 70 */ 71 class KEP_TOOLBOX_DLL_PUBLIC leg 72 { 73 friend std::ostream &operator<<(std::ostream &s, const leg &in); 74 75 public: 76 std::string human_readable() const; 77 /// Constructor. 78 /** 79 * Default constructor. Constructs a meaningless leg that will need to be properly initialized 80 * using the various setters.... 81 */ leg()82 leg() : t_i(), x_i(), throttles(), t_f(), x_f(), m_sc(), m_mu(0), m_hf(false), m_tol(-10) {} 83 84 /// Constructs the leg from epochs, sc_states and cartesian components of throttles 85 /** 86 * Constructs entirely a leg assuming high fidelity propagation switched off and equally spaced segments 87 * 88 * \param[in] epoch_i Inital epoch 89 * \param[in] state_i Initial sc_state (spacecraft state) 90 * \param[in] throttle sequence of doubles (\f$ x_1,y_1,z_1, ..., x_N,y_N,z_N \f$) representing the cartesian 91 * components of the throttles 92 * \param[in] state_f Final sc_state (spacecraft state) 93 * \param[in] sc Spacecraft object 94 * \param[in] mu Primary body gravitational constant 95 * 96 */ leg(const epoch & epoch_i,const sc_state & state_i,const std::vector<double> & thrott,const epoch & epoch_f,const sc_state & state_f,const spacecraft & sc,const double mu)97 leg(const epoch &epoch_i, const sc_state &state_i, const std::vector<double> &thrott, const epoch &epoch_f, 98 const sc_state &state_f, const spacecraft &sc, const double mu) 99 : m_sc(sc), m_hf(false), m_tol(-10) 100 { 101 set_leg(epoch_i, state_i, thrott.begin(), thrott.end(), epoch_f, state_f, mu); 102 } 103 104 /// Initialize a leg 105 /** 106 * Initialize a leg assuming that the user has or will initialize separately the spacecraft. 107 * The throttles are provided via two iterators pointing 108 * to the beginning and to the end of a throttle sequence. 109 * 110 * \param[in] epoch_i Inital epoch 111 * \param[in] state_i Initial sc_state (spacecraft state) 112 * \param[in] throttles_start iterator pointing to the beginning of a cartesian throttle sequence. 113 * \param[in] throttles_end iterator pointing to the end+1 of a cartesian throttle sequence. 114 * \param[in] epoch_f Final epoch. Needs to be later than epoch_i 115 * \param[in] state_f Final sc_state (spacecraft state) 116 * \param[in] mu_ Primary body gravitational constant 117 * 118 & \throws value_error if final epoch is before initial epoch, if mu_ not positive 119 */ 120 template <typename it_type> set_leg(const epoch & epoch_i,const sc_state & state_i,it_type throttles_start,it_type throttles_end,const epoch & epoch_f,const sc_state & state_f,double mu_,typename boost::enable_if<boost::is_same<typename std::iterator_traits<it_type>::value_type,throttle>>::type * =0)121 void set_leg(const epoch &epoch_i, const sc_state &state_i, it_type throttles_start, it_type throttles_end, 122 const epoch &epoch_f, const sc_state &state_f, double mu_, 123 typename boost::enable_if< 124 boost::is_same<typename std::iterator_traits<it_type>::value_type, throttle>>::type * = 0) 125 { 126 if (epoch_f.mjd2000() <= epoch_i.mjd2000()) { 127 throw_value_error("Final epoch is before initial epoch"); 128 } 129 130 t_i = epoch_i; 131 x_i = state_i; 132 t_f = epoch_f; 133 x_f = state_f; 134 135 throttles.assign(throttles_start, throttles_end); 136 137 if (mu_ <= 0) { 138 throw_value_error("Gravitational constant is less or equal to zero"); 139 } 140 m_mu = mu_; 141 } 142 143 /// Initialize a leg 144 /** 145 * Initialize a leg assuming that the user has or will initialize separately the spacecraft. 146 * The throttles are provided via two iterators pointing to the beginning and to the end of 147 * a sequence of doubles (\f$ x_1,y_1,z_1, ..., x_N,y_N,z_N \f$ containing the 148 * cartesian components of each throttle \f$ x_i,y_i,z_i \in [0,1]\f$. The constructed leg 149 * will have by default equally spaced segments. 150 * 151 * \param[in] epoch_i Inital epoch 152 * \param[in] state_i Initial sc_state (spacecraft state) 153 * \param[in] throttles_start iterator pointing to the beginning of a cartesian throttle sequence. 154 * \param[in] throttles_end iterator pointing to the end+1 of a cartesian throttle sequence. 155 * \param[in] epoch_f Final epoch. Needs to be later than epoch_i 156 * \param[in] state_f Final sc_state (spacecraft state) 157 * \param[in] mu_ Primary body gravitational constant 158 * 159 & \throws value_error if final epoch is before initial epoch, if mu_ not positive 160 */ 161 162 template <typename it_type> 163 void set_leg(const epoch & epoch_i,const sc_state & state_i,it_type throttles_start,it_type throttles_end,const epoch & epoch_f,const sc_state & state_f,double mu_,typename boost::enable_if<boost::is_same<typename std::iterator_traits<it_type>::value_type,double>>::type * =0)164 set_leg(const epoch &epoch_i, const sc_state &state_i, it_type throttles_start, it_type throttles_end, 165 const epoch &epoch_f, const sc_state &state_f, double mu_, 166 typename boost::enable_if<boost::is_same<typename std::iterator_traits<it_type>::value_type, double>>::type 167 * = 0) 168 { 169 // if (epoch_f.mjd2000() <= epoch_i.mjd2000()) 170 // { 171 // throw_value_error("Final epoch is before initial epoch"); 172 // } 173 if (std::distance(throttles_start, throttles_end) % 3 || std::distance(throttles_start, throttles_end) <= 0) { 174 throw_value_error("The length of the throttles list must be positive and a multiple of 3"); 175 } 176 if (mu_ <= 0) { 177 throw_value_error("Gravitational constant is less or equal to zero"); 178 } 179 m_mu = mu_; 180 181 t_i = epoch_i; 182 x_i = state_i; 183 t_f = epoch_f; 184 x_f = state_f; 185 186 auto throttles_vector_size = (throttles_end - throttles_start) / 3; 187 throttles.resize(throttles_vector_size); 188 const double seg_duration = (epoch_f.mjd() - epoch_i.mjd()) / static_cast<double>(throttles_vector_size); 189 for (decltype(throttles_vector_size) i = 0u; i < throttles_vector_size; ++i) { 190 kep_toolbox::array3D tmp 191 = {{*(throttles_start + 3 * i), *(throttles_start + 3 * i + 1), *(throttles_start + 3 * i + 2)}}; 192 throttles[i] = throttle(epoch(epoch_i.mjd() + seg_duration * static_cast<double>(i), epoch::MJD), 193 epoch(epoch_i.mjd() + seg_duration * (static_cast<double>(i) + 1), epoch::MJD), tmp); 194 } 195 } 196 197 /// Initialize a leg set_leg(const epoch & epoch_i,const sc_state & state_i,const std::vector<double> & thrott,const epoch & epoch_f,const sc_state & state_f)198 void set_leg(const epoch &epoch_i, const sc_state &state_i, const std::vector<double> &thrott, const epoch &epoch_f, 199 const sc_state &state_f) 200 { 201 set_leg(epoch_i, state_i, thrott.begin(), thrott.end(), epoch_f, state_f, m_mu); 202 } 203 204 /** @name Setters*/ 205 //@{ 206 207 /// Sets the leg's spacecraft 208 /** 209 * 210 *In order for the trajectory leg to be able to propagate the states, information on the 211 * low-thrust propulsion system used needs to be available. This is provided by the object 212 * spacecraft private member of the class and can be set using this setter. 213 * 214 * \param[in] sc The spacecraft object 215 */ set_spacecraft(const spacecraft & sc)216 void set_spacecraft(const spacecraft &sc) 217 { 218 m_sc = sc; 219 } 220 221 /// Sets the leg's primary body gravitational parameter 222 /** 223 * 224 * Sets the leg's central body gravitational parameter 225 * 226 * \param[in] mu_ The gravitational parameter 227 */ set_mu(const double & mu_)228 void set_mu(const double &mu_) 229 { 230 m_mu = mu_; 231 } 232 233 /// Sets the throttles 234 /** 235 * 236 * \param[in] b iterator pointing to the begin of a throttles sequence 237 * \param[in] e iterator pointing to the end of a throttles sequence 238 */ 239 template <typename it_type> set_throttles(it_type b,it_type e)240 void set_throttles(it_type b, it_type e) 241 { 242 throttles.assign(b, e); 243 } 244 245 /// Sets the throttles size 246 /** 247 * Resizes the throttles vector to a new size. 248 * 249 * \param[in] size The new size of the throttles vector (number of segments) 250 */ set_throttles_size(const int & size)251 void set_throttles_size(const int &size) 252 { 253 throttles.resize(size); 254 } 255 256 /** 257 * Sets the throttles 258 * 259 * \param[in] t std::vector of throttles 260 */ set_throttles(const std::vector<throttle> & t)261 void set_throttles(const std::vector<throttle> &t) 262 { 263 throttles = t; 264 } 265 266 /** 267 * Sets the ith throttle 268 * 269 * \param[in] index the index of the throttle 270 * \param[in] t the throttle 271 */ set_throttles(int index,const throttle & t)272 void set_throttles(int index, const throttle &t) 273 { 274 throttles[index] = t; 275 } 276 277 /// Sets the final sc_state 278 /** 279 * Sets the spacecraft state at the end of the leg 280 * 281 */ set_x_f(const sc_state & s)282 void set_x_f(const sc_state &s) 283 { 284 x_f = s; 285 } 286 287 /// Sets the initial sc_state 288 /** 289 * Sets the spacecraft state at the beginning of the leg 290 * 291 */ set_x_i(const sc_state & s)292 void set_x_i(const sc_state &s) 293 { 294 x_i = s; 295 } 296 //@} 297 298 /// Sets the initial epoch 299 /** 300 * Sets the epoch at the beginning of the leg 301 * 302 */ set_t_i(epoch e)303 void set_t_i(epoch e) 304 { 305 t_i = e; 306 } 307 308 /// Sets the final epoch 309 /** 310 * Sets the epoch at the end of the leg 311 * 312 */ set_t_f(epoch e)313 void set_t_f(epoch e) 314 { 315 t_f = e; 316 } 317 318 /// Sets the leg high fidelity state 319 /** 320 * Activates the evaluation of the state-mismatches using a high-fidelity model. Resulting leg 321 * is a real low-thrust trajectory (i.e. no impulses approximation) 322 * 323 */ set_high_fidelity(bool state)324 void set_high_fidelity(bool state) 325 { 326 m_hf = state; 327 } 328 329 /** @name Getters*/ 330 //@{ 331 332 /// Gets the leg's spacecraft 333 /** 334 * Returns the spacecraft 335 * 336 * @return sc const reference to spacecraft object 337 */ get_spacecraft() const338 const spacecraft &get_spacecraft() const 339 { 340 return m_sc; 341 } 342 343 /// Gets the gravitational parameter 344 /** 345 * @return the gravitational parameter 346 */ get_mu() const347 double get_mu() const 348 { 349 return m_mu; 350 } 351 352 /// Gets the throttle vector size 353 /** 354 * Returns the throttle vector size (number of segments) 355 * 356 * @return size_t containing the throttle vector size. 357 */ get_throttles_size() const358 size_t get_throttles_size() const 359 { 360 return throttles.size(); 361 } 362 363 /// Gets the i-th throttle 364 /** 365 * Returns the i-th throttle 366 * 367 * @return const ref to the i-th throttle 368 */ get_throttles(int index)369 const throttle &get_throttles(int index) 370 { 371 return throttles[index]; 372 } 373 374 /// Gets the throttles 375 /** 376 * Returns all throttles 377 * 378 * @return const ref to a vector of throttle 379 */ get_throttles()380 const std::vector<throttle> &get_throttles() 381 { 382 return throttles; 383 } 384 385 /// Gets the leg's initial epoch 386 /** 387 * Gets the epoch at the beginning of the leg 388 * 389 * @return const reference to the initial epoch 390 */ get_t_i() const391 const epoch &get_t_i() const 392 { 393 return t_i; 394 } 395 396 /// Gets the leg's final epoch 397 /** 398 * Gets the epoch at the end of the leg 399 * 400 * @return const reference to the final epoch 401 */ get_t_f() const402 const epoch &get_t_f() const 403 { 404 return t_f; 405 } 406 407 /// Gets the sc_state at the end of the leg 408 /** 409 * Gets the spacecraft state at the end of the leg 410 * 411 * @return const reference to the final sc_state 412 */ get_x_f() const413 const sc_state &get_x_f() const 414 { 415 return x_f; 416 } 417 418 /// Gets the initial sc_state 419 /** 420 * Gets the spacecraft state at the beginning of the leg 421 * 422 * @return const reference to the initial sc_state 423 */ get_x_i() const424 const sc_state &get_x_i() const 425 { 426 return x_i; 427 } get_high_fidelity() const428 bool get_high_fidelity() const 429 { 430 return m_hf; 431 } 432 //@} 433 434 /** @name Leg Feasibility*/ 435 //@{ 436 437 /// Evaluate the state mismatch 438 /** 439 * This is the main method of the class leg as it performs the orbital propagation from the initial sc_state, and 440 * accounting for all the throttles, up to a mid-point. The same is done starting from the final sc_state up to 441 * the same mid-point. The difference between the obtained values is then recorded at the memory location pointed by 442 * the iterators 443 * If not all zero the leg is unfeasible. The values stored are \f$\mathbf r, \mathbf v, m\f$ 444 * 445 * @param begin iterator pointing to the beginning of the memory where the mismatches will be stored 446 * @param begin iterator pointing to the end of the memory where the mismatches will be stored 447 */ 448 449 template <typename it_type> get_mismatch_con(it_type begin,it_type end) const450 void get_mismatch_con(it_type begin, it_type end) const 451 { 452 if (m_hf) { 453 get_mismatch_con_low_thrust(begin, end); 454 } else { 455 get_mismatch_con_chemical(begin, end); 456 } 457 } 458 459 protected: 460 template <typename it_type> get_mismatch_con_chemical(it_type begin,it_type end) const461 void get_mismatch_con_chemical(it_type begin, it_type end) const 462 { 463 assert(end - begin == 7); 464 (void)end; 465 size_t n_seg = throttles.size(); 466 auto n_seg_fwd = (n_seg + 1) / 2, n_seg_back = n_seg / 2; 467 468 // Aux variables 469 double max_thrust = m_sc.get_thrust(); 470 double isp = m_sc.get_isp(); 471 double norm_dv; 472 array3D dv; 473 474 // Initial state 475 array3D rfwd = x_i.get_position(); 476 array3D vfwd = x_i.get_velocity(); 477 double mfwd = x_i.get_mass(); 478 479 // Forward Propagation 480 double current_time_fwd = t_i.mjd2000() * ASTRO_DAY2SEC; 481 for (decltype(n_seg_fwd) i = 0u; i < n_seg_fwd; ++i) { 482 double thrust_duration 483 = (throttles[i].get_end().mjd2000() - throttles[i].get_start().mjd2000()) * ASTRO_DAY2SEC; 484 double manouver_time 485 = (throttles[i].get_start().mjd2000() + throttles[i].get_end().mjd2000()) / 2. * ASTRO_DAY2SEC; 486 propagate_lagrangian(rfwd, vfwd, manouver_time - current_time_fwd, m_mu); 487 current_time_fwd = manouver_time; 488 489 for (unsigned j = 0u; j < 3; j++) { 490 dv[j] = max_thrust / mfwd * thrust_duration * throttles[i].get_value()[j]; 491 } 492 493 norm_dv = norm(dv); 494 sum(vfwd, vfwd, dv); 495 mfwd *= exp(-norm_dv / isp / ASTRO_G0); 496 // Temporary solution to the creation of NaNs when mass gets too small (i.e. 0) 497 if (mfwd < 1) mfwd = 1; 498 } 499 500 // Final state 501 array3D rback = x_f.get_position(); 502 array3D vback = x_f.get_velocity(); 503 double mback = x_f.get_mass(); 504 505 // Backward Propagation 506 double current_time_back = t_f.mjd2000() * ASTRO_DAY2SEC; 507 for (decltype(n_seg_back) i = 0; i < n_seg_back; i++) { 508 double thrust_duration = (throttles[throttles.size() - i - 1].get_end().mjd2000() 509 - throttles[throttles.size() - i - 1].get_start().mjd2000()) 510 * ASTRO_DAY2SEC; 511 double manouver_time = (throttles[throttles.size() - i - 1].get_start().mjd2000() 512 + throttles[throttles.size() - i - 1].get_end().mjd2000()) 513 / 2. * ASTRO_DAY2SEC; 514 // manouver_time - current_time_back is negative, so this should propagate backwards 515 propagate_lagrangian(rback, vback, manouver_time - current_time_back, m_mu); 516 current_time_back = manouver_time; 517 518 for (int j = 0; j < 3; j++) { 519 dv[j] = -max_thrust / mback * thrust_duration * throttles[throttles.size() - i - 1].get_value()[j]; 520 } 521 norm_dv = norm(dv); 522 sum(vback, vback, dv); 523 mback *= exp(norm_dv / isp / ASTRO_G0); 524 } 525 526 // finally, we propagate from current_time_fwd to current_time_back with a keplerian motion 527 propagate_lagrangian(rfwd, vfwd, current_time_back - current_time_fwd, m_mu); 528 529 // Return the mismatch 530 diff(rfwd, rfwd, rback); 531 diff(vfwd, vfwd, vback); 532 533 std::copy(rfwd.begin(), rfwd.end(), begin); 534 std::copy(vfwd.begin(), vfwd.end(), begin + 3); 535 begin[6] = mfwd - mback; 536 } 537 538 template <typename it_type> get_mismatch_con_low_thrust(it_type begin,it_type end) const539 void get_mismatch_con_low_thrust(it_type begin, it_type end) const 540 { 541 assert(end - begin == 7); 542 (void)end; 543 auto n_seg = throttles.size(); 544 auto n_seg_fwd = (n_seg + 1) / 2, n_seg_back = n_seg / 2; 545 546 // Aux variables 547 double max_thrust = m_sc.get_thrust(); 548 double veff = m_sc.get_isp() * ASTRO_G0; 549 array3D thrust; 550 551 // Initial state 552 array3D rfwd = x_i.get_position(); 553 array3D vfwd = x_i.get_velocity(); 554 double mfwd = x_i.get_mass(); 555 556 // Forward Propagation 557 for (decltype(n_seg_fwd) i = 0; i < n_seg_fwd; ++i) { 558 double thrust_duration 559 = (throttles[i].get_end().mjd2000() - throttles[i].get_start().mjd2000()) * ASTRO_DAY2SEC; 560 561 for (unsigned j = 0u; j < 3; j++) { 562 thrust[j] = max_thrust * throttles[i].get_value()[j]; 563 } 564 propagate_taylor(rfwd, vfwd, mfwd, thrust, thrust_duration, m_mu, veff, m_tol, m_tol); 565 } 566 567 // Final state 568 array3D rback = x_f.get_position(); 569 array3D vback = x_f.get_velocity(); 570 double mback = x_f.get_mass(); 571 572 // Backward Propagation 573 for (decltype(n_seg_back) i = 0; i < n_seg_back; i++) { 574 double thrust_duration = (throttles[throttles.size() - i - 1].get_end().mjd2000() 575 - throttles[throttles.size() - i - 1].get_start().mjd2000()) 576 * ASTRO_DAY2SEC; 577 for (unsigned j = 0u; j < 3; j++) { 578 thrust[j] = max_thrust * throttles[throttles.size() - i - 1].get_value()[j]; 579 } 580 propagate_taylor(rback, vback, mback, thrust, -thrust_duration, m_mu, veff, m_tol, m_tol); 581 } 582 583 // Return the mismatch 584 diff(rfwd, rfwd, rback); 585 diff(vfwd, vfwd, vback); 586 587 std::copy(rfwd.begin(), rfwd.end(), begin); 588 std::copy(vfwd.begin(), vfwd.end(), begin + 3); 589 begin[6] = mfwd - mback; 590 } 591 592 public: 593 /// Evaluate the state mismatch 594 /** 595 * This method overloads the same method using iterators but the mismatch values are stored 596 * in a sc_state object 597 * 598 * @param retval the state mismatch structured as a spacecraft state 599 */ get_mismatch_con(sc_state & retval) const600 void get_mismatch_con(sc_state &retval) const 601 { 602 array7D tmp; 603 get_mismatch_con(tmp.begin(), tmp.end()); 604 retval.set_state(tmp); 605 } 606 607 /// Evaluate the throttles magnitude 608 /** 609 * This methods loops on the vector containing the throttles \f$ (x_1,y_1,z_1,x_2,y_2,z_2,...,x_n,y_n,z_n) \f$ 610 * and stores the magnitudes \f$ x_i^2 + y_i^2 + z_i^2 - 1\f$ at the locations indicated by the iterators. The 611 * iterators must have a distance of \f$ n\f$. If the stored values are not all \f$ \le 0 \f$ then the trajectory 612 * is unfeasible. 613 * 614 * @param[out] start std::vector<double>iterator from the first element where to store the magnitudes 615 * @param[out] start std::vector<double>iterator to the last+1 element where to store the magnitudes 616 */ 617 template <typename it_type> get_throttles_con(it_type start,it_type end) const618 void get_throttles_con(it_type start, it_type end) const 619 { 620 if ((end - start) != (int)throttles.size()) { 621 throw_value_error("Iterators distance is incompatible with the throttles size"); 622 } 623 int i = 0; 624 while (start != end) { 625 const array3D &t = throttles[i].get_value(); 626 *start = std::inner_product(t.begin(), t.end(), t.begin(), -1.); 627 ++i; 628 ++start; 629 } 630 } 631 //@} 632 633 /// Approximate the leg dv 634 /** 635 * This method returns the leg dv assuming a constant mass throughout the leg. Useful when 636 * mass propagation is switched off by setting an infinite specific impulse 637 * 638 * @return the leg dv (mass is considered constant) 639 */ evaluate_dv() const640 double evaluate_dv() const 641 { 642 double tmp = 0; 643 for (std::vector<double>::size_type i = 0; i < throttles.size(); ++i) { 644 tmp += (throttles[i].get_end().mjd2000() - throttles[i].get_start().mjd2000()) * ASTRO_DAY2SEC 645 * throttles[i].get_norm() * m_sc.get_thrust() / m_sc.get_mass(); 646 } 647 return tmp; 648 } 649 650 private: 651 // Serialization code 652 friend class boost::serialization::access; 653 template <class Archive> serialize(Archive & ar,const unsigned int)654 void serialize(Archive &ar, const unsigned int) 655 { 656 ar &t_i; 657 ar &x_i; 658 ar &throttles; 659 ar &t_f; 660 ar &x_f; 661 ar &m_sc; 662 ar &m_mu; 663 ar &m_hf; 664 ar &m_tol; 665 } 666 // Serialization code (END) 667 epoch t_i; 668 sc_state x_i; 669 std::vector<throttle> throttles; 670 epoch t_f; 671 sc_state x_f; 672 spacecraft m_sc; 673 double m_mu; 674 bool m_hf; 675 int m_tol; 676 }; 677 678 KEP_TOOLBOX_DLL_PUBLIC std::ostream &operator<<(std::ostream &s, const leg &in); 679 } // namespace sims_flanagan 680 } // namespace kep_toolbox 681 #endif // KEP_TOOLBOX_LEG_H 682