Tag Archives: SatelliteToolbox.jl

O pacote SatelliteToolbox.jl para Linguagem Julia

Olá!

Neste post, gostaria de apresentar o SatelliteToolbox.jl, que é um pacote para a linguagem Julia com muitas opções para analisar missões espaciais. Ele é usado diariamente no Instituto Nacional de Pesquisas Espaciais (INPE). Primeiramente, apresento um breve histórico sobre o pacote e, em seguida, mostro algumas análises interessantes que podem ser feitas com ele.

Histórico

Em 2013, ingressei no INPE como Engenheiro de Sistemas Espaciais Júnior. Fui designado para a divisão de sistemas espaciais onde tive que trabalhar com o subsistema de controle de atitude e órbita (AOCS). Como eu tinha apenas um conhecimento intermediário sobre órbitas, decidi me aprofundar nesse assunto codificando algoritmos e comparando os resultados com a herança do INPE e com a literatura no meu tempo livre.

O primeiro passo foi selecionar a linguagem! Em meu Doutorado, utilizei o MATLAB para simular um sistema de navegação inercial, mas as simulações de Monte Carlo eram tão lentas que precisei reescrever muitas partes em C usando CMEX. Por outro lado, na minha pesquisa de pós-doutorado, onde estudei estimação em sistemas distribuídos não-lineares, decidi usar FORTRAN (com os padrões de 2008 para ter pelo menos um código legível …) para que a velocidade de execução não fosse um problema. Sim, o desempenho foi muito bom, mas demorei muito tempo para codificar. Então ouvi falar de uma nova linguagem que prometia o melhor dos dois mundos: algo que se assemelhava a uma linguagem interpretada com a velocidade de uma linguagem compilada! E foi assim que conheci Julia.

Naquela época (usando a v0.2, eu acho), Julia era uma linguagem realmente nova. Mas decidi aceitar as dificuldades e tentar codificar meus algoritmos com ela. De qualquer forma, era apenas um pequeno projeto pessoal para aprender mais sobre órbitas. Eu enfrentei muitos bugs e até tive que usar o branch master (pre-v0.3) devido a alguns problemas e funcionalidades que não existiam, mas foi divertido 🙂

Depois de alguns anos (e várias re-escritas devido a mudanças na linguagem), a v0.4 da linguagem Julia foi lançada. Neste momento, dada a quantidade de código que eu tinha e o estado da linguagem, comecei a ver que esse monte de algoritmos poderia, de fato, ser utilizado para auxiliar minhas atividades no INPE. Por isso, decidi criar um pacote privado, chamado SatToolbox.jl, para organizar tudo o que fiz.

Depois de algum tempo, esse pequeno projeto pessoal acabou se tornando o núcleo de um simulador do conceito operacional de missões espaciais chamado Forplan, que está sendo desenvolvido no Centro de Projeto Integrado de Missões Espaciais (CPRIME) do INPE. Dado o bom feedback que recebi, decidi renomear o pacote para SatelliteToolbox.jl e lançá-lo como um pacote oficial da linguagem Julia em março de 2018.

Neste post, gostaria de descrever brevemente o SatelliteToolbox.jl e como ele pode ser utilizado para algumas análises relativas a missões espaciais. Todos os recursos disponíveis podem ser vistos na documentação (disponível apenas em inglês). Uma lista breve dos algoritmos implementados até o momento deste post (na v0.5.0) é:

  • Modelos atmosféricos terrestre:
  • Modelos do campo geomagnético:
  • Índices espaciais:
    • Capacidade de se obter automaticamente diversos índices espaciais, como F10.7, Ap, Kp, etc.
  • Funções para realizar análises gerais relacionadas com órbitas, como conversão de anomalias, cálculo de perturbações, etc.
  • Propagadores orbitais:
    • Two body;
    • J2; e
    • SGP4/SDP4.
  • Funções para converter entre sistemas de referência ECI e ECEF:
    • Toda a teoria IAU-76/FK5 é suportada. Então, a conversão entre quaisquer sistemas de referência a seguir está disponível:
      • ITRF: International Terrestrial Reference Frame;
      • PEF: Pseudo-Earth Fixed reference frame;
      • MOD: Mean-Of-Date reference frame;
      • TOD: True-Of-Date reference frame;
      • GCRF: Geocentric Celestial Reference Frame;
      • J2000: J2000 reference frame;
      • TEME: True Equator, Mean Equinox reference frame.
  • Funções para converter entre referências Geocêntricas e Geodésicas (WGS-84).

