Skip to content
Snippets Groups Projects
Commit 34a1050e authored by Yvan's avatar Yvan
Browse files

Try to make coherence between updates

parent 8f805d0c
No related branches found
No related tags found
No related merge requests found
Showing
with 116 additions and 142 deletions
......@@ -28,4 +28,6 @@ recursive-include sempy/core/init/kernels/spk *
# Get abacuse
recursive-include sempy/core/orbits/abacuses *
# Get data from orekit
recursive-include sempy/core/utils/data *
......@@ -11,50 +11,30 @@ Created on Mon Jul 22 15:09:46 2019
import numpy as np
def generate_plane_event(plane, direction, terminal=True):
"""This function generates an event to detect plane crossings.
Parameters
----------
plane : str
coordinating plane for which the event should be created.
Must be one of ``xy``, ``xz``, ``yz``.
direction : int, optional
event direction. Must be one of ``1``, ``-1``, ``0`` where:
* | 1 -- detect events going from negative to positive w.r.t. the third coordinating axis
* | -1 -- detect events going from positive to negative w.r.t. the third coordinating axis
* | 0 -- detect all plane crossing events independently from their direction.
Default to ``0``.
terminal : bool, optional
if ``True``, propagation is stopped at the first event occurrence. If ``False``, no action
is taken and the propagation continues until the target time. Default to ``False``.
Returns
-------
fcn : function
an event function to be passed as keyword arguments to the various propagators.
def xz_plane_event(_, state):
"""Type of event: reaching the y = 0 plane."""
return state[1]
"""
if plane == 'xz':
coord = 1
elif plane == 'xy':
coord = 2
elif plane == 'yz':
coord = 0
if direction == 'neg_pos':
direct = +1 # Direction from negative to positive
elif direction == 'pos_neg':
direct = -1 # Direction from positive to positive
else:
raise Exception("Direction of event must be either 'exit' or 'enter'")
xz_plane_event.terminal = True
xz_plane_event.direction = -1
def xz_plane_event_positive(_, state):
"""Type of event: reaching the y = 0 plane, with a postive direction of the velociti alog y."""
return state[1]
xz_plane_event_positive.terminal = True
xz_plane_event_positive.direction = 1
def xz_plane_event_neutral(_, state):
"""Type of event: reaching the y = 0 plane, the direction does not matter."""
return state[1]
xz_plane_event_neutral.terminal = False
fun = lambda t, state: state[coord]
fun.direction = direct
fun.type = plane + ' plane'
fun.terminal = terminal # Do not stop integration if encountered
return fun
def generate_soi_event(cr3bp, direction, terminal=False):
......@@ -93,63 +73,55 @@ def generate_soi_event(cr3bp, direction, terminal=False):
return fun
def generate_impact_event(cr3bp, prim, alt_lim=None):
def generate_impact_event(cr3bp, prim):
""" This function generates the event functions to be passed to the solve_ivp integrator. It
is necessary to add m1_impact_event = generate_impact_event(cr3bp, "m1") in the main
script, and pass m1_impact_event to solve_ivp. Depending on the value of prim, this function
generates an impact event with m1 or m2. The alt_lim value allows to define an atmospheric
limit. For example, for the Earth it is often set at 120 km. """
generates an impact event with m1 or m2. """
def primary_impact(state, prim_pos, prim_rad):
prim_distance = np.linalg.norm(state[0:3] - prim_pos) - prim_rad
return prim_distance
if not alt_lim:
alt_lim = 0
if prim == "m1":
# Type of event: impact with M1. Terminal.
fun = lambda t, state: primary_impact(state,
cr3bp.m1_pos,
(cr3bp.R1+alt_lim) / cr3bp.L)
fun = lambda t, state: primary_impact(state, cr3bp.m1_pos, cr3bp.R1 / cr3bp.L)
fun.type = 'm1 impact'
elif prim == "m2":
# Type of event: impact with M2. Terminal.
fun = lambda t, state: primary_impact(state,
cr3bp.m2_pos,
(cr3bp.R2+alt_lim) / cr3bp.L)
fun = lambda t, state: primary_impact(state, cr3bp.m2_pos, cr3bp.R2 / cr3bp.L)
fun.type = 'm2 impact'
else:
raise Exception("Primary has to be 'm1' or 'm2'")
fun.type = prim + ' impact'
fun.terminal = True # Stops integration if met
fun.direction = -1 # Direction from positive to negative
return fun
def generate_x_event(cr3bp, x_value, direction, terminal=False):
def generate_li_gate_event(cr3bp, l_i, direction, terminal=False):
""" This function generates the event functions to be passed to the solve_ivp integrator. It
is necessary to add L1_exit_event = generate_x_event(cr3bp, cr3bp.l1.position[0], "exit") in the
main script, and pass L1_exit_event to solve_ivp. Direction can be either 'exit' or 'enter',
with reference to the x axis.
is necessary to add L1_exit_event = generate_li_gate_event(cr3bp, cr3bp.l1, "exit") in the
main script, and pass L1_exit_event to solve_ivp. Depending on the value of cr3bp.l1,
this function generates an exit event from the L1 or L2 gates (L3, L4 and L5 are not
considered). Direction can be either 'exit' or 'enter'.
"""
def x_crossing(state, _x):
return state[0] - _x
def gate_crossing(state, li_x):
return state[0] - li_x
if direction == 'exit':
direct = +1 # From negative to positive if exit
direct = +1 # From negative (x<L1 or x<L2) to positive (x>L1 or x>L2) if enter
elif direction == 'enter': # The opposite as exit
direct = -1 # From positive to negative if enter
direct = -1 # From positive (x>L1 or x>L2) to negative (x<L1 or x<L2) if enter
else:
raise Exception("Direction of event must be either 'exit' or 'enter'")
# Type of event: x crossing.
fun = lambda t, state: x_crossing(state, x_value)
if l_i not in (cr3bp.l1, cr3bp.l2):
raise Exception("Libration point has to be cr3bp.l1 or cr3bp.l2")
# Type of event: Li gate crossing.
fun = lambda t, state: gate_crossing(state, l_i.position[0])
fun.direction = direct
fun.terminal = terminal # Stops integration if met
if x_value == cr3bp.l1.position[0]:
fun.type = "L1 " + direction
elif x_value == cr3bp.l2.position[0]:
fun.type = "L2 " + direction
else:
fun.type = "x " + direction
fun.type = f"l{l_i.number} gate " + direction
return fun
No preview for this file type
No preview for this file type
No preview for this file type
File added
......@@ -27,9 +27,9 @@ class TestODEEvent(unittest.TestCase):
cr3bp = Cr3bp(Primary.EARTH, Primary.MOON)
# Events initialisation
xz_plane_event = generate_plane_event('xz', -1, True)
xy_plane_event = generate_plane_event('xy', -1, True)
yz_plane_event = generate_plane_event('yz', -1, True)
xz_plane_event = generate_plane_event('xz', 'pos_neg', True)
xy_plane_event = generate_plane_event('xy', 'pos_neg', True)
yz_plane_event = generate_plane_event('yz', 'pos_neg', True)
# Propagate initial state:
prop = Cr3bpSynodicPropagator(cr3bp.mu)
......
File added
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment