diff --git a/_posts/2021-04-01-air-propulsion-simulation/air-propulsion-simulation.Rmd b/_posts/2021-04-01-air-propulsion-simulation/air-propulsion-simulation.Rmd index 7a71b12..6282d4b 100644 --- a/_posts/2021-04-01-air-propulsion-simulation/air-propulsion-simulation.Rmd +++ b/_posts/2021-04-01-air-propulsion-simulation/air-propulsion-simulation.Rmd @@ -24,7 +24,7 @@ For Capstone my team was tasked with designing a system capable of moving mining ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, results = 'hide') library(JuliaCall) -#julia_setup(JULIA_HOME = "/opt/julia-1.6.0/bin/") +julia_setup(JULIA_HOME = "/opt/julia-1.6.0/bin/") ``` ```{julia, code_folding=TRUE} diff --git a/_posts/2021-04-01-air-propulsion-simulation/air-propulsion-simulation.html b/_posts/2021-04-01-air-propulsion-simulation/air-propulsion-simulation.html index 21f980c..276f3bb 100644 --- a/_posts/2021-04-01-air-propulsion-simulation/air-propulsion-simulation.html +++ b/_posts/2021-04-01-air-propulsion-simulation/air-propulsion-simulation.html @@ -110,6 +110,8 @@ code span.wa { color: #5e5e5e; font-style: italic; } /* Warning */ + + @@ -1456,7 +1458,7 @@ code span.wa { color: #5e5e5e; font-style: italic; } /* Warning */ @@ -1473,13 +1475,13 @@ code span.wa { color: #5e5e5e; font-style: italic; } /* Warning */
Determining how much sunlight a body is receiving.
+Determining how much sunlight a body orbiting a planet is receiving.
The above image is a simple representation of what an eclipse is. You’ll notice there is the Umbra which is complete darkness, then the Penumbra which is a shadow of varying darkness, and then the rest of the orbit is in complete sunlight. For this example I will be using the ISS which has a very low orbit so the Penumbra isn’t much of a problem. You can tell by looking at the diagram that higher altitude orbits would spend more time in the Penumbra.
Here is a more detailed view of the eclipse that will make it easier to explain the code. There are 2 Position vectors and 2 radius’s that need to be known for simple eclipse determination. There are more advanced cases where the atmosphere of the body your orbiting can greatly affect the Umbra and Penumbra, and other bodies could also potentially block the Sun, but for this example we will keep it simple since those have very little affect for the ISS’s orbit. Rsun
and Rbody
are the radius’s of the Sun and Body (In this case Earth) respectively. r_sun_body
is a vector from the center of the Sun to the center of the target body. For this example I will only be using one vector, but for more rigorous eclipse determination its important to calculate this at least once a day since it does significantly change over the course of a year. The reason that I am ignoring it at the moment is because there is currently no good way to calculate Ephemerides in Julia but the package is being worked on so I may revisit this and do a more rigorous analysis in the future. r_body_sc
is a position vector from the center of the body being orbitted, to the center of our spacecraft.
In order to get the orbit for the ISS I used a Two-Line Element which is a data format for explaining orbits. the US Joint Space Operations Center makes these widely available, but https://live.ariss.org/tle/ makes the TLE for the ISS way more accessible (“ARISS TLE,” n.d.). The Julia Package SatelliteToolbox.jl makes it super easy to turn a TLE into an orbit that can be propagated. Simply putting the TLE in a string and using the tle
string macro like below and now we have access to the information to start making our ISS orbit.
= tle"""
ISS ISS (ZARYA)
1 25544U 98067A 21103.84943184 .00000176 00000-0 11381-4 0 9990
2 25544 51.6434 300.9481 0002858 223.8443 263.8789 15.48881793278621
"""
1-element Vector{TLE}:
+ TLE: ISS (ZARYA) (Epoch = 2021-04-13T20:23:10.911)
+Now that we have the TLE we can pass that into SatelliteToolbox’s orbit propagator. Before we can propagate the orbit we need to have a range of time steps to pass into the propagator. The TLE gives the mean motion, n, which is the revolutions per day so using that we can calculate the amount of time required for one orbit which is all that were worried about for this analysis. The propagator returns a tuple containing the Orbital elements, a position vector with units meters, and a velocity vector with units meters per second. For this analysis were only worried about the position vector.
+ISS[1].n
= init_orbit_propagator(Val(:twobody), ISS[1]);
- orbit = 0:0.1:((24 / ISS[1].n) .* 60 * 60);
- time , r, v = propagate!(orbit, time); o
= init_orbit_propagator(Val(:twobody), ISS[1]);
+ orbit = 0:0.1:((24 / ISS[1].n) .* 60 * 60); # ISS[1].n gives the mean motion, or orbits per day.
+ time , r, v = propagate!(orbit, time); o
Now we just need way to use the radii and vectors discussed earlier to determine if the ISS is in the penumbra or umbra. This is a lot of pretty basic trigonometry and vector math.
+add more discussion about the math
function sunlight(Rbody, r_sun_body, r_body_sc)
-= 695_700u"km"
- Rsun
- = Rbody * norm(r_sun_body) / (Rsun - Rbody)
- hu
- = acos((r_sun_body ⋅ r_body_sc) / (norm(r_sun_body) * norm(r_body_sc)))
- θe
-= atan(Rbody / hu)
- θu = hu * sin(θu) / sin(θe + θu)
- du
-= π - atan(norm(r_sun_body) / (Rsun + Rbody))
- θp = Rbody * sin(θp) / cos(θe - θp)
- dp
-= 1
- S if (θe < π / 2) && (norm(r_body_sc) < du)
- = 0
- S end
- if (θe < π / 2) && ((du < norm(r_body_sc)) && (norm(r_body_sc) < dp))
- = (norm(r_body_sc .|> u"km") - du) / (dp - du) |> ustrip
- S end
-
-return S
- end
function sunlight(Rbody, r_sun_body, r_body_sc)
+= 695_700u"km"
+ Rsun
+ = Rbody * norm(r_sun_body) / (Rsun - Rbody)
+ hu
+ = acos((r_sun_body ⋅ r_body_sc) / (norm(r_sun_body) * norm(r_body_sc)))
+ θe
+= atan(Rbody / hu)
+ θu = hu * sin(θu) / sin(θe + θu)
+ du
+= π - atan(norm(r_sun_body) / (Rsun + Rbody))
+ θp = Rbody * sin(θp) / cos(θe - θp)
+ dp
+= 1
+ S if (θe < π / 2) && (norm(r_body_sc) < du)
+ = 0
+ S end
+ if (θe < π / 2) && ((du < norm(r_body_sc)) && (norm(r_body_sc) < dp))
+ = (norm(r_body_sc .|> u"km") - du) / (dp - du) |> ustrip
+ S end
+
+return S
+ end
Then we can pass all the values we’ve gathered into the function we just made.
= r .|> R -> sunlight(6371u"km", [0.5370, 1.2606, 0.5466] .* 1e8u"km", R .* u"m") S
= r .|> R -> sunlight(6371u"km", [0.5370, 1.2606, 0.5466] .* 1e8u"km", R .* u"m"); S
The sunlight
function returns values from 0 to 1, 0 being complete darkness, 1 being complete sunlight, and anything between being the fraction of light being received. Again since the ISS has a very low orbit, the amount of time spend in the penumbra is almost insignificant.
= range(colorant"black", stop = colorant"yellow", length = 101);
- light_range = [light_range[unique(round(Int, 1 + s * 100))][1] for s in S];
- light_colors
-
- plot(LinRange(0, 24, length(S)),
- .* 100,
- S = 5,
- linewidth = false,
- legend = light_colors,
- color ;
- )
-!("Time (hr)");
- xlabel!("Sunlight (%)");
- ylabel!("ISS Sunlight Over a Day") title
# Get fancy with the line color.
+= range(colorant"black", stop = colorant"orange", length = 101);
+ light_range = [light_range[unique(round(Int, 1 + s * 100))][1] for s in S];
+ light_colors
+
+ plot(LinRange(0, 24, length(S)),
+ .* 100,
+ S = 5,
+ linewidth = false,
+ legend = light_colors,
+ color ;
+ )
+!("Time (hr)");
+ xlabel!("Sunlight (%)");
+ ylabel!("ISS Sunlight Over a Day") title
-Figure 1: Rocket Motor Data: [@thrustcurve] +Figure 1: ISS Sunlight
Looking at the plot its pretty easy to see by the vertical transition from 0% to 100% that the time in the penumbra is limited, but almost counterintutively it also looks like the ISS gets more sunlight than it does darkness. Using the raw sunlight data we can actually calculate almost exactly how much time is spent in each region.
+Time in Sun: += length(S[S.==1])/length(S) * 100 sun
62.03323593209401
+= length(S[S.==0])/length(S) * 100 umbra
37.64408511553699
+= 100 - umbra - sun penumbra
0.322678952369003
+The ISS spends about 62% of its time in the sun, this is because if you go back and reference the diagram at the beginning of this post you can see that the umbra is actually a cone. This is mainly due to the fact that the Sun is massive compared to the Earth, but this effect is also stronger with orbits of higher altitudes.
+
If you see mistakes or want to suggest changes, please create an issue on the source repository.
The above image is a simple representation of what an eclipse is. You’ll notice there is the Umbra which is complete darkness, then the Penumbra which is a shadow of varying darkness, and then the rest of the orbit is in complete sunlight. For this example I will be using the ISS which has a very low orbit so the Penumbra isn’t much of a problem. You can tell by looking at the diagram that higher altitude orbits would spend more time in the Penumbra.
Here is a more detailed view of the eclipse that will make it easier to explain the code. There are 2 Position vectors and 2 radius’s that need to be known for simple eclipse determination. There are more advanced cases where the atmosphere of the body your orbiting can greatly affect the Umbra and Penumbra, and other bodies could also potentially block the Sun, but for this example we will keep it simple since those have very little affect for the ISS’s orbit. Rsun
and Rbody
are the radius’s of the Sun and Body (In this case Earth) respectively. r_sun_body
is a vector from the center of the Sun to the center of the target body. For this example I will only be using one vector, but for more rigorous eclipse determination its important to calculate this at least once a day since it does significantly change over the course of a year. The reason that I am ignoring it at the moment is because there is currently no good way to calculate Ephemerides in Julia but the package is being worked on so I may revisit this and do a more rigorous analysis in the future. r_body_sc
is a position vector from the center of the body being orbitted, to the center of our spacecraft.
In order to get the orbit for the ISS I used a Two-Line Element which is a data format for explaining orbits. the US Joint Space Operations Center makes these widely available, but https://live.ariss.org/tle/ makes the TLE for the ISS way more accessible (“ARISS TLE,” n.d.). The Julia Package SatelliteToolbox.jl makes it super easy to turn a TLE into an orbit that can be propagated. Simply putting the TLE in a string and using the tle
string macro like below and now we have access to the information to start making our ISS orbit.
= tle"""
ISS ISS (ZARYA)
1 25544U 98067A 21103.84943184 .00000176 00000-0 11381-4 0 9990
2 25544 51.6434 300.9481 0002858 223.8443 263.8789 15.48881793278621
"""
1-element Vector{TLE}:
+ TLE: ISS (ZARYA) (Epoch = 2021-04-13T20:23:10.911)
+Now that we have the TLE we can pass that into SatelliteToolbox’s orbit propagator. Before we can propagate the orbit we need to have a range of time steps to pass into the propagator. The TLE gives the mean motion, n, which is the revolutions per day so using that we can calculate the amount of time required for one orbit which is all that were worried about for this analysis. The propagator returns a tuple containing the Orbital elements, a position vector with units meters, and a velocity vector with units meters per second. For this analysis were only worried about the position vector.
+ISS[1].n
= init_orbit_propagator(Val(:twobody), ISS[1]);
- orbit = 0:0.1:((24 / ISS[1].n) .* 60 * 60);
- time , r, v = propagate!(orbit, time); o
= init_orbit_propagator(Val(:twobody), ISS[1]);
+ orbit = 0:0.1:((24 / ISS[1].n) .* 60 * 60); # ISS[1].n gives the mean motion, or orbits per day.
+ time , r, v = propagate!(orbit, time); o
Now we just need way to use the radii and vectors discussed earlier to determine if the ISS is in the penumbra or umbra. This is a lot of pretty basic trigonometry and vector math.
+add more discussion about the math
function sunlight(Rbody, r_sun_body, r_body_sc)
-= 695_700u"km"
- Rsun
- = Rbody * norm(r_sun_body) / (Rsun - Rbody)
- hu
- = acos((r_sun_body ⋅ r_body_sc) / (norm(r_sun_body) * norm(r_body_sc)))
- θe
-= atan(Rbody / hu)
- θu = hu * sin(θu) / sin(θe + θu)
- du
-= π - atan(norm(r_sun_body) / (Rsun + Rbody))
- θp = Rbody * sin(θp) / cos(θe - θp)
- dp
-= 1
- S if (θe < π / 2) && (norm(r_body_sc) < du)
- = 0
- S end
- if (θe < π / 2) && ((du < norm(r_body_sc)) && (norm(r_body_sc) < dp))
- = (norm(r_body_sc .|> u"km") - du) / (dp - du) |> ustrip
- S end
-
-return S
- end
function sunlight(Rbody, r_sun_body, r_body_sc)
+= 695_700u"km"
+ Rsun
+ = Rbody * norm(r_sun_body) / (Rsun - Rbody)
+ hu
+ = acos((r_sun_body ⋅ r_body_sc) / (norm(r_sun_body) * norm(r_body_sc)))
+ θe
+= atan(Rbody / hu)
+ θu = hu * sin(θu) / sin(θe + θu)
+ du
+= π - atan(norm(r_sun_body) / (Rsun + Rbody))
+ θp = Rbody * sin(θp) / cos(θe - θp)
+ dp
+= 1
+ S if (θe < π / 2) && (norm(r_body_sc) < du)
+ = 0
+ S end
+ if (θe < π / 2) && ((du < norm(r_body_sc)) && (norm(r_body_sc) < dp))
+ = (norm(r_body_sc .|> u"km") - du) / (dp - du) |> ustrip
+ S end
+
+return S
+ end
Then we can pass all the values we’ve gathered into the function we just made.
= r .|> R -> sunlight(6371u"km", [0.5370, 1.2606, 0.5466] .* 1e8u"km", R .* u"m") S
= r .|> R -> sunlight(6371u"km", [0.5370, 1.2606, 0.5466] .* 1e8u"km", R .* u"m"); S
The sunlight
function returns values from 0 to 1, 0 being complete darkness, 1 being complete sunlight, and anything between being the fraction of light being received. Again since the ISS has a very low orbit, the amount of time spend in the penumbra is almost insignificant.
= range(colorant"black", stop = colorant"yellow", length = 101);
- light_range = [light_range[unique(round(Int, 1 + s * 100))][1] for s in S];
- light_colors
-
- plot(LinRange(0, 24, length(S)),
- .* 100,
- S = 5,
- linewidth = false,
- legend = light_colors,
- color ;
- )
-!("Time (hr)");
- xlabel!("Sunlight (%)");
- ylabel!("ISS Sunlight Over a Day") title
# Get fancy with the line color.
+= range(colorant"black", stop = colorant"orange", length = 101);
+ light_range = [light_range[unique(round(Int, 1 + s * 100))][1] for s in S];
+ light_colors
+
+ plot(LinRange(0, 24, length(S)),
+ .* 100,
+ S = 5,
+ linewidth = false,
+ legend = light_colors,
+ color ;
+ )
+!("Time (hr)");
+ xlabel!("Sunlight (%)");
+ ylabel!("ISS Sunlight Over a Day") title
-Figure 1: Rocket Motor Data: [@thrustcurve] +Figure 1: ISS Sunlight
Looking at the plot its pretty easy to see by the vertical transition from 0% to 100% that the time in the penumbra is limited, but almost counterintutively it also looks like the ISS gets more sunlight than it does darkness. Using the raw sunlight data we can actually calculate almost exactly how much time is spent in each region.
+Time in Sun: += length(S[S.==1])/length(S) * 100 sun
62.03323593209401
+= length(S[S.==0])/length(S) * 100 umbra
37.64408511553699
+= 100 - umbra - sun penumbra
0.322678952369003
+The ISS spends about 62% of its time in the sun, this is because if you go back and reference the diagram at the beginning of this post you can see that the umbra is actually a cone. This is mainly due to the fact that the Sun is massive compared to the Earth, but this effect is also stronger with orbits of higher altitudes.
+