A seguir, eu forneço alguns exemplos de como o SatelliteToolbox.jl pode ser utilizado para analisar missões espaciais.

Instalação

A primeira coisa a se fazer (assumindo que você já instalou Julia, que pode ser obtida aqui) é instalar o pacote. Isso pode ser feito digitando:

Para carregar o pacote, que deve ser feito toda a vez que Julia é iniciada, digite:

Informação

Aqui e a seguir, o comando que você deve digitar é o que está apresentado depois de julia>, como você verá no Julia REPL. Todo o resto é o que você deve ver impresso na tela.

Exemplos

Agora, mostrarei algumas análises que podem ser feitas com as funções que estão disponíveis.

Requisitos necessários

Para manter este post pequeno, assumirei que você possui conhecimento sobre linguagem Julia e um pouco de experiência em assuntos relacionados a satélites e órbitas.

Ano Novo na ISS

Vamos ver como podemos calcular onde os astronautas a bordo da ISS estavam durante o Ano Novo em Greenwich! A primeira coisa que devemos fazer é obter a informação sobre a órbita da ISS. Nesse caso, devemos obter o TLE (Two-Line Element), que é um formato de dados consistindo em duas linhas com 70 carácteres cada e que contém toda a informação relacionada à órbita. O TLE a seguir foi obtido do site Celestrak em 4 de Janeiro de 2019, às 12:25 (horário de Brasília).

Esse TLE deve ser carregado em uma variável em Julia. Existem inúmeros métodos para se fazer isso utilizando o SatelliteToolbox.jl. Aqui, utilizaremos uma string especial:

Esse código carrega o primeiro TLE especificado dentro da string entre tle"""...""" para a variável iss_tle.

Agora, devemos inicializar um propagador orbital utilizando o TLE carregado. Nesse caso, vamos usar o SGP4:

A variável orbp agora possui a estrutura do propagador orbital do tipo SGP4 com a órbita especificada pelo TLE iss_tle. Esse TLE foi gerado no dia Juliano 2458487.75253 (2019-01-04 06:03:38.592 +0000). Portanto, devemos retro-propagar a órbita para o instante desejado 2019-01-01 00:00:00.000 +0000 (Ano Novo em Greenwich). Isso pode ser feito através da função propagate!. No entanto, essa função requer o número de segundos a partir da época em que os elementos orbitais foram definidos (nesse caso, a época do TLE) para o qual você deseja propagar a órbita. O código a seguir pode ser utilizado para fazer isso:

A função propagate! retorna três valores. O primeiro o são os elementos orbitais osculantes no instante da propagação, o segundo r_teme é o vetor posição e o último v_teme é o vetor velocidade. Esses vetores são representados no mesmo sistema de referência utilizado para descrever os elementos orbitais quando o propagador foi inicializado. Como estamos utilizando o TLE, então esses vetores estão representados no sistema de referência TEME (True Equator, Mean Equinox).

TEME é um sistema de referência inercial (ECI, Earth-Centered Inertial). Então, devemos converter o vetor posição para um sistema de referência fixo à Terra (ECEF, Earth-Centered, Earth-Fixed) para que possamos calcular a posição da ISS (latitude, longitude e altitude) no instante desejado. O SatelliteToolbox.jl possui toda a teoria IAU-76/FK5 relacionada com a conversão entre sistemas de referência. Neste exemplo, converteremos o TEME para o International Terrestrial Reference Frame (ITRF) para um cálculo mais preciso. Para esse tipo de conversão, deve-se obter os dados da orientação da Terra (EOP, Earth Orientation Data) que são fornecidos pelo IERS. O SatelliteToolbox.jl pode facilmente carregar e usar esses dados através do comando:

