-
-
Notifications
You must be signed in to change notification settings - Fork 276
Add CR3BP functions created by Dhruv Jain #1475
base: main
Are you sure you want to change the base?
Changes from all commits
04c2d69
7295787
49cd924
6a74aae
ac3852f
f077f39
6cb460b
5485b8e
1305e56
48e6240
985bbe8
fb21132
b825dc0
3022263
2861210
83aaa65
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,27 @@ | ||
# Dhruv Jain - Circular Restircted Three Body Problem | ||
|
||
I hope to add the capability to investiage the dynamical flows in CR3BP for Poliastro users. The current files are capable of: | ||
|
||
1. cr3bp_char_quant.py: Computing the characterisitic quantities of a system in CR3BP: mu,l*, t* | ||
2. cr3bp_lib_JC_calc.py: Computing the 5 libration points position and Jacobi Constant | ||
3. cr3bp_master.py: Numerical Integration of CR3BP EOMs (with the option of using events functions) | ||
4. cr3bp_master.py: Numerical Integration of CR3BP + State Transition Matrix EOMs (with the option of using events functions) | ||
5. cr3bp_master.py: Computing the first derivative of pseudo potenial, second-derivative of pseudo potenial and acceleration terms | ||
|
||
I will also be working to leverage the preexisisting functions of poliastro to further streamline these features and tweak cr3bp_master.py to be a class. | ||
|
||
The current work requires the widely used numpy and matplotlib libraries. | ||
|
||
## Future Work: | ||
I want to make the entire setup more robust. Create a good starting point for users to compute families of periodic orbits and understand the general dynamical flows in CR3BP | ||
1. General purpose Multiple shooter to compute collinear libration point periodic orbit families (Planar and Spatial) | ||
2. Numerical Continuation schemes: Natural parameter continuation, Pseudo-arc length continuation | ||
3. Stability analysis: Sorting Eigenevalues and Eigenvectors of Monodromy matrix, Stability Index plots, Broucke Stability Diagram | ||
4. Manifold generation | ||
5. Heteroclinic transfer design using Tau-Alpha method | ||
6. Method to transition and corret in N-body ephemeris model | ||
|
||
My dream is to build enough features to enable me to include a targeter that I have created to compute Quasi-Periodic Orbits(QPOs) in CR3BP, so that users can access the less investaged but more useful QPOs. | ||
|
||
## About Me: | ||
I am Dhruv Jain, [email protected], pursuing Master of Science in Aeroanautics and Astronautics Engienering at Purdue University. I am part of the Multi-Body Dynamics Research Group and I am working on designing Quasi-Periodic Orbits in CR3BP. When I am not working, I like playing Badminton and cooking! |
Original file line number | Diff line number | Diff line change | ||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
@@ -0,0 +1,123 @@ | ||||||||||||||||||
""" | ||||||||||||||||||
Created on Mon Feb 21 20:37:16 2022 | ||||||||||||||||||
|
||||||||||||||||||
@author: Dhruv Jain, Multi-Body Dynamics Research Group, Purdue University | ||||||||||||||||||
[email protected] | ||||||||||||||||||
|
||||||||||||||||||
Objective: This file contains functions required to target a Periodic Orbit in the Circular Restricted | ||||||||||||||||||
Three Body Problem (CR3BP) model | ||||||||||||||||||
|
||||||||||||||||||
Features: | ||||||||||||||||||
1. Setup multiple nodes to use a Multiple Shooter Targeter | ||||||||||||||||||
2. (will be added)Periodic Orbit Multiple Shooter Targeter (Acts as a single shooter if n_node == 1) | ||||||||||||||||||
3. (will be added)Newton-Raphson Solver | ||||||||||||||||||
Comment on lines
+10
to
+13
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
|
||||||||||||||||||
|
||||||||||||||||||
References | ||||||||||||||||||
____________ | ||||||||||||||||||
Comment on lines
+15
to
+16
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
|
||||||||||||||||||
This work heavily relies on the work done by the various past and current members of the Multi-Body Dynamics Research Group and Prof. Kathleen C. Howell | ||||||||||||||||||
These are some of the referneces that provide a comprehensive brackground and have been the foundation for the work: | ||||||||||||||||||
1. E. Zimovan, "Characteristics and Design Strategies for Near Rectilinear Halo Orbits Within the Earth-Moon System," M.S., August 2017 | ||||||||||||||||||
2. E. Zimovan Spreen, "Trajectory Design and Targeting for Applications to the Exploration Program in Cislunar Space," Ph.D., May 2021 | ||||||||||||||||||
3. V. Szebehely, "Theory of Orbits: The Restricted Problem of Three Bodies", 1967 | ||||||||||||||||||
4. W. Koon, M. Lo, J. Marsden, S. Ross, "Dynamical Systems, The Three-Body Problem, and Space Mission Design", 2006 | ||||||||||||||||||
""" | ||||||||||||||||||
import numpy as np | ||||||||||||||||||
from cr3bp_master import prop_cr3bp | ||||||||||||||||||
|
||||||||||||||||||
|
||||||||||||||||||
def multi_shooter_nodes_setup_cr3bp( | ||||||||||||||||||
mu, ic, tf, n_node, node_place_opt="time", int_tol=1e-12 | ||||||||||||||||||
): | ||||||||||||||||||
""" | ||||||||||||||||||
Calculate states of n_nodes and time between each node | ||||||||||||||||||
Patch point placement strategy: | ||||||||||||||||||
1) Computes the n_nodes placed after NEARLY equal time intervals | ||||||||||||||||||
2) Computes the n_nodes placed after NEARLY equal time history INDEX | ||||||||||||||||||
(This might be better as more integration steps are taken in sensitive | ||||||||||||||||||
regions and placing the nodes by using index will place more points | ||||||||||||||||||
in the sensitive regions) | ||||||||||||||||||
|
||||||||||||||||||
First Node = I.C. of trajectory | ||||||||||||||||||
If symmetry is to be leverage to target a P.O.: last Node propagated by time 't' reaches the vicinity of the desired final state | ||||||||||||||||||
|
||||||||||||||||||
Dhruv Jain, Feb 18 2022 | ||||||||||||||||||
Comment on lines
+31
to
+43
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Just a nitpick:
Thus, make sure your docstrings have a one-line summary ending with period, followed by a blank line that introduces an additional summary. On top of that, there is no need to indent lists within docstrings. |
||||||||||||||||||
|
||||||||||||||||||
Parameters | ||||||||||||||||||
---------- | ||||||||||||||||||
mu : float, M2/(M1+M2) | ||||||||||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Only supported types can be included within docstrings. If you would like to include the mathematics behind the variable, you can do so in the description of the argument, right below this line:
Suggested change
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Apply to the rest of docstrings in your code. |
||||||||||||||||||
M1 and M2 are mass of Primary Bodies and M2<M1 | ||||||||||||||||||
ic : numpy ndarray (6x1), {Can handle all 42 states for CR3BP+STM integration} | ||||||||||||||||||
States are defined about the barycenter of the two primaries, P1 and P2 | ||||||||||||||||||
Initial condition: 6 states to compute a trajectory]; | ||||||||||||||||||
[0:x0, 1:y0, 2:z0, 3:vx0, 4:vy0, 5:vz0] [non-dimensional] [nd] | ||||||||||||||||||
tf : float | ||||||||||||||||||
Integration time [nd] | ||||||||||||||||||
Can be negative or positive, negative => Integration in backwards time | ||||||||||||||||||
n_node : int, Number of nodes | ||||||||||||||||||
IC is node 1 and nth node is node when propagated by tn value should reach the vicinity of the desired final state | ||||||||||||||||||
node_place_opt : string, optional | ||||||||||||||||||
'time': Computes the n_nodes placed after NEARLY equal time intervals | ||||||||||||||||||
'index': Computes the n_nodes placed after NEARLY equal time history INDEX | ||||||||||||||||||
int_tol : float, optional | ||||||||||||||||||
Absolute = Relative Integration Tolerance | ||||||||||||||||||
The default is 1e-12. | ||||||||||||||||||
|
||||||||||||||||||
Returns | ||||||||||||||||||
------- | ||||||||||||||||||
ic_node : list, ndarray, float64 | ||||||||||||||||||
Stores the IC of the n_node, where ic_node[0] = ic0 | ||||||||||||||||||
t_node : list, float | ||||||||||||||||||
Stores the time from one node to the next, t_node[-1]: time to reach the desired state (somekind of corrsing) | ||||||||||||||||||
""" | ||||||||||||||||||
|
||||||||||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. First line of code must start right below the docstring. No need to have this blank line. |
||||||||||||||||||
if node_place_opt != "time" and node_place_opt != "index": | ||||||||||||||||||
print( | ||||||||||||||||||
"Incorrect node placement option passed. Allowable options: time and index" | ||||||||||||||||||
) | ||||||||||||||||||
return 0 | ||||||||||||||||||
|
||||||||||||||||||
results_stm = prop_cr3bp( | ||||||||||||||||||
mu, ic, tf, stm_bool=0, xcross_cond=0 | ||||||||||||||||||
) # Propagate I.C. till tf | ||||||||||||||||||
ic_node = [] | ||||||||||||||||||
t_node = [] | ||||||||||||||||||
|
||||||||||||||||||
ic_node.append(ic) | ||||||||||||||||||
|
||||||||||||||||||
if node_place_opt == "time": | ||||||||||||||||||
# Time of each segment | ||||||||||||||||||
ti = np.linspace(0, tf, n_node + 1) | ||||||||||||||||||
ti = ti[ | ||||||||||||||||||
1: | ||||||||||||||||||
] # As ti is the time from node_i to next node, omit the first time, i.e. 0 | ||||||||||||||||||
|
||||||||||||||||||
for i in range(1, n_node): | ||||||||||||||||||
index = np.argmin( | ||||||||||||||||||
abs(results_stm["t"] - ti[i - 1]) | ||||||||||||||||||
) # compute index when index from time history when time is nearly | ||||||||||||||||||
ic_node.append(results_stm["states"][index, :]) | ||||||||||||||||||
t_node.append( | ||||||||||||||||||
results_stm["t"][index] - sum(t_node) | ||||||||||||||||||
) # t_node runs one node behind as it is the time to next node, subtract sum as ti+t2+t3 = tf | ||||||||||||||||||
|
||||||||||||||||||
t_node.append( | ||||||||||||||||||
tf - sum(t_node) | ||||||||||||||||||
) # Time from final node to vicinity of desired state | ||||||||||||||||||
|
||||||||||||||||||
elif node_place_opt == "index": | ||||||||||||||||||
num_index = len(results_stm["t"]) | ||||||||||||||||||
indices = np.linsapce( | ||||||||||||||||||
0, num_index, n_node, dtype="int" | ||||||||||||||||||
) # Linearly spaced indices | ||||||||||||||||||
|
||||||||||||||||||
for i in range(1, n_node): | ||||||||||||||||||
ic_node.append(results_stm["states"][indices[i], :]) | ||||||||||||||||||
t_node.append( | ||||||||||||||||||
results_stm["t"][indices[i]] - sum(t_node) | ||||||||||||||||||
) # t_node runs one node behind as it is the time to next node, subtract sum as ti+t2+t3 = tf | ||||||||||||||||||
|
||||||||||||||||||
t_node.append( | ||||||||||||||||||
tf - sum(t_node) | ||||||||||||||||||
) # Time from final node to vicinity of desired state | ||||||||||||||||||
|
||||||||||||||||||
return ic_node, t_node |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,173 @@ | ||
""" | ||
Created on 10 Nov 2021 22:43:38 2022 | ||
|
||
@author: Dhruv Jain, Multi-Body Dynamics Research Group, Purdue University | ||
[email protected] | ||
|
||
Objective: This file contains functions to compute various characteristic quantities | ||
for a Cirular Restircted Three Body Problem (CR3BP) model setup | ||
|
||
Features: | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. There is no need to indent these lines. |
||
1. Compute CR3BP 'mu', mass ratio of P1-P2 primary bodies | ||
2. Compute CR3BP 'l*', characterisitic lenght of P1-P2 primary bodies | ||
3. Comptue CR3BP 't*', charcterisitic time of P1-P2 primary bodies | ||
|
||
Useful ideas: | ||
1. To non-dimensioanlize a physical position vector: divide by l* | ||
2. To dimensioanlize a [nd] position vector: multiply by l* | ||
3. To non-dimensionalize a physical veloctiy vecotr, multiply by t*/l* | ||
4. To dimensionalize a [nd] veloctiy vecotr, multiply by l*/t* | ||
5. Position of P1 = [-mu, 0, 0] [nd] {x,y,z} | ||
5. Position of P2 = [1-mu, 0, 0] [nd] {x,y,z} | ||
|
||
Physical Qunatities obtained from: JPL’s ephemerides file de405.spk and https://ssd.jpl.nasa.gov/?planet_pos | ||
""" | ||
from poliastro.bodies import ( | ||
Earth, | ||
Jupiter, | ||
Mars, | ||
Mercury, | ||
Moon, | ||
Neptune, | ||
Pluto, | ||
Saturn, | ||
Sun, | ||
Uranus, | ||
Venus, | ||
) | ||
|
||
|
||
def mu_calc(mu_pi): | ||
"""Calculate mu of CR3BP | ||
Parameters | ||
Comment on lines
+41
to
+42
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Remember to add a blank line between the one-summary line and the rest of the sections. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Applies to the rest of the docstrings in this module. |
||
---------- | ||
mu_pi : ndarray, float | ||
mu_pi[0] = mu of P1 | ||
mu_pi[1] = mu of P2 | ||
|
||
Returns | ||
------- | ||
mu: float, M2/(M1+M2) | ||
M1 and M2 are mass of Primary Bodies and M2<M1 | ||
""" | ||
return mu_pi[1] / (mu_pi[0] + mu_pi[1]) | ||
|
||
|
||
def tstar_calc(mu_pi, dist): | ||
"""Calculate t* of CR3BP | ||
Parameters | ||
---------- | ||
dist: float, [nd] | ||
Non-dimensional distance between P1 and P2 | ||
mu_pi: ndarray, float | ||
mu_pi[0] = mu of P1 | ||
mu_pi[1] = mu of P2 | ||
|
||
Returns | ||
------- | ||
t*: float, [nd] | ||
sqrt(dist^3/m*) | ||
Non-dimensional time of P1-P2 system | ||
""" | ||
return (dist**3 / (mu_pi[0] + mu_pi[1])) ** 0.5 | ||
|
||
|
||
def bodies_char(b1, b2): | ||
"""Calculates mu, dist[km], t* [nd] of the 'body_i - body_j system' | ||
Parameters | ||
---------- | ||
b1 : string | ||
'body_i' | ||
b2 : string | ||
'body_j' | ||
|
||
Returns | ||
------- | ||
mu: float | ||
mu2/(mu1+mu2), mu2<mu1 | ||
dist: float, [nd] | ||
Non-dimensional distance between P1 and P2 | ||
tstar: float, [nd] | ||
sqrt(dist^3/m*) | ||
Non-dimensional time of P1-P2 system | ||
""" | ||
# Body values | ||
mu = {} # km^3 kg^-1 s^-2 | ||
mu["Sun"] = Sun.k.value * 1e-9 | ||
mu["Mercury"] = Mercury.k.value * 1e-9 | ||
mu["Venus"] = Venus.k.value * 1e-9 | ||
mu["Moon"] = Moon.k.value * 1e-9 | ||
mu["Earth"] = Earth.k.value * 1e-9 | ||
|
||
mu["Mars"] = Mars.k.value * 1e-9 | ||
mu["Jupiter"] = Jupiter.k.value * 1e-9 | ||
mu["Saturn"] = Saturn.k.value * 1e-9 | ||
mu["Uranus"] = Uranus.k.value * 1e-9 | ||
mu["Neptune"] = Neptune.k.value * 1e-9 | ||
mu["Pluto"] = Pluto.k.value * 1e-9 | ||
|
||
mu["Phobos"] = 0.0007112 # Phobos, GM | ||
mu["Titan"] = 8978.1382 # Titan, GM | ||
mu["Ganymede"] = 9887.834 # Ganymede, GM | ||
mu["Titania"] = 228.2 # Titania, GM | ||
mu["Triton"] = 1427.598 # Triton, GM | ||
mu["Charon"] = 102.30 # Charon, GM | ||
|
||
############# | ||
distances = {} # km | ||
distances["EarthMoon"] = 384400 | ||
distances["SunEarth"] = 149600000 | ||
|
||
distances["SunMars"] = 227944135 | ||
distances["SunJupiter"] = 778279959 | ||
distances["SunSaturn"] = 1427387908 | ||
distances["SunUranus"] = 2870480873 | ||
distances["SunNeptune"] = 4498337290 | ||
distances["SunPluto"] = 5907150229 | ||
|
||
distances["MarsPhobos"] = 9376 | ||
distances["JupiterGanymede"] = 1070400 | ||
distances["SaturnTitan"] = 1221865 | ||
distances["UranusTitania"] = 436300 | ||
distances["NeptuneTriton"] = 354759 | ||
distances["PlutoCharon"] = 17536 | ||
Comment on lines
+116
to
+133
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I suppose these are mean distances? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Yes. they are all the mean distances used for CR3BP models. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Thanks, @JackCrusoe47! Yes, these are mean distances obtained from the JPL references stated at the top of the code. @astrojuanlu does poliastro have the mean distances saved anywhere? If not, do you think we should look into adding these constants to constants.py like file.
Comment on lines
+117
to
+133
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Some of these distances are defined in |
||
######### | ||
|
||
mu_pi = [] | ||
|
||
try: | ||
temp1 = mu[b1] | ||
temp2 = mu[b2] | ||
if temp1 == temp2: | ||
print( | ||
"Same bodies passed as P1 and P2. Please pass the secondary body of P1 as P2" | ||
) | ||
return 0, 0, 0 | ||
except KeyError: | ||
print("KeyError-> Incorrect/Does Not Exist Input Bodies name") | ||
return 0, 0, 0 | ||
Comment on lines
+141
to
+148
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Another nitpick. Although very useful for quick debugging, printing is not the most "professional" way of retrieving information about the status of a program. As an alternative, consider raising logging, raising warnings or even errors. Remember that you can customize the message of an error in Python by doing something like: raise ValueError("Please, pass the right parameters...") |
||
|
||
# Sort P1 and P2 based on their mu values | ||
if mu[b1] >= mu[b2]: # b1 is the bigger primary | ||
mu_pi.append(mu[b1]) | ||
mu_pi.append(mu[b2]) | ||
bodies = b1 + b2 | ||
else: # b2 is the bigger primary | ||
mu_pi.append(mu[b2]) | ||
mu_pi.append(mu[b1]) | ||
|
||
# To create string to get distance | ||
bodies = b2 + b1 | ||
|
||
mu = mu_calc(mu_pi) | ||
|
||
try: | ||
dist = distances[bodies] | ||
except KeyError: | ||
print( | ||
"KeyError-> Error in combination of bodies P1-P2, typo/DNE/combo not created" | ||
) | ||
return 0, 0, 0 | ||
else: | ||
tstar = tstar_calc(mu_pi, dist) | ||
return mu, dist, tstar |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That's great! ❤️