AstroToolbox/PowSwingByInv.cpp

   1 // ------------------------------------------------------------------------ //
   2 // This source file is part of the 'ESA Advanced Concepts Team's                        //
   3 // Space Mechanics Toolbox' software.                                       //
   4 //                                                                          //
   5 // The source files are for research use only,                              //
   6 // and are distributed WITHOUT ANY WARRANTY. Use them on your own risk.     //
   7 //                                                                          //
   8 // Copyright (c) 2004-2007 European Space Agency                            //
   9 // ------------------------------------------------------------------------ //
  10
  11 #include "PowSwingByInv.h"
  12 #include <math.h>
  13
  14 void PowSwingByInv (const double Vin,const double Vout,const double alpha,
  15                                         double &DV,double &rp)
  16 {
  17         const int maxiter = 30;
  18         int i = 0;
  19         double err = 1.0;
  20         double f,df;                    // function and its derivative
  21         double rp_new;
  22         const double tolerance = 1e-8;
  23
  24         double aIN  = 1.0/pow(Vin,2);     // semimajor axis of the incoming hyperbola
  25         double aOUT = 1.0/pow(Vout,2);    // semimajor axis of the incoming hyperbola
  26
  27         rp = 1.0;
  28         while ((err > tolerance)&&(i < maxiter))
  29         {
  30                 i++;
  31                 f = asin(aIN/(aIN + rp)) + asin(aOUT/(aOUT + rp)) - alpha;
  32                 df = -aIN/sqrt((rp + 2 * aIN) * rp)/(aIN+rp) -
  33                          aOUT/sqrt((rp + 2 * aOUT) * rp)/(aOUT+rp);
  34                 rp_new = rp - f/df;
  35                 if (rp_new > 0)
  36                 {
  37                         err = fabs(rp_new - rp);
  38                         rp = rp_new;
  39                 }
  40                 else
  41                         rp /= 2.0;
  42         }
  43
  44         // Evaluation of the DV
  45         DV = fabs (sqrt(Vout*Vout + (2.0/rp)) - sqrt(Vin*Vin + (2.0/rp)));
  46 }