end of work on the mga problem .. notebook added

doc/conf.py

@@ -74,6 +74,7 @@
 nb_execution_excludepatterns = [
  "udp_point2point*",
  "udp_pl2pl*",
+ "/udp_mga*" 
 ]
 
 latex_engine = "xelatex"

doc/tut_trajopt.rst

@@ -28,6 +28,7 @@ Trajectory Optimization
 
  notebooks/udp_point2point
  notebooks/udp_pl2pl
+ notebooks/udp_mga
 
 
 

pykep/plot/__init__.py

@@ -41,7 +41,7 @@ def propagate_lagrangian_theta_v(rv=[[1, 0, 0], [0, 1, 0]], thetas=[0.0, _np.pi
 
  Propagates (Keplerian) the state for an assigned difference in True anomalies. It does not compute the State Transition Matrix
  A similar API is offered as that of :func:`~pykep.propagate_lagrangian`, but not identical. The function is essentially offered only for plotting
- purposes as to avoid (or anyway alleviate) the non uniform grid effect deriving when plotting at equal time intervals.
+ purposes as to avoid (or anyway alleviate) the non uniform grid effect when plotting at equal time intervals.
 
  Args:
  *rv* (2D array-like): Cartesian components of the initial position vector and velocity [[x0, y0, z0], [v0, vy0, vz0]]. Defaults to [[1,0,0], [0,1,0]].
@@ -75,9 +75,9 @@ def propagate_lagrangian_theta_v(rv=[[1, 0, 0], [0, 1, 0]], thetas=[0.0, _np.pi
  v0 = _np.array(rv[1])
  R0 = _np.linalg.norm(r0)
  V02 = _np.dot(v0, v0)
- energy = V02 / 2 - mu / R0
+ #energy = V02 / 2 - mu / R0
  # energy will be negative for hyperbolae
- a = -mu / 2.0 / energy
+ #a = -mu / 2.0 / energy
  sigma0 = _np.dot(r0, v0) / _np.sqrt(mu)
  h = _np.linalg.norm(_np.cross(r0, v0))
  p = h * h / mu

pykep/plot/_lambert.py

@@ -40,15 +40,30 @@ def add_lambert(ax, lp, N: int = 60, sol: int = 0, units=_pk.AU, **kwargs):
  v0 = lp.v0[sol]
  r1 = lp.r1
  mu = lp.mu
- costheta = _np.dot(r0, r1) / _np.linalg.norm(r0) / _np.linalg.norm(r1)
- theta = _np.arccos(costheta)
-
+
+# Compute the magnitudes of the vectors
+ norm_r1 = _np.linalg.norm(r0)
+ norm_r2 = _np.linalg.norm(r1)
+
+ # Compute the dot product
+ dot_product = _np.dot(r0, r1)
+
+ # Compute the angle using the dot product (gives the cosine of the angle)
+ cos_theta = dot_product / (norm_r1 * norm_r2)
+
+ # Compute the angle in radians (clipped to avoid domain errors due to floating point precision)
+ theta = _np.arccos(_np.clip(cos_theta, -1.0, 1.0))
+
+ # Compute the cross product to determine the direction
+ cross_product = _np.cross(r0, r1)
+
+ # Assuming motion is in the xy-plane (i.e., third component z gives the direction)
+ if cross_product[2] < 0: # If the z-component is negative, the motion is clockwise
+ theta = 2 * _np.pi - theta # Adjust the angle to account for counterclockwise motion
+
  # We define the integration grid
  if sol == 0:
- if _np.cross(r0, r1)[0] > 0:
- thetagrid = _np.linspace(0, theta, N)
- else:
- thetagrid = _np.linspace(0, 2 * _np.pi - theta, N)
+ thetagrid = _np.linspace(0, theta, N)
  else:
  thetagrid = _np.linspace(0, 2 * _np.pi, N)
 

pykep/trajopt/_mga.py

@@ -303,7 +303,7 @@ def _compute_dvs(self, x: List[float]) -> Tuple[
  r[i], v[i] = self.seq[i].eph(ep[i])
 
  # 3 - we solve the lambert problems (and store trajectory r,v)
- l = list()
+ lps = list()
  for i in range(self._n_legs):
  lp = _pk.lambert_problem(
  r0=r[i],
@@ -313,20 +313,20 @@ def _compute_dvs(self, x: List[float]) -> Tuple[
  cw=False,
  multi_revs=0,
  )
- l.append(lp)
+ lps.append(lp)
 
  # 4 - we compute the various dVs needed at fly-bys to match incoming
  # and outcoming
  DVfb = list()
- for i in range(len(l) - 1):
- v_rel_in = [a - b for a, b in zip(l[i].v1[0], v[i + 1])]
- v_rel_out = [a - b for a, b in zip(l[i + 1].v0[0], v[i + 1])]
+ for i in range(len(lps) - 1):
+ v_rel_in = [a - b for a, b in zip(lps[i].v1[0], v[i + 1])]
+ v_rel_out = [a - b for a, b in zip(lps[i + 1].v0[0], v[i + 1])]
  DVfb.append(_pk.fb_dv(v_rel_in, v_rel_out, self.seq[i + 1]))
 
  # 5 - we add the departure and arrival dVs
- DVlaunch_tot = _np.linalg.norm([a - b for a, b in zip(v[0], l[0].v0[0])])
+ DVlaunch_tot = _np.linalg.norm([a - b for a, b in zip(v[0], lps[0].v0[0])])
  DVlaunch = max(0, DVlaunch_tot - self.vinf)
- DVarrival = _np.linalg.norm([a - b for a, b in zip(v[-1], l[-1].v1[0])])
+ DVarrival = _np.linalg.norm([a - b for a, b in zip(v[-1], lps[-1].v1[0])])
  if self.orbit_insertion:
  # In this case we compute the insertion DV as a single pericenter
  # burn
@@ -339,7 +339,7 @@ def _compute_dvs(self, x: List[float]) -> Tuple[
  - MU_SELF / self.rp_target * (1.0 - self.e_target)
  )
  DVarrival = _np.abs(DVper - DVper2)
- return (DVlaunch, DVfb, DVarrival, l, DVlaunch_tot, ep, T)
+ return (DVlaunch, DVfb, DVarrival, lps, DVlaunch_tot, ep, T)
 
  # Objective function
  def fitness(self, x):
@@ -459,7 +459,7 @@ def pretty(self, x):
 
  print("\nTotal DV: ", DVlaunch + _np.sum(DVfb) + DVarrival)
 
- def plot(self, x, ax=None, units=_pk.AU, N=60, figsize=(5, 5)):
+ def plot(self, x, ax=None, units=_pk.AU, N=60, c_orbit = 'dimgray', c = 'indianred', leg_ids = [], figsize=(5, 5), **kwargs):
  """
  Plots the trajectory leg 3D axes.
 
@@ -473,6 +473,10 @@ def plot(self, x, ax=None, units=_pk.AU, N=60, figsize=(5, 5)):
  *N* (:class:`int`, optional): The number of points to use when plotting the trajectory. Defaults to 60.
 
  *figsize* (:class:`tuple`): The figure size (only used if a*ax* is None and axis ave to be created.), Defaults to (5, 5).
+ 
+ *leg_ids* (:class:`list`): selects the legs to plot. Optional, defaults to all legs.
+ 
+ *\*\*kwargs*: Additional keyword arguments to pass to the trajectory plot (all Lambert arcs)
 
  Returns:
  :class:`mpl_toolkits.mplot3d.axes3d.Axes3D`: The 3D axis where the trajectory was plotted.
@@ -481,21 +485,20 @@ def plot(self, x, ax=None, units=_pk.AU, N=60, figsize=(5, 5)):
 
  if ax is None:
  ax = _pk.plot.make_3Daxis(figsize=figsize)
- _, _, _, _, _, mjd2000s, _ = self._compute_dvs(x)
+ _, _, _, lps, _, mjd2000s, _ = self._compute_dvs(x)
  for i, item in enumerate(self.seq):
- _pk.plot.add_planet(pla=item, ax=ax, when=mjd2000s[i], c="r", units=units)
- _pk.plot.add_planet_orbit(pla=item, ax=ax, units=units, N=N)
+ _pk.plot.add_planet(pla=item, ax=ax, when=mjd2000s[i], c=c, units=units)
+ _pk.plot.add_planet_orbit(pla=item, ax=ax, units=units, N=N, c=c_orbit)
 
  _pk.plot.add_sun(ax=ax)
+
+ if len(leg_ids) == 0:
+ for lp in lps:
+ _pk.plot.add_lambert(ax, lp, N = 60, sol = 0, units=units, c=c, **kwargs)
+ else:
+ for leg_id in leg_ids:
+ _pk.plot.add_lambert(ax, lps[leg_id], N = 60, sol = 0, units=units, c=c, **kwargs)
 
- _pk.plot.add_planet_orbit(
- pla=self.to_planet(x),
- ax=ax,
- plot_range=[mjd2000s[0], mjd2000s[-1]],
- c="r",
- units=units,
- N=N,
- )
  return ax