Aviso

Os dados EOP são medidos e publicados pelo IERS. Isso significa que o sistema de referência ITRF não pode ser utilizado aqui para prever qual será a posição da ISS, uma vez que os dados EOP não estarão disponíveis. Nesse caso, se um cálculo preciso não é necessário, então pode-se utilizar o sistema de referência Pseudo-Earth Fixed (PEF). A conversão entre o TEME e o PEF utilizando a teoria IAU-76/FK5 não requer dados externos.

A DCM (matriz de cosenos diretores) que gira o TEME para se alinhar ao ITRF é calculada por:

Então, o vetor posição representado no ITRF é:

Finalmente, considerando o elipsóide de referência WGS-84, a latitude, longitude e altitude da ISS durante o Ano Novo em Greenwich pode ser obtida pela função ECEFtoGeodetic conforme mostrado a seguir:

ou seja, latitude 46.189° S, longitude 112.319° W e altitude 419.859 km. Isso está de acordo com o histórico disponível no website I.S.S. Tracker:

Informação

A pequena diferença pode ser explicada pelo TLE que utilizamos para calcular a posição da ISS e o que é utilizado pelo histórico do I.S.S. Tracker. Uma vez que o TLE possui erros e a propagação orbital os aumenta, é melhor obter um TLE gerado o mais próximo possível do instante desejado. Isso não foi buscado aqui.

Perfil de densidade atmosférica

Neste segundo exemplo, usaremos as funções disponíveis no SatelliteToolbox.jl para calcular o perfil de densidade atmosférica. Existem muitos modelos disponíveis na literatura. O SatelliteToolbox.jl implementa quatro deles: o modelo atmosférico exponencial, o Jacchia-Roberts 1971, o Jacchia-Bowman 2008 e o NRLMSISE-00. Todos eles, exceto o primeiro, requerem como entrada alguns índices espaciais, como o F10.7, que mede a atividade do Sol, e o Ap, que mede a atividade geomagnética. O SatelliteToolbox.jl está preparado para obter todos os arquivos necessários da Internet a fim de que esses índices possam ser facilmente calculados. Isso pode ser feito por:

Informação

A keyword wdcfiles_newest_year não é necessária, mas foi utilizada aqui para evitar um erro dado que o arquivo kp2019.wdc não estava disponível quando esse tutorial foi escrito. Para maiores informações, veja a documentação.

Vamos calcular o perfil de densidade atmosférica de 100 km a 1000 km (passo de 1 km), utilizando todos os quatro modelos, em 2018-11-1 00: 00: 00 +0000 sobre a cidade de São José dos Campos, SP, Brasil (Latitude 23.2237 ° S, longitude 45,9009 ° W).

O modelo atmosférico exponencial, o mais simples, não depende nem dos índices espaciais nem da localização, apenas da altitude. Assim, o perfil atmosférico é calculado por:

Cada elemento é a densidade atmosférica [kg / m³] relacionada a uma altitude.

Informação

Aqui, utilizamos o operador de broadcast . para calcular a densidade em todo o intervalo de altitudes utilizando apenas uma linha de código. Para maiores informações, veja a documentação.

Para o modelo Jacchia-Robert 2008, devemos especificar a latitude geodésica [rad], a longitude [rad] e a altitude [m]. Observe que, como já inicializamos os índices espaciais, todas as informações necessárias serão obtidas automaticamente:

Cada elemento é uma instância da estrutura JB2008_Output que contém a densidade dos constituintes da atmosfera em [kg / m³] relacionadas a uma altitude.

O modelo NRLMSISE-00 requer as mesmas informações, mas em uma ordem diferente. Mais uma vez, como já inicializamos os índices espaciais, todas as informações necessárias são obtidas automaticamente:

Cada elemento é uma instância da estrutura NRLMSISE00_Output que contém a densidade dos constituintes da atmosfera em [kg / m³] relacionadas a uma altitude.

NRLMSISE-00 online

