AstroToolbox/mga.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 <iomanip>
  12 #include <fstream>
  13 #include <math.h>
  14 #include <vector>
  15 #include <iostream>
  16 #include "Pl_Eph_An.h"
  17 #include "mga.h"
  18 #include "Lambert.h"
  19 #include "PowSwingByInv.h"
  20 #include "Astro_Functions.h"
  21 #define MAX(a, b) (a > b ? a : b)
  22
  23 using namespace std;
  24
  25 //the function return 0 if the input is right or -1 it there is something wrong
  26
  27 int MGA(vector<double> t,       // it is the vector which provides time in modified julian date 2000.
  28                                                                 // The first entry is launch date, the next entries represent the time needed to
  29                                                                 // fly from last swing-by to current swing-by.
  30                         mgaproblem problem,
  31
  32                         /* OUTPUT values: */
  33                         vector<double>& rp,  // periplanets radius
  34                         vector<double>& DV,  // final delta-Vs
  35                         double &obj_funct)   //objective function
  36
  37 {
  38         const int n = problem.sequence.size();
  39         const vector<int> sequence = problem.sequence;
  40         const vector<int> rev_flag = problem.rev_flag;// array containing 0 clockwise, 1 un-clockwise
  41         customobject cust_obj = problem.asteroid;
  42
  43         double MU[9] = {//1.32712440018e11, //SUN = 0
  44                                         1.32712428e11,
  45                                         22321,          // Gravitational constant of Mercury    = 1
  46                                         324860,         // Gravitational constant of Venus              = 2
  47                                         398601.19,      // Gravitational constant of Earth              = 3
  48                                         42828.3,        // Gravitational constant of Mars               = 4
  49                                         126.7e6,        // Gravitational constant of Jupiter    = 5
  50                                         37.9e6,         // Gravitational constant of Saturn             = 6
  51                                         5.78e6,         // Gravitational constant of Uranus             = 7
  52                                         6.8e6           // Gravitational constant of Neptune    = 8
  53                                     };
  54         double penalty[9] = {0,
  55                                 0,        // Mercury
  56                                                 6351.8,   // Venus
  57                                                 6778.1,   // Earth
  58                                                 6000,     // Mars
  59                                                 //671492, // Jupiter
  60                         600000,   // Jupiter
  61                                                 70000,    // Saturn
  62                         0,        // Uranus
  63                                                 0         // Neptune
  64         };
  65
  66         double penalty_coeffs[9] = {0,
  67                                                                 0,      // Mercury
  68                                                                 0.01,   // Venus
  69                                                                 0.01,   // Earth
  70                                                                 0.01,   // Mars
  71                                                                 0.001,  // Jupiter
  72                                                                 0.01,   // Saturn
  73                                                                 0,      // Uranus
  74                                                                 0       // Neptune
  75         };
  76
  77         double DVtot = 0;
  78         double Dum_Vec[3],Vin,Vout;
  79         double V_Lamb[2][2][3],dot_prod;
  80         double a,p,theta,alfa;
  81         double DVrel, DVarr=0;
  82
  83         //only used for orbit insertion (ex: cassini)
  84                 double DVper, DVper2;
  85                 const double rp_target = problem.rp;
  86                 const double e_target = problem.e;
  87                 const double DVlaunch = problem.DVlaunch;
  88
  89         //only used for asteroid impact (ex: gtoc1)
  90                 const double initial_mass = problem.mass;   // Satellite initial mass [Kg]
  91                 double final_mass;                                                      // satelite final mass
  92                 const double Isp = problem.Isp;             // Satellite specific impulse [s]
  93                 const double g = 9.80665 / 1000.0;          // Gravity
  94
  95
  96
  97
  98         double *vec, *rec;
  99         vector<double*> r;    // {0...n-1} position
 100         vector<double*> v;    // {0...n-1} velocity
 101
 102         double T = 0.0;         // total time
 103
 104         int i_count, j_count, lw;
 105
 106         int iter = 0;
 107
 108         if (n >= 2)
 109         {
 110                 for ( i_count = 0; i_count < n; i_count++)
 111                 {
 112                         vec = new double [3];  // velocity and position are 3 D vector
 113                         rec = new double [3];
 114                         r.push_back(vec);
 115                         v.push_back(rec);
 116
 117                         DV [i_count] = 0.0;
 118                 }
 119
 120                 T = 0;
 121                 for (i_count = 0; i_count < n; i_count++)
 122                 {
 123                         T += t[i_count];
 124                         if (sequence[i_count]<10)
 125                                 Planet_Ephemerides_Analytical (T, sequence[i_count],
 126                                         r[i_count], v[i_count]); //r and  v in heliocentric coordinate system
 127                         else
 128                         {
 129                                 Custom_Eph(T+2451544.5, cust_obj.epoch, cust_obj.keplerian, r[i_count], v[i_count]);
 130                         }
 131                 }
 132
 133                 vett(r[0], r[1], Dum_Vec);
 134
 135                 if (Dum_Vec[2] > 0)
 136                         lw = (rev_flag[0] == 0) ? 0 : 1;
 137             else
 138                         lw = (rev_flag[0] == 0) ? 1 : 0;
 139
 140                 LambertI(r[0],r[1],t[1]*24*60*60,MU[0],lw,          // INPUT
 141                              V_Lamb[0][0],V_Lamb[0][1],a,p,theta,iter); // OUTPUT
 142                 DV[0] = norm(V_Lamb[0][0], v[0]);                // Earth launch
 143
 144                 for (i_count = 1; i_count <= n-2; i_count++)
 145                 {
 146                         vett(r[i_count], r[i_count+1], Dum_Vec);
 147
 148                         if (Dum_Vec[2] > 0)
 149                                 lw = (rev_flag[i_count] == 0) ? 0 : 1;
 150                         else
 151                                 lw = (rev_flag[i_count] == 0) ? 1 : 0;
 152
 153                         /*if (i_count%2 != 0)   {*/
 154                         LambertI(r[i_count],r[i_count+1],t[i_count + 1]*24*60*60,MU[0],lw, // INPUT
 155                                    V_Lamb[1][0],V_Lamb[1][1],a,p,theta,iter);                  // OUTPUT
 156
 157                         // norm first perform the subtraction of vet1-vet2 and the evaluate ||...||
 158                         Vin  = norm(V_Lamb[0][1], v[i_count]);
 159                         Vout = norm(V_Lamb[1][0], v[i_count]);
 160
 161                         dot_prod = 0.0;
 162                         for (int i = 0; i < 3; i++)
 163                         {
 164                                 dot_prod += (V_Lamb[0][1][i] - v[i_count][i]) * (V_Lamb[1][0][i] - v[i_count][i]);
 165                         }
 166                         alfa = acos ( dot_prod /(Vin * Vout) );
 167
 168                         // calculation of delta V at pericenter
 169                         PowSwingByInv(Vin, Vout, alfa, DV[i_count], rp[i_count - 1]);
 170
 171                         rp[i_count - 1] *= MU[sequence[i_count]];
 172
 173                         if (i_count != n-2)  //swap
 174                         for (j_count = 0; j_count < 3; j_count++)
 175                                 {
 176                                         V_Lamb[0][0][j_count] = V_Lamb[1][0][j_count];  // [j_count];
 177                                         V_Lamb[0][1][j_count] = V_Lamb[1][1][j_count];  // [j_count];
 178                                 }
 179                         }
 180         }
 181         else
 182         {
 183                 return -1;
 184         }
 185
 186         for (i_count = 0; i_count < 3; i_count++)
 187                 Dum_Vec[i_count] = v[n-1][i_count] - V_Lamb[1][1][i_count];
 188
 189         DVrel = norm2(Dum_Vec);
 190
 191         if (problem.type == total_DV_orbit_insertion){
 192
 193                 DVper  = sqrt(DVrel*DVrel + 2*MU[sequence[n-1]]/rp_target);
 194                 DVper2 = sqrt(2*MU[sequence[n-1]]/rp_target - MU[sequence[n-1]]/rp_target*(1-e_target));
 195                 DVarr = fabs(DVper - DVper2);
 196         }
 197
 198         else if (problem.type == asteroid_impact){
 199
 200                 DVarr = DVrel;
 201         }
 202
 203
 204
 205         DVtot = 0;
 206
 207         for (i_count = 1; i_count < n-1; i_count++)
 208                 DVtot += DV[i_count];
 209
 210         if (problem.type == total_DV_orbit_insertion){
 211
 212                 DVtot += DVarr;
 213         }
 214
 215
 216
 217         // Build Penalty
 218         for (i_count = 0;i_count < n-2; i_count++)
 219                 if (rp[i_count] < penalty[sequence[i_count+1]])
 220                         DVtot += penalty_coeffs[sequence[i_count+1]]*fabs(rp[i_count] - penalty[sequence[i_count+1]]);
 221
 222         // Launcher Constraint
 223         if (DV[0] > DVlaunch)
 224                 DVtot += (DV[0] - DVlaunch);
 225
 226         if (problem.type == total_DV_orbit_insertion){
 227
 228                 obj_funct = DVtot;
 229         }
 230
 231         else if (problem.type == asteroid_impact){
 232
 233                 // Evaluation of satellite final mass
 234                 obj_funct = final_mass = initial_mass * exp(- DVtot/ (Isp * g));
 235
 236                 // V asteroid - V satellite
 237                 for (i_count = 0; i_count < 3; i_count++)
 238                         Dum_Vec[i_count] = v[n-1][i_count] - V_Lamb[1][1][i_count];// arrival relative velocity at the asteroid;
 239
 240                 dot_prod = 0;
 241                 for (i_count = 0; i_count < 3 ; i_count++)
 242                         dot_prod += Dum_Vec[i_count] * v[n-1][i_count];
 243
 244                 obj_funct = 2000000 - (final_mass)* fabs(dot_prod);
 245         }
 246
 247
 248         // final clean
 249         for ( i_count = 0;i_count < n;i_count++)
 250         {
 251                 delete [] r[i_count];
 252                 delete [] v[i_count];
 253         }
 254         r.clear();
 255         v.clear();
 256         return 0;
 257 }
 258
 259
 260