PaGMO  1.1.5
PowSwingByInv.cpp
1 /*****************************************************************************
2  * Copyright (C) 2004-2015 The PaGMO development team, *
3  * Advanced Concepts Team (ACT), European Space Agency (ESA) *
4  * *
5  * https://github.com/esa/pagmo *
6  * *
7  * act@esa.int *
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  * This program is distributed in the hope that it will be useful, *
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
17  * GNU General Public License for more details. *
18  * *
19  * You should have received a copy of the GNU General Public License *
20  * along with this program; if not, write to the *
21  * Free Software Foundation, Inc., *
22  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
23  *****************************************************************************/
24 
25 #include "PowSwingByInv.h"
26 #include <math.h>
27 
28 void PowSwingByInv (const double Vin,const double Vout,const double alpha,
29  double &DV,double &rp)
30 {
31  const int maxiter = 30;
32  int i = 0;
33  double err = 1.0;
34  double f,df; // function and its derivative
35  double rp_new;
36  const double tolerance = 1e-8;
37 
38  double aIN = 1.0/pow(Vin,2); // semimajor axis of the incoming hyperbola
39  double aOUT = 1.0/pow(Vout,2); // semimajor axis of the incoming hyperbola
40 
41  rp = 1.0;
42  while ((err > tolerance)&&(i < maxiter))
43  {
44  i++;
45  f = asin(aIN/(aIN + rp)) + asin(aOUT/(aOUT + rp)) - alpha;
46  df = -aIN/sqrt((rp + 2 * aIN) * rp)/(aIN+rp) -
47  aOUT/sqrt((rp + 2 * aOUT) * rp)/(aOUT+rp);
48  rp_new = rp - f/df;
49  if (rp_new > 0)
50  {
51  err = fabs(rp_new - rp);
52  rp = rp_new;
53  }
54  else
55  rp /= 2.0;
56  }
57 
58  // Evaluation of the DV
59  DV = fabs (sqrt(Vout*Vout + (2.0/rp)) - sqrt(Vin*Vin + (2.0/rp)));
60 }