Os valores calculados aqui podem ser comparados à versão online do modelo NRLMSISE-00. Observe, no entanto, que a versão online no momento em que este tutorial foi escrito permite apenas o cálculo entre 1960/02/14 e 2018/03/17.

O modelo Jacchia-Roberts 1971 não possui suporte à obtenção automática dos índices espaciais ainda. Portanto, precisaremos fazer isso manualmente. São necessários três índices: o F10.7 diário, a média do F10.7 (janela de 81 dias, centrada no tempo de entrada) e o índice geomagnético Kp (com um atraso de 3 horas). Essa informação pode ser obtida por:

Assim, o perfil atmosférico computado por JR1971 é obtido por:

Cada elemento é uma instância da estrutura JR1971_Output que contém a densidade dos constituintes da atmosfera em [kg / m³] relacionadas a uma altitude.

Finalmente, utilizando o pacote PyPlot.jl, os perfis atmosféricos (altitude vs. densidade) em escala semi-log podem ser plotados usando:

que leva a:

Atmospheric density: 2018-11-01 00:00:00+0000

Para mais informações sobre as muitas opções para calcular a densidade atmosférica, consulte a documentação.

Atenção

Se a data desejada não estiver disponível nos arquivos dos índices espaciais, um erro será gerado ao calcular a densidade atmosférica.

Conclusão

Espero que este tutorial tenha ajudado você a entender um pouco de como o pacote SatelliteToolbox.jl pode ser usado para realizar análises relacionadas a satélites e órbitas. Se você tiver alguma dúvida, por favor, fique à vontade para deixar um comentário abaixo!

Atualização 2018-01-08: Correção do algoritmo que plota os dados, pois a variável h não estava corretamente definida (obrigado Bernard_GODARD).

The Satellite Toolbox for Julia

Hi!

In this post, I would like to introduce the SatelliteToolbox.jl, which is a package for Julia language with many options to analyze space missions. It is used on a daily basis at the Brazilian National Institute for Space Research (INPE). First, it is presented a brief history about the package, and then I show some interesting analysis that can be done with it.

History

In 2013, I joined INPE as a Junior Space Systems Engineer. I was assigned to the space systems division in which I had to work with the attitude and orbit control subsystem (AOCS). Since I had only an intermediate knowledge about orbits, I decided to go into this subject by coding algorithms and comparing the results to the INPE heritage and to the literature in my spare time.

The very first thing was to select the language! In my Ph.D., I used MATLAB to simulate an inertial navigation system, but the Monte Carlo simulations were so slow that I had to rewrite many parts in C using CMEX. On the other hand, in my post-doctoral research where I studied estimation in distributed, non-linear systems, I decided to go with FORTRAN (using the 2008 standards to have at least a readable code…) so that the execution speed is not an issue. Yeah, the performance was very good, but it took me too much time to code. Then I heard about a new language that was promising the best of both worlds: something that resembles an interpreted language with speed of a compiled one! And that’s how I met Julia.

By that time (using v0.2 I think), Julia was a really new language. But I decided to accept the rough edges and try to code my algorithms using it. Anyway, it was just a personal side project to learn more about orbits. I did face many bugs, I had to use master (pre-v0.3) due to some bugs and missing features, but it was fun 🙂

After some years (and huge rewrites due to breaking changes), Julia released its v0.4. In this time, given the amount of code I had and the state of the language, I started to see that this bunch of algorithms can indeed be used for something at INPE to help in my activities. Hence, I decided to create a private package, which was called SatToolbox.jl, to organize everything I have done.

After some time, this little side project turned out to be a very important core for a simulator, called Forplan, of the mission operational concept we started to develop at INPE’s Space Mission Integrated Design Center (CPRIME). Given the good feedback I received, I decided to rename the toolbox to SatelliteToolbox.jl and it was released as an official Julia package in March 2018.

