---
title: "ISS Eclipse Determination"
description: |
  Determining how much sunlight a body is receiving.
draft: false
author:
  - name: Anson Biggs
    url: https://ansonbiggs.com
repository_url: https://gitlab.com/lander-team/air-prop-simulation
date: 04-01-2021
fig_width: 6
fig_align: "center"
output:
  distill::distill_article:
    self_contained: false
categories:
  - Julia
  - Astrodynamics
#bibliography: ../../citations.bib
creative_commons: CC BY
preview: preview.png
---

Determining the eclipses a satellite will encounter is a major driving factor when designing a mission in space. Thermal and power budgets have to be made with the fact that a satellite will periodically be in the complete darkness of space with no solar radiation to power the solar panels and keep the spacecraft from freezing. 

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, results = 'hide')
library(JuliaCall)
#julia_setup(JULIA_HOME = "/opt/julia-1.6.1/bin/")
julia_setup(installJulia = TRUE)
```

## What is an Eclipse

![Geometry of an Eclipse](geometry.svg)


![Body Radius's and Position Vectors](vectors_radiuss.svg)


## The Code

```{julia, code_folding=TRUE}
using Unitful
using LinearAlgebra
using SatelliteToolbox
using Plots
using Colors
theme(:ggplot2)
```


```{julia}
ISS = tle"""
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
"""
```

```{julia}
orbit = init_orbit_propagator(Val(:twobody), ISS[1]);
time = 0:0.1:((24 / ISS[1].n) .* 60 * 60);
o, r, v = propagate!(orbit, time);
```


```{julia}
function sunlight(Rbody, r_sun_body, r_body_sc)
    Rsun = 695_700u"km"
    
    hu = Rbody * norm(r_sun_body) / (Rsun - Rbody)
    
    θe = acos((r_sun_body ⋅ r_body_sc) / (norm(r_sun_body) * norm(r_body_sc)))

    θu = atan(Rbody / hu)
    du = hu * sin(θu) / sin(θe + θu)

    θp = π - atan(norm(r_sun_body) / (Rsun + Rbody))
    dp = Rbody * sin(θp) / cos(θe - θp)

    S = 1
    if (θe < π / 2) && (norm(r_body_sc) < du)
        S = 0
    end
    if (θe < π / 2) && ((du < norm(r_body_sc)) && (norm(r_body_sc) < dp))
        S = (norm(r_body_sc .|> u"km") - du) / (dp - du) |> ustrip
    end

    return S
end
```

```{julia}
S = r .|> R -> sunlight(6371u"km", [0.5370, 1.2606, 0.5466] .* 1e8u"km", R .* u"m")
```

## Plotting the Results

```{julia, code_folding=TRUE, results='show',layout="l-body-outset",fig.cap= "Rocket Motor Data: [@thrustcurve]", preview=TRUE }
light_range = range(colorant"black", stop = colorant"yellow", length = 101);
light_colors = [light_range[unique(round(Int, 1 + s * 100))][1] for s in S];

plot(
    LinRange(0, 24, length(S)),
    S .* 100,
    linewidth = 5,
    legend = false,
    color = light_colors,
);

xlabel!("Time (hr)");
ylabel!("Sunlight (%)");
title!("ISS Sunlight Over a Day")
```