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