In this post, I would like to briefly describe SatelliteToolbox.jl and how it can be used to analyze space missions. The entire set of the features can be seen in the documentation. A brief, non-exhaustive list of the algorithms implemented by the time of this post (in v0.5.0) is:

  • Earth atmospheric models:
  • Earth geomagnetic field models:
  • Space indices:
    • Capability to automatically fetch many space indices, e.g. F10.7, Ap, Kp, etc.
  • Functions to perform general analysis related to orbits, e.g. converting anomalies, computing perturbations, etc.
  • Orbit propagators:
    • Two body;
    • J2; and
    • SGP4/SDP4.
  • Functions to convert between ECI and ECEF references frames:
    • All the IAU-76/FK5 theory is supported. Hence, the conversion between any of the following frames is available:
      • ITRF: International Terrestrial Reference Frame;
      • PEF: Pseudo-Earth Fixed reference frame;
      • MOD: Mean-Of-Date reference frame;
      • TOD: True-Of-Date reference frame;
      • GCRF: Geocentric Celestial Reference Frame;
      • J2000: J2000 reference frame;
      • TEME: True Equator, Mean Equinox reference frame.
  • Functions to convert between Geocentric and Geodetic (WGS-84) references.

In the following, I provide a few examples of how SatelliteToolbox.jl can be used to analyze space missions.

Installation

The very first thing (provided that you have Julia installed, which can be obtained here) is to install the package. This can be done by typing:

To load the package, which must be done every time Julia is restarted, type:

Information

Here and in the following, the command you have to type is what is presented after julia>, as you will see in the Julia REPL. Everything else is what you should see printed on the screen.

Examples

Now, I will show some analysis that can be done with the functions that are available.

Necessary background

To keep this post short, I will assume that you have knowledge about Julia and a little background in satellite and orbits.

New Year in ISS

Let’s see how can we calculate where the astronauts onboard the ISS were during the New Year in Greenwich! The first thing we need to do is to get information about the ISS orbit. In this case, we must obtain the TLE (Two-Line Element) file, which is a data format consisting of two lines with 70 characters each that contain all information related to the orbit. The following TLE was obtained from Celestrak at January 4, 2019, 12:25 +0200.

This TLE must be loaded into a variable inside Julia. There are a number of methods to do this using SatelliteToolbox.jl. Here, we will use a special string type:

This code loads the first TLE specified inside the string enclosed by tle"""...""" into the variable iss_tle.

Now, we have to initialize an orbit propagator using the loaded TLE. In this case, we will use the SGP4:

The variable orbp now holds the orbit propagator structure of type SGP4 with the orbit specified by the TLE iss_tle. This TLE was generated at the Julian day 2458487.75253 (2019-01-04 06:03:38.592 +0000). Thus, we have to backpropagate the orbit to the desired instant 2019-01-01 00:00:00.000 +0000 (New Year in Greenwich). This can be accomplished by the function propagate!. However, this function requires a number of seconds from the orbital elements epoch (in this case, the TLE epoch) to which you want to propagate the orbit. Thus, the following code can be used to accomplish this:

The function propagate! returns three values. The first one o is the osculating orbit elements at the propagation instant, the second one r_teme is the position vector, and the last one v_teme is the velocity vector. Those vectors are represented on the same reference frame that was used to describe the orbit elements when the propagator was initialized. Since we are using the TLE, those vectors are represented in the TEME (True Equator, Mean Equinox) reference frame.

TEME is an Earth-Centered Inertial (ECI) reference frame. Hence, we must convert the position vector to an Earth-Centered, Earth-Fixed (ECEF) reference frame so that we can compute what was the ISS position (latitude, longitude, and altitude) at the desired instant. SatelliteToolbox.jl has the entire IAU-76/FK5 theory related to the conversion between reference frames. For this example, we will convert TEME into the International Terrestrial Reference Frame (ITRF) for more accurate computation. This kind of conversion requires the Earth Orientation Data (EOP) that is provided by IERS. SatelliteToolbox.jl can easily load and use this data as follows:

Warning

The EOP data is measured and published by IERS. This means that the ITRF cannot be used here to predict what will be the ISS position in the future since the EOP data will not be available. In this scenario, if high accuracy is not needed, then the Pseudo-Earth Fixed (PEF) reference frame can be used. The conversion between TEME and PEF using IAU-76/FK5 does not require external data.

