Add CR3BP functions created by Dhruv Jain

contrib/cr3bp_DhruvJ/README.md

@@ -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!

contrib/cr3bp_DhruvJ/cr3bp_PO_targeter.py

@@ -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
+
+References
+____________
+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
+
+ Parameters
+ ----------
+ mu : float, M2/(M1+M2)
+ 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)
+ """
+
+ 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

contrib/cr3bp_DhruvJ/cr3bp_char_quant.py

@@ -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:
+ 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
+ ----------
+ 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
+ #########
+
+ 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
+
+ # 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