Implement sun eclipse duration rollup [feature request]

[SoRobby]

Want to start off by saying thanks for providing such an amazing astrodynamics library and your great YouTube videos!

Eclipse Duration Calculator
Implement functionality to calculate the average and max eclipse time. Furthermore you can add the ability to plot the eclipse durations for a given orbit. Potentially add the ability to add other eclipse objects like the moon.

[alfonsogonzalez]

@SoRobby Thank you for reaching out! This is a good idea. Here are my initial thoughts:

This could be implemented in 2 ways.

Just like the stop conditions functions, the ODE solver can take in calc_umbra and calc_penumbra functions that equal -1 if False, and 1 if True, thus when the spacecraft goes in and out of eclipse, and event will be triggered (although won't stop the propagation), and then afterwards the ode_sol object can be parsed to calculate the entry and exit times, and other calculations that you've mentioned.

There can also be a function in the orbit_calculations.py that takes in position and eclipsing body as inputs, and outputs whether the spacecraft is in umbra, penumbra, or neither. This would be more of a post processing implementation that would be similar to JPL's MONTE event finder that can take in a trajectory and output entrances and exits of eclipse.

There are a few SPICE functions that deal with lighting calculations (et2lst, subsol, and others) but I don't think they will be very useful in this case. So I'll have to do a little bit of reading on how to implement umbra / penumbra geometry and then implement it.

Again, thank you for creating this Issue! Its an interesting problem and I'm going to learn lots from solving it.

[SoRobby]

Both are good ways it could be implemented!

I'd personally push for having it as a post processor. This keeps the ODE solver more lightweight, especially if you want to implement other functions such as link budget, attitude control, thermal simulations, power / battery modeling, modeling spacecraft operational states, etc... I believe this approach offers greater flexibility, below are some examples on some use cases:

• Scenario 1 - Positional data can be post processed and fed directly into post processing functions like calc_max_eclipse or calc_average_eclipse or into other post processor functions (i.e. link budget).
• Scenario 2 - Data can be output to a json / txt file. You can load that information into post processor functions like an eclipse calculator, link budget, or power analysis tool. Since these functions would be post processed, you wouldn't have to re-run the entire simulation if you change for example, the link budgets frequency or transmitting input power.

There's certainty a lot of complexities associated with calculating eclipse times, from what I've seen in your library a lot of the groundwork is there! Keep up the awesome work!!

Here is another good diagram that illustrates the geometry of this a bit better:

[alfonsogonzalez]

Thank you for this info! My knowledge of full system analysis (like link budgets and thermal requirements) are limited so info like this is very useful.

The scenario 1 functions you mentioned make a lot of sense. I think I am going to make the next few videos in the orbital mechanics with python series explaining the geometry and then the software for those functions.

For scenario 2, my initial thoughts are that it would be best to save / read the trajectory as a SPICE kernel, since there are SPICE functions that take care of interpolation between data points.

I'm wondering though how precise the eclipse start/stop times should be. If one is okay with whatever the output state/time that the ODE solver used, this function can be more straightforward. But the function can get more precise by using a root solver method (using finite differences to estimate the derivative) to more precisely find the entrance/exit times. The SciPy ODE solver already has this implemented for event functions.

Overall, I think that both implementations would make sense because it depends on the specific needs of the user. And it doesn't necessarily add any overhead when implemented into the ODE, because its an optional argument to add event functions to the ODE solver (part of the config dictionary / stop conditions implementation in the Spacecraft class).

I'm going to get started on it this weekend and I'll keep you updated!

[alfonsogonzalez]

@SoRobby I'm likely going to be able to finish this up tomorrow, but for now check out this NASA paper (old but gold): https://ntrs.nasa.gov/api/citations/19950023025/downloads/19950023025.pdf

That paper does a good job of going over the geometry of this problem. I'll be using the assumptions that Earth is spherical and not considering atmospheric diffraction effects. However, I also found this PhD thesis that goes deep into the math when you don't make those limiting assumptions: https://www.researchgate.net/publication/323184678_Highly_Physical_Solar_Radiation_Pressure_Modeling_During_Penumbra_Transitions

[SoRobby]

@alfonsogonzalez Excellent references, it really illustrates how complicated calculating the eclipse duration can be depending on the assumptions you make. You can potentially use these sources to help validate it too.

I think those are excellent assumptions, especially if you want it to be applicable to most celestial bodies as a first order estimate / calculation. This is going to be amazing, once implemented I'd love to contribute and add some models to define the spacecraft side of things from power, attitude control, spacecraft operational states (system modes), telecom links, etc... There are lots of spacecraft related models that can be implemented from SMAD (Space Mission Engineering: The New SMAD) with some validation cases (note: important to note that some of SMAD's models are only applicable to Earth-based missions, models like SMAD's solar array become a bit unrealistic for sizing scenarios beyond 2-3 AU).

Keep me posted and let me know if you want to discuss any of this - it can be difficult developing / planning something like this alone.

[alfonsogonzalez]

@SoRobby Preliminary software has been added to this branch: https://github.com/alfonsogonzalez/AWP/tree/eclipses

I will be refining / updating as I make the video this weekend

[alfonsogonzalez]

For reference, here are some plots to explain the geometry (of course will be discussed in detail in the video):

[SoRobby]

@alfonsogonzalez Awesome work on developing this so quickly! I went through the code a bit this afternoon - currently no immediate feedback, but one question (below). Looking forward to the video!

I forked the build and added a very basic calculation tool to calculate solar panel power generation - feel free to integrate - eventually want to add more complex analysis and models to model the entire spacecraft and it's subsystems.
Forked Repo: https://github.com/SoRobby/AWP/blob/eclipses/src/python_tools/spacecraft_toos.py

I do have one question,
2 = Umbra (total eclipse)
1 = Penumbra (partial eclipse) --> should this actually be between 0 and 1?
0 to -1 = In direct sun?

Was reading through another paper on how they denote umbra, penumbra, and direct sun (full phase).

paper: https://discovery.ucl.ac.uk/id/eprint/10059978/1/eclipse_model_revision3.pdf

Looking forward to implementing more complex spacecraft models!

[alfonsogonzalez]

The reason I implemented it as 2,1,-1 is because if you do it this way, each event has a unique number value assigned to it (no eclipse --> umbra, no eclipse --> penumbra, umbra --> umbra, umbra --> no eclipse, etc.). This is not the same as modeling for SRP where the shadow function is a continuous value. This method is useful for the eclipse finding function, which looks for values that indicate specific events (for example, no eclipse to umbra = 3, no eclipse to penumbra = 2, umbra to penumbra = -1, etc.). I found that the logic for doing this was a bit more complex than I initially expected, and these values ended up solving it.

Also, through testing I found that its most common that the finder just finds no eclipse --> penumbra. Also, the difference for comparing to the SPICE function gfoclt is 2-5 seconds, which for thermal modeling probably doesn't really make that much difference. So you can just use SRP = 0 for umbra for a good first order approximation using this tool.

I'll explain this in detail in the software video. For now, here is the first video that goes over the math and lunar, geostationary, sun synchronour terminator, and ISS orbits (software video probably coming next weekend): https://www.youtube.com/watch?v=LeJpWfyXJPw

For now, I only want to merge software that I've written and made videos covering for this repository. The videos, this repo, and the book are all meant to be as learning tools for astrodynamics, so I want to implement the 3 of them together (I have a lot of catching up to do with the book since this was a recent idea). I do eventually plan on implementing ADCS but I don't want to post it before making a video. I already have a version of the Spacecraft class with it implemented, but I have found that I don't truly and wholly understand a problem / math until I am able to teach it, so as of now I don't feel comfortable posting it.

I'm going to keep this Issue open until I finish with the software video, where I will likely post more scripts / tests