The DCM (Direction Cosine Matrix) that rotates TEME into alignment with ITRF is computed by:

Thus, the position vector represented in ITRF is:

Finally, considering the WGS-84 reference ellipsoid, the latitude, longitude, and altitude of the ISS during the New Year in Greenwich can be obtained by the function ECEFtoGeodetic as follows:

i.e., latitude 46.189° S, longitude 112.319° W, altitude 419.859 km. This is in agreement with the historical information on I.S.S. Tracker website:

Information

The minor difference can be explained by the TLE we used to computed the ISS position and the one used by the historical data of I.S.S. Tracker. Since a TLE has errors and the propagation increases them, it is desirable to obtain a TLE as closest as possible to the desired instant, which was not pursued here.

Atmospheric density profile

In this second example, we will use the built-in functions of SatelliteToolbox.jl to compute the atmospheric density profile. There are many models available in the literature. SatelliteToolbox.jl implements four of them: Exponential atmospheric model, Jacchia-Roberts 1971, Jacchia-Bowman 2008, and NRLMSISE-00. All of them but the former requires as input some space indices, like the F10.7 that measures the Sun activity and the Ap that measures the geomagnetic activity. SatelliteToolbox.jl is prepared to download all the required files from the Internet so that those indices can be easily obtained. This can be accomplished by:

Information

The keyword wdcfiles_newest_year is not necessary, but it is added to avoid an error since the kp2019.wdc file was not available by the time this tutorial was written. For more information, see the documentation.

We will compute the atmospheric density profile from 100 km to 1000 km (steps of 1 km) using all the four models at 2018-11-1 00:00:00+0000 over São José dos Campos, SP, Brazil (Latitude 23.2237° S, Longitude 45.9009° W).

The exponential atmospheric model, the simpler one, depends neither on space indices nor on the location, only on the altitude. Thus, the atmospheric profile is computed by:

Each element is the atmospheric density [kg/m³] related to one altitude.

Information

Here we use the Julia broadcast operator . to compute for the entire altitude interval using only one line. For more information, see the documentation.

For the Jacchia-Robert 2008, we must specify the geodetic latitude [rad], longitude [rad], and altitude. Notice that, since we have already initialized the space indices, all the required information will be gathered automatically:

Each element is an instance of the structure JB2008_Output that contains the density of the atmospheric species in [kg/m³] related to one altitude.

The NRLMSISE-00 requires the same information, but in a different order. One more time, since we have already initialized the space indices, all the required information is fetched automatically:

Each element is an instance of the structure NRLMSISE00_Output that contains the density of the atmospheric species in [kg/m³] related to one altitude.

Online NRLMSISE-00

The values computed here can be compared to the online NRLMSISE-00 model. Notice, however, that the online version by the time this tutorial was written only allows computation between 1960/02/14 and 2018/03/17.

The Jacchia-Roberts 1971 model does not support fetching the space indices automatically yet. Hence, we will need to do this manually. It requires three indices: the daily F10.7, the averaged F10.7 (81-day window, centered on input time), and the Kp geomagnetic index (with a delay of 3 hours). That information can be fetched by:

Thus, the atmospheric profile computed by JR1971 is obtained by:

Each element is an instance of the structure JR1971_Output that contains the density of the atmospheric species in [kg/m³] related to one altitude.

Finally, using the PyPlot.jl package, the atmospheric profiles (altitude vs. density) in semi-log scale can be plotted using:

which leads to:

Atmospheric density: 2018-11-01 00:00:00+0000

For more information about the many options to compute the atmospheric density, please see the documentation.

Warning

If the desired date is not available on the space indices files, then an error will be thrown when computing the atmospheric density.

Conclusion

I hope that this tutorial has helped you to understand a little bit how SatelliteToolbox.jl can be used to perform analysis related to satellites and orbits. If you have any question, please, feel free to leave a comment below!

Update 2018-01-08: Fix algorithm that plots the data, because the variable h was not defined (thanks Bernard_GODARD).