Análise da órbita de um veículo espacial

Prévia do material em texto

UNIVERSIDADE FEDERAL DO ABC 
 
 
 
 
 
 
RELATÓRIO 1 DA DISCIPLINA LABORATÓRIO DE NAVEGAÇÃO, 
GUIAGEM E CONTROLE 
 
 
 
 
Caio Paradizo Benedetti, RA: 11094714 
João Paulo Silva Souza, RA: 21054015 
 
 
 
 
Analisar a órbita que um veículo espacial descreve em torno da 
terra pelo problema de dois corpos 
 
 
 
 
 
 
 
 
São Bernardo - SP 
2019 
 
 
1. Introdução 
Uma trajetória é definida como o caminho que um objeto percorre pelo espaço. De 
forma a lançar um veículo espacial, o veículo lançador busca seguir uma trajetória de 
máxima eficiência. Uma órbita é onde um objeto reside no espaço. Robert Giffen [1] faz um 
paralelo entre órbitas e pistas de corridas, em que o corpo viaja ao redor de um planeta ou 
outro corpo celestial, ou seja, pode ser entendido como um movimento padrão regular e 
repetido que um objeto no espaço percorre ao redor de outro. 
A princípio, existem infinitas possíveis órbitas, sendo necessário escolher a que 
melhor atenda os objetivos de uma determinada missão. Tamanho da órbita, formato e 
orientação são fatores comumente utilizados para determinar a órbita a ser utilizada. Fácil 
de se pensar que a altura da órbita também interfere em sua escolha, uma vez que requer 
maior energia para o lançamento do corpo, um maior veículo lançador, resultando em 
maiores gastos. Por outro lado, uma órbita mais alta permite a visualização de uma maior 
área. 
Uma vez sabido o que é uma órbita, convém ressaltar como podemos realizar a 
análise de seu movimento. Motion Analysis Process (MAP) é um sistema criado para a 
análise e entendimento de todos os tipos de movimento de um objeto pelo espaço. Consiste 
na definição de um sistema de coordenadas e um ponto de referência, relação de seu 
movimento e as forças envolvidas, assumir modelos simplificados, definição de condições 
iniciais, realizar análise de erros e testar o modelo. Este raciocínio será utilizado para a 
modelagem do problema de dois corpos, a ser introduzido no capítulo de Fundamentação 
Teórica. 
O objetivo deste trabalho é determinar numericamente dada órbita, utilizando o 
software MatLab, baseando-se em simplificações e assumindo o problema restrito de dois 
corpos. Além disso, de forma a comparar e verificar os resultados, a órbita também será 
representada a partir da utilização do software STK.O integrador numérico ODE45 será 
utilizado com o objetivo de solucionar a EDO de segunda ordem obtida a partir do problema 
de dois corpos. 
 
2. Fundamentação Teórica 
No passado o ser humano percebeu que alguns astros descrevem um movimento 
regular, resultando em noção do tempo e épocas do ano. A primeira conclusão foi o modelo 
Geocêntrico, em que acreditava-se que o Sol e os outros planetas giravam ao redor da 
terra. Posteriormente, Copérnico apresentou o modelo Heliocêntrico, em que o Sol estava 
no centro do universo e os planetas descreviam órbitas ao seu redor. Anos mais tarde, 
Johanes Kepler enunciou 3 leis que descrevem o movimento planetário. 
Primeira Lei de Kepler → ​As órbitas dos planetas são elipses com o Sol no foco. 
Basicamente, à partir desta lei Kepler definiu que as órbitas não eram circunferências como 
acreditava-se, mas elipses. 
Segunda Lei de Kepler → ​Em sua segunda lei, Kepler estudou a velocidade com 
que um planeta orbita em torno do Sol, assumindo uma relação de área com os períodos. 
De acordo com Halliday [5], uma linha unindo um planeta ao Sol varre áreas iguais em 
períodos de tempos iguais. 
Terceira Lei de Kepler → Lei dos períodos. A partir dessa lei, Kepler apresenta uma 
relação diretamente proporcional entre o período de revolução de um planeta ao redor do 
 
 
Sol e o raio médio da órbita do planeta. Segundo ele, “Os quadrados dos períodos de 
revolução dos planetas ao redor do Sol são diretamente proporcionais aos cubos dos raios 
médios de suas órbitas”. 
A partir de suas 3 leis que descrevem o movimento dos planetas, Kepler foi capaz de 
dar uma nova ênfase em como entender e quantificar as causas físicas de movimento, se 
tornando um dos precursores da astrodinâmica. Seu trabalho, somado a de outros grandes, 
a citar: Brahe e Galileo, permitiu grandes avanços no estudo de dinâmica orbital. 
Isaac Newton, filósofo e matemático, a partir do trabalho de Galileo em dinâmica, e 
“apoiado no ombro de gigantes” publicou suas 3 leis do movimento e a gravitação universal. 
A partir de suas leis, nos tornamos capazes de prever não apenas o movimento da Terra, 
mas também de cometas, luas e etc. 
Primeira Lei de Newton → Na ausência de forças, um corpo em repouso tende a 
ficar em repouso e um corpo em movimento tende a permanecer em movimento retilíneo e 
uniforme. É sabido que além disso, um corpo em repouso possui resistência ao movimento, 
assim como um objeto em movimento resiste a uma mudança de velocidade ou direção. É 
importante observar que apesar de ambos resistirem ao movimento, a quantidade de 
resistência diverge devido a existência do conceito de Momento. Existem 2 tipos de 
momentos: momento linear e momento angular. 
Momento Linear → Basicamente trata-se de uma grandeza capaz de medir a 
quantidade de movimento presente em um determinado corpo. Pode ser calculado como: 
 
 mVp = (1) 
 
Onde p é o vetor momento linear (kgm/s), m é a massa (kg) e V é o vetor velocidade 
(m/s). 
Momento Angular ​→ O momento angular, por sua vez, trata-se da resistência de 
um objeto em rotação de alterar a velocidade de rotação ou direção. Pode ser calculado 
como: 
 IΩH = (2) 
 
Sendo, H o vetor momento angular (kgm²/s), I é o momento de Inércia (kg m²) e o Ω 
vetor velocidade angular (rad/s). 
Uma vez definido os conceitos de momento, seja angular ou linear, podemos 
determinar a segunda lei de Newton, a citar: 
Segunda Lei de Newton → “A força resultante que atua sobre um corpo é 
proporcional ao produto da massa pela aceleração por ele adquirida.” ou ainda, a taxa 
temporal de mudança de momento de um determinado objeto relaciona-se com a força 
aplicada. Isso quer dizer que, se desejamos alterar rapidamente o momento de um objeto, 
devemos aplicar uma grande força. Caso queiramos alterar o momento de um objeto 
aplicando uma força menor, devemos fazê-la por um tempo maior. Matematicamente, 
temos: 
F maΔt
Δp = variação tempo
variação momento → ∑ = (3) 
 
Onde F é a somatória das forças (N), m é a massa (kg) e a é a aceleração (m/s²).∑ 
 
 
Terceira Lei de Newton ​→ Ação e Reação. Newton nos diz que: “Quando um objeto 
A exerce uma força sobre um outro objeto B, este outro objeto B vai exercer uma força de 
mesma intensidade, mesma direção e sentido contrário sobre o objeto A”. 
Lei da Gravitação Universal → “A força gravitacional entre dois corpos é 
diretamente proporcional ao produto das duas massas e inversamente proporcional ao 
quadrado da distância entre eles”. Matematicamente, podemos representar como: 
 
 F = R²
Gm1m2 (4) 
 
Para F a força devido a gravidade (N), G a constante universal gravitacional 
(6,67e-11Nm²/kg²), m1 e m2 a massa dos dois corpos (kg) e R a distância entre os dois 
corpos (m). 
 
Leis de Conservação:Sabe-se que as leis de Newton, apresentadas anteriormente, conduzem ao 
princípio da conservação de energia e de momento (angular e linear). As leis da 
conservação nos mostram que que um sistema ao passar por uma transformação apresenta 
as mesmas propriedades antes, durante e depois do processo. 
Conservação de Momento ​→ A terceira lei de Newton implica que o momento em 
um determinado sistema é conservado. Importante ressaltar que a resultante das forças 
externas deve ser nula, ou até mesmo não existirem forças externas atuando sobre o 
sistema. 
Conservação de Energia ​→ Energia pode ter diversas formas, a citar: mecânica, 
nuclear, química e elétrica. Para o estudo do movimento convém enfatizar o estudo de 
energia mecânica. A energia total mecânica surge a partir da posição e movimento de um 
objeto e é composta de Energia Cinética e Energia Potencial. Para entender melhor, 
convém ressaltar que o campo gravitacional é conservativo. Dessa forma a soma de 
energias em um campo conservativo é constante. Matematicamente, 
 
 KE PEE = + (5) 
 
Onde E é a energia mecânica total (kg m²/s²), PE a energia potencial (kg m²/s²) e KE 
a energia cinética (kg m²/s²). A energia potencial relaciona-se com a posição de um 
determinado objeto e é dada por: 
E ghP = m ​(6) 
 
Onde m é a massa (kg), g é a aceleração da gravidade (m/s²) e h é a altura do ponto 
de referência (m). 
Quando estudamos a energia potencial de um satélite sobre a terra, não podemos 
assumir que a gravidade é constante, tampouco utilizar a superfície da terra como um ponto 
de referência. Dessa forma, devemos derivar uma equação que relaciona a quantidade de 
trabalho realizado para mover um satélite desde o centro da terra até a localização da 
órbita. Essa derivação pode ser encontrada em [1] ​apêndice C.10. Aqui assumimos que a 
energia potencial orbital pode ser dada por 
 
 
 
E μ/RP = − m (7) 
Para PE a energia potencial da espaçonave (kg km²/s²), m a massa da espaçonave 
(kg), o parâmetro gravitacional (3,986e5 km³/s²) e R a distância da espaçonave até o μ 
centro da Terra (km). 
A energia cinética, por outro lado, está relacionada com a velocidade e a massa de 
um determinado corpo. Pode ser escrita como: 
 
E mV ²K = 2
1 (8) 
KE é a energia cinética (kg km²/s²), m a massa (kg) e V a velocidade (km/s). Sendo 
assim, a Energia total de um sistema pode ser escrito como: 
 
 mV ² μ/R E = 2
1 − m (9) 
 
Geometria Orbital: 
Uma vez que estamos interessados em estudar órbitas de espaçonaves, que são 
elípticas, convém apresentar os principais parâmetros geométricos da órbita, de forma a 
facilitar o entendimento das equações que serão apresentadas. 
 
R é o raio desde o foco da elipse, para este caso desde o centro da terra até a espaçonave. 
Rp é o raio do perigeu. Ra é o raio do apogeu. Abaixo seguem algumas relações 
importantes: 
 
a Ra Rp2 = + (10) c Ra Rp2 = − (11) 2c/2ae = (12) 
 
Uma vez introduzida as 4 leis de Newton, conceito de conservação e apresentados 
os parâmetros de uma órbita elíptica, somos capazes de entender o problema restrito de 
dois corpos. 
 
 
 
 
 
 
Problema Restrito de Dois Corpos: 
Inicialmente, de forma a sermos capazes de aplicar as leis de Newton, é importante 
ressaltar que devemos expressá-las em um referencial inercial, ou seja, que não possui 
aceleração. Quando estudamos equações de movimento, o sistema de coordenada 
escolhido é de crucial importância, uma vez que a depender da escolha podemos facilitar a 
resolução ou até mesmo torná-la impossível. Para o estudo de aeronaves orbitando a terra, 
escolhe-se o sistema de coordenada Geocêntrico-Equatorial. Suas principais 
características são: 
● Origem: Centro da Terra; 
● Plano fundamental: Equador Terrestre. Perpendicular ao plano: Direção polo 
norte; 
● Direção Principal: Equinócio vernal; 
● Terceira direção: Regra da mão direita. 
 
Uma vez escolhido o sistema de coordenadas, pode-se aplicar a segunda lei de 
Newton e analisar as forças que agem no sistema, como por exemplo: Empuxo, arrasto, 
gravidade, força de um terceiro corpo, entre outros. Observa-se que caso optemos por não 
simplificar esse estudo, nos deparamos com infinitas forças que causam perturbação no 
sistema, dificultando a resolução do problema. Nesse contexto que surge o problema 
restrito de dois corpos. Neste caso desprezamos praticamente todas as forças e levamos 
em consideração apenas a gravitacional. A partir dele somos capazes de estudar o 
movimento relativo de um corpo em relação a outro, considerando-se dois pontos de massa 
que interagem entre si apenas devido à atração gravitacional. Sendo assim, a segunda lei 
de Newton se torna: 
ext Fgrav ma∑
 
 
F = = (12) 
Na forma vetorial, pode-se escrever (12) como: 
(13) 
 
Chegando então a equação do problema de dois corpos: 
Onde d²R/dt² é a aceleração da espaçonave (km/s²), é o parâmetro gravitacional, μ 
R é o vetor posição da espaçonave (km) e R é a magnitude da posição da espaçonave 
(km). 
(14) 
 
Verifica-se que a equação (14) ​trata-se de uma equação diferencial vetorial, não 
linear de segunda ordem. Dessa forma, de forma a facilitar o entendimento, abaixo é 
possível observar sua solução analítica disposta unidimensionalmente. 
(15) 
Importante citar que a equação 15​ ​trata-se de uma equação das cônicas. 
 
 
 
Constantes de movimento e período orbital: 
Lembrando que espaçonaves conservam momento e energia, somos capazes de 
descobrir constantes de movimento orbital. 
Energia Mecânica Específica​ → Da conservação de energia, temos: 
 
 E/m ε V ²/2 μ/R V ε = → = − → = √2(μ/R )+ ε (16) 
Sabendo que 
 /(2a)ε = − μ (17) 
 
Somos capazes de determinar a energia específica a partir do semi eixo maior da 
órbita. Além disso, podemos entender o tipo de trajetória a partir dos valores de energia. Ou 
seja, 
Se <0 → órbita circular ou elíptica;ε 
Se =0 → órbita parabólica;ε 
Se >0 → órbita hiperbólica.ε 
 
Além disso, uma vez sabido a energia específica, somos capazes de calcular o 
período orbital. Período orbital pode ser definido como o tempo que uma espaçonave leva 
para dar uma volta completa em sua órbita. Matematicamente, 
 
 2π P = √a³/μ (18) 
 
Formas e Representação da Terra: Geóide Elipsóide: 
Quando desejamos descrever a posição de objetos utilizamos sistemas de 
referência. Na terra, utiliza-se o sistema de referência terrestre ou geodésico. Um elipsóide 
é uma superfície definida matematicamente que se aproxima da forma da terra sobre o 
qual são realizados os cálculos relacionados a suas coordenadas. Normalmente, cada país 
ou grupo de países adotam um elipsóide como referência para os trabalhos geodésicos e 
topográficos que mais se aproxima do geóide na região considerada. Geralmente define-se 
um elipsóide de referência para especificar o raio do semi-eixo equatorial e a relação de 
achatamento. Para este trabalho utilizou-se o elipsóide GRS80. Na tabela a seguir é 
possível observar as principais diferenças presentes entre os elipsóides GRS80 e WGS84. 
 
 
Verifica-se que o WGS84 é 0,21 mm maior que o GRS80. Enquanto o WGS84 é 
utilizado por todo o mundo, o elipsóide GRS80 é utilizado principalmente nos Estados 
Unidose Canadá. WGS84 possui um erro menor que 2 cm do centro de massa, melhor 
resultado que o GRS80. Outra principal diferença entre os dois modelos elipsóides é 
relacionado aos pontos ao norte da américa, onde para o WGS84 define-se de acordo com 
a média das estações ao redor do mundo, enquanto para o GRS80 não muda. 
 
 
 
3. Método Numérico 
O ODE45 é o integrador padrão de EDO’s do software Matlab. É baseado em um 
algoritmo de Dormand and Prince. Utiliza a estratégia FSAL e utiliza métodos de 
Runge-Kutta de quarta e quinta ordem. Possui 6 estágios. 
ODE23 utiliza também métodos de Runge-Kutta, porém de segunda e terceira 
ordem. O 23 no nome indica que dois ​single-steps ​funcionam simultaneamente, sendo que 
estão envolvidos um de terceira ordem e um de segunda ordem. Possui apenas 3 estágios. 
Quando comparado com o ODE23, ODE45 demora mais por step, mas é capaz de 
trabalhar com steps maiores. Para trabalhar com EDO’s que necessitam de resultados mais 
acurados, em geral utiliza-se o ODE45. 
Também existe o ODE113 que se baseia no método de Adams-Bashforth-Moulton 
de ordem variável. Para problemas com pouca tolerância a erros ou problemas 
computacionalmente pesados. Ou ainda pode ser usado o ODE15s para situações em que 
o ODE45 se torna muito lento, provavelmente porque a equação diferencial deve ser de um 
problema stiff. Além desse, podemos usar o ode23s, ode23t ou ode23tb. 
De maneira geral, entre os operadores citados é preferível usar o ODE45 como 
primeira escolha para resolver uma equação diferencial, uma vez que ele consegue 
minimizar o erro e manter um baixo tempo de programação. 
No nosso caso ao rodarmos o programa com o ODE23 observamos um menor peso 
computacional do programa e um menor tempo de programação comparado a todos os 
outros, mas o ODE45 acaba sendo a melhor opção devido seu erro ser o menor entre os 
integradores e o seu tempo de programação ser razoável. 
 
Reduzindo EDO para equação de primeira ordem: 
Vimos anteriormente que a partir do problema restrito de dois corpos, podemos 
escrever a equação do movimento como: 
 ​ (19) 
Por se tratar de uma EDO de segunda ordem, devemos reduzir para uma equação 
de primeira ordem, para que seja possível aplicarmos o integrador ODE45. Observa-se que 
por R se tratar de um vetor, com coordenadas x, y e z, podemos abrir a Eq. 14 em outras 3 
Eq. 
 
Assumindo que 
 
 
 
Conclui-se que 
, e 
Dessa forma somos capazes de reduzir a ordem de uma EDO e obter a solução desse 
sistema de EDO’s a partir da função ODE45. 
 
 
 
Formulação MatLab: 
No código DoisCorpos.m foi feita uma função para ser chamada no código 
IntegradorODE45. Essa função tem como objetivo fazer uma transformação de variáveis 
para transformar a nossa equação (19), que é uma EDO de segunda ordem, em uma EDO 
de primeira ordem, já que o método de Runge-Kutta só pode ser usado para EDOs de 
primeira ordem no máximo. Como solução foi feito um vetor vet 6x1 e um vetor w 6x1 para 
armazenar as equações remodeladas. 
No código IntegradorODE45 foi inicialmente calculado o valor do período orbital 
partindo das informações iniciais fornecidas e usando a equação (18). Depois, inserimos na 
posição do tempo final o valor do período na função ODE45, e também foram inseridos os 
dados iniciais nas posições adequadas dentro do operador. Como mencionado antes, o 
operador ODE45 faz o procedimento de integração numérica, e com o auxílio dos dados 
iniciais é possível plotar a órbita que o satélite desenvolve. 
Após o uso do ODE45, foi inserida uma esfera aproximada representando a terra, foi 
usado o GRS80 como parâmetros para isso. Depois que o globo GRS80 foi construído e 
posicionado, foi plotado o globo junto com a órbita anteriormente criada, mas agora ambos 
no mesmo eixo de referência, no caso o da terra. 
Também foi analisada a órbita desenvolvida pelo satélite para valores de período 
maiores, no caso para 3 e 10 períodos, e isso foi feito com o intuito de medir o erro dentro 
do integrador numérico. 
Abaixo segue o fluxograma criado para o desenvolvimento do código. 
 
 
Os códigos desenvolvidos serão dispostos ao final deste trabalho, no apêndice 1. 
 
4. Resultados Numéricos 
Utilizando os conceitos definidos anteriormente, foi possível o cálculo do período 
orbital, no caso foi de aproximadamente 86502s, e como o dia sideral terrestre é de 
aproximadamente 86164s, podemos aproximar a órbita encontrada para um geossíncrona 
ou até mesmo geoestacionária, outros pontos que corroboram com essa aproximação serão 
apresentados nessa seção, por exemplo a leve inclinação e excentricidade calculadas. Uma 
vez implementado o código e com o período definido, foi possível observar que quanto 
maior o valor de input do período menor o raio orbital, e isso implica que o satélite pode 
acabar entrando na atmosfera ou acabar se chocando com o solo. 
 
 
Num primeiro momento acreditamos que isso era devido a ações de forças externas 
,mas não poderia ser verdade porque estamos considerando um problema restrito de dois 
corpos, ou seja, somente a força gravitacional foi considerada e computada. Portanto, essa 
anomalia pode ser explicada em virtude de erros de arredondamento devido ao próprio 
operador. 
É importante ressaltar que é normal que ocorra erros de arredondamento em 
métodos numéricos, o importante é saber quando são significativos ou não, e se caso forem 
significativos é importante empregar um método adequado para corrigir ou controlar esse 
erro. Como no nosso caso constatamos que em poucos períodos já ocorria a entrada 
aparente do satélite na atmosfera e como a altitude encontrada foi de aproximadamente 
35896 km, ou seja, o satélite estava inicialmente pouco acima de uma órbita 
geoestacionária, julgamos que esse erro numérico representa algo significativo, lembrando 
que estamos falando do problema restrito de dois corpos, e portanto deveria ser controlado. 
Para corrigir esse problema de erro de arredondamento dentro do nosso integrador 
foi adicionado a seguinte linha ao código IntegradorODE45.m: “options = odeset('RelTol', 
0.00001);” [2]. As imagens (1), (2) e (3) mostram, respectivamente, os resultados para 1 
período, 3 períodos e 10 períodos sem considerar o erro. Já a imagem (4) considera erro de 
0,00001 para 10 períodos, repare que fica igual a imagem (1) porque o satélite permaneceu 
na mesma órbita depois de 10 voltas. 
 
(1) (2) 
(3) (4) 
 
Utilizando o matlab, fomos capazes também de calcular algum dos principais 
elementos orbitais, que posteriormente serão utilizados para plotar a órbita no software 
STK. Importante ressaltar que nem sempre podemos definir os elementos orbitais clássicos. 
Por exemplo, quando lidamos com uma órbita circular, neste caso o argumento do perigeu 
ou a anomalia verdadeira não estão definidos. Quando lidamos também com uma órbita 
equatorial, ou seja, seu ângulo de inclinação é 0º ou 180º o nodo ascendente não existe. 
Ao calcular o módulo da excentricidade obtivemos um valor de 0,0024, 
aproximando-se de uma órbita circular, como esperado. Além disso, a inclinação 
encontrada foi de 0,0186 graus, se aproximando de uma órbitaequatorial. Vale ressaltar 
 
 
que, assim como estudado na teoria, ao tentarmos calcular a anomalia verdadeira e o 
argumento do perigeu, encontramos números imaginários, reforçando a ideia de uma órbita 
circula, já que a mesma não possui ambos (não tem perigeu), mas sim argumento de 
latitude. A partir do software STK, foi possível plotar a órbita, restando calcular a altitude. 
Uma vez que a órbita estudada se assemelha com uma órbita circular, foi possível descobrir 
a altitude a partir da diferença entre o raio terrestre e o raio da órbita. Os resultados obtidos 
no STK foram: 
 
 
 
5. Conclusão 
 
A órbita descrita pelo veículo espacial com as condições iniciais propostas é uma 
elipse, semelhante a uma órbita circular-equatorial, uma vez possuidora de excentricidade e 
inclinação de aproximadamente 0. Além disso, devido a aproximação do período sideral 
com o período calculado para essa órbita somado com os dois resultados anteriormente 
citados, podemos aproximar essa órbita para uma geossíncrona ou até mesmo 
geoestacionária. A determinação dessa órbita foi feita considerando o sistema geocêntrico 
inercial para uma melhor comparação visual e para a realização dos cálculos mantendo o 
satélite e a terra no mesmo sistema. A utilização do problema de dois corpos foi utilizada 
com o intuito de simplificar o nosso modelo de tal maneira que obtivemos uma EDO de 
segunda ordem vetorial, para solucionar essa equação fizemos uso dos integradores 
numéricos via Matlab, e o ODE45 se mostrou o melhor para nosso problema. Entretanto, 
nossa solução para a órbita apresentou um erro significativo por parte do integrador 
numérico que foi solucionado. Pela análise dos dados e comparação com a órbita obtida no 
STK confirmamos que o nosso modelo analítico e numérico estão condizentes. 
 
6. Referências 
 
[1] SELLERS, J. J.; ASTORE, W. J.; GRIFFEN, R. B.; LARSON, W. Understanding Space: 
An Introduction to Astronautics. 3. ed. New York: McGraw-Hill, 2000. 
[2]​What Every Engineer Should Know about MATLAB® and Simulink® 
[3] Orbital Mechanics for Engineering Students, Elsevier Aerospace Engineering Series, 
CURTIS, H., 2005. 
[4] Matlab Documentation. 
[5]David Halliday, Robert Resnik, and Denneth S. Krane, Fisica 2, 5th ed. Rio de 
Janeiro: LTC, 2004, vol. 1. 
 
 
 
 
 
7. Apêndice 
 
APENDICE A - Códigos 
 
Código DoisCorpos.m​: 
 
%A ideia foi fazer uma transformacao da nossa eq do problema de DoisCorpos, 
%que é uma EDO de segunda ordem, em uma EDO de primeira ordem para usar o 
%método de Ruge-Kutta através da função ode45. 
 
function [vet] = DoisCorpos (t, w) 
%criando uma funçao envolvendo o vetor vet e o vetor w(qualquer), 
%onde os 3 primeiros termos são referentes a posição e os 3 últimos a 
%velocidade. 
mi = 3.986*10^5; %[(km^3)*(s^-2)], é o parâmetro gravitacional da terra. 
r = (w(1)*w(1) + w(2)*w(2) + w(3)*w(3))^(1/2); %raio(em módulo). 
vet(1, 1)= w(4); %velocidade 
vet(2, 1)= w(5); %velocidade 
vet(3, 1)= w(6); %velocidade 
vet(4, 1)= -(mi/(r)^3)*w(1);%aceleração 
vet(5, 1)= -(mi/(r)^3)*w(2);%aceleração 
vet(6, 1)= -(mi/(r)^3)*w(3);%aceleração 
 
IntegradorODE45.m​: 
 
clear all 
clc 
 
%cálculos: 
mi = 3.986*10^5; %[(km^3)*(s^-2)], é o parâmetro gravitacional da terra. 
si = [3010.33 -42067.38 -0.59];%[km], é a posição inicial. 
vi = [3.07 0.22 0.001]; %[km*(s^-1)], é a velocidade inicial. 
r = norm(si); %[km], é o raio da órbita(em módulo). 
v = norm(vi); %[km*(s^-1)], é o módulo da velocidade. 
epsilon = ((v^2)/2)-(mi/r); %[(km^2)*(s^-2)], é a energia total específica. 
a = -mi/(2*epsilon); %[km], é o semieixo maior da elipse 
periodo = 2*pi*sqrt((a^3)/mi); %[s], é o tempo que um obj. leva para completar um volta na 
órbita. 
% 
 
t = [0:0.1:periodo]; 
 
 
 
options = odeset('RelTol', 0.00001); %erro adotado para evitar que o satélite 
%O vetor T está associado ao vetor t do DoisCorpos que representa o tempo. 
%O vetor W está associado ao vetor w do DoisCorpos que representa as 
%condições iniciais(as 3 primeiras da posição, as 3 últimas da velociade). 
[t,W]= ode45 (@DoisCorpos, t,[3010.33 -42067.38 -0.59 3.07 0.22 0.001], options); 
 
 
%Inserindo o globo GRS80: 
grs80 = referenceEllipsoid('grs80','km'); 
domeRadius = 6378.137;%[km], é o raio da terra. 
domeLat = 0; %[graus], posição de um obj em relação a terra(latitude). 
domeLon = 0; %[graus], posição de um obj em relação a terra(longitude). 
domeAlt = 0; %[km], posição de um obj em relação a terra(altura partindo do centro). 
 
[al,bl,cl] = sphere(20); 
xEast = domeRadius * al; 
yNorth = domeRadius * bl; 
zUp = domeRadius * cl; 
zUp(0<zUp < 0) = 0; 
 
figure('Renderer','opengl') 
surf(xEast, yNorth, zUp,'FaceColor','blue','FaceAlpha',0.5) 
axis equal 
% 
 
%O comando seguinte é usado para plotar duas geometrias na mesma figura, 
%no caso plotamos a órbita junto com o planeta. 
hold on 
 
%Plotando a órbita junto com o globo. 
plot3(W(:,1),W(:,2),W(:,3)); 
xlabel ('x[Km]'); 
ylabel ('y[Km]'); 
zlabel ('z[Km]'); 
grid3 
% 
 
ev = 1/mi*((v^2 - mi/r)*si - (dot(si, vi))*vi); %calculo excentricidade vetor 
e = norm(ev); %modulo excentricidade 
hv = cross (si, vi); %vetor momento angular específico 
h = norm(hv); %modulo momento angular 
i = acos(dot([0,0,1],hv)/h); %(inclinação radianos) 
igraus = (i*180/pi); 
nv = cross([0,0,1],hv);%vetor n passando pela linha dos nodos 
n = norm(nv); %modulo vetor nv 
w = acos((dot(nv,ev))/e*n);% argumento do perigeu 
 
 
anom = acos(dot(ev,si)/e*r); %anomalia verdadeira 
no = acos(dot([1,0,0],nv)/n); 
nograu = no*180/pi; 
h = a - 6378; %altitude 
 
%mostando o tempo que um satélite leva para completar uma volta(período): 
disp(periodo) 
cputime