Dev/orbit propagator - #89
Conversation
…; added constant e in constants.py
CHANGE: perturbation calculation in encke_motion and cowell_motion
| """ | ||
| x, y, z = r_eci_m | ||
| r_norm_m = np.linalg.norm(r_eci_m) | ||
| mu_m3_s2 = 3.986004418e14 |
There was a problem hiding this comment.
Should have these kinds of constants available in a file somewhere or at least defined at the top level. I don't want to have floating definitions in different places because you can end up using slightly different values. I think Jenny already implemented this somewhere.
| if altitude_m < 0.0: | ||
| raise ValueError("Altitude must be nonnegative.") | ||
|
|
||
| longitude_deg = np.degrees(np.arctan2(y, x)) |
There was a problem hiding this comment.
This evaluates to a pair of spherical coordinates in ECI, but lat/lon usually refer to ECEF. I'm not sure off hand if it matters for NRLMSIS, please double check and apply changes if needed. If it doesn't just let me know in a reply and re-submit the review.
If you need true lat/lon, we should implement a separate function to retrieve that (it's going to be an involved one).
There was a problem hiding this comment.
This has been fixed in the updated environment.py code fyi
| z2_over_r2 = (z**2) / (r_norm_m**2) | ||
| factor = - 3 / 2 * (mu_m3_s2 * j2 * r_eq_m**2) / (r_norm_m**5) | ||
|
|
||
| ax = factor * x * (1 - 5 * z2_over_r2) |
There was a problem hiding this comment.
J2 perturbation is symmetric about the planet's pole, which will not be aligned with the ECI z-axis at any given time. I think we should keep this as a fallback for now, but for better fidelity we should implement a version that accounts for the Earth orientation parameters and transforms the state into an appropriate coordinate system with a pole-aligned z-axis. Resultant acceleration should still be output in ECI for simplicity/compatibility.
| from environment import j2_acceleration_m_s2, aerodynamic_drag_perturbation_m_s2 | ||
|
|
||
|
|
||
| class simulation_config(): |
There was a problem hiding this comment.
Class appears to be misconfigured.
I would probably just decorate this as a dataclass and forego the init(). Treat it like a struct
| """ | ||
| t0: datetime # Simulation start time (UTC) | ||
| tf: datetime # Simulation end time (UTC) | ||
| time_steps: int = 1000 # Number of time steps |
There was a problem hiding this comment.
Our preference will probably be to use an adaptive step size integrator like RK45. It's good to have this field in case we do use a fixed step integrator, but make sure it's properly handled either way.
| r_mag = np.linalg.norm(r_vec) # magnitude of r vector | ||
|
|
||
| kep = cartesian2keplerian(r, mu) # for now assume we are only doing encke motion for orbit around sun | ||
| kep = cartesian2keplerian(r, mu) # for now assume we are only doing encke motion for orbit around sun (but mu was calculated with earth's mass??????????) |
There was a problem hiding this comment.
What's this about orbiting the Sun? We're not orbiting the Sun
| Returns: | ||
| x: (np.ndarray) (time_steps, 6) Array of orbital states at each time step. | ||
| """ | ||
| if config.propagator_method == "cowell": |
There was a problem hiding this comment.
This is the kind of messiness I was referring to when using str config. If you provide the propagator function directly, then all you need to do is ensure the signature is consistent across models and then you can call solve_ivp(fun=propagator, ...).
Realistically, though, we're going to implement a higher level function that encompasses all aspects of the simulation, including both attitude and orbit. So you would have something like motion_model(t, x) which goes over perturbations, control inputs, attitude motion, orbit motion, and anything else that needs to be calculated per timestep and you feed that into solve_ivp(fun=motion_model, ...).
| tf: datetime # Simulation end time (UTC) | ||
| time_steps: int = 1000 # Number of time steps | ||
| propagator_method: str = "cowell" # Propagator to use: "cowell", "encke", or "sgp4" | ||
| x0: np.ndarray = np.array([0., 0., 0., 0., 0., 0.]) # Initial state vector (6x1) (must be in ECI (for now)) |
There was a problem hiding this comment.
Also needs attitude state and an attitude motion model
| """ | ||
| Configuration for the simulation. | ||
| """ | ||
| t0: datetime # Simulation start time (UTC) |
There was a problem hiding this comment.
Should include a spacecraft object. We haven't defined this yet, but it will contain all the information about the sensors and actuators, controller, state estimator, physical properties, etc.
| from environment import j2_acceleration_m_s2, aerodynamic_drag_perturbation_m_s2 | ||
|
|
||
|
|
||
| class simulation_config(): |
There was a problem hiding this comment.
Please extract this and the orbit motion model into something like simulator.py. Remember that the simulator covers more than just the orbital mechanics
ADD: cowell motion unit test
Implemented code for issue #78.
ADD: test for aerodynamic_drag_perturbation_m_s2
5 participants
No description provided.