1. Составить математическую модель невозмущённого движения низкоорбитального космического аппарата (КА). 2. Реализовать модель на языке C++. 3. Провести моделирование 100 витков траектории и определить накопленную погрешность параметров полёта. Исходные

Формулировка задания
1. Составить математическую модель невозмущённого движения низкоорбитального космического аппарата.
2. Реализовать модель на языке C++.
3. Провести моделирование 100 витков траектории и определить накопленную погрешность параметров полёта.
Исходные данные
В качестве исходной орбиты принять круговую околоземную орбиту высотой ℎ и наклонением 𝑖 по варианту. Гравитационный параметр Земли 𝜇з принять равным 𝜇з=398600,45км3𝑐2, фигуру Земли считать сферической радиусом 6 371 км.
Обязательные требования к модели
1. Угловая ориентация КА не учитывается. Моделируется движение только центра масс. Модель движения должна быть пространственной.
2. Дифференциальные уравнения ММД интегрировать с шагом по времени не более Δ𝑡𝑖𝑛𝑡=0,01 при использовании алгоритма с постоянным шагом или с абсолютной точностью 0,1 и относительной точностью 1е-06 при использовании алгоритма с адаптивным шагом.
3. Вращение Земли не учитывается. Исходные данные: Высота, км - 226, угол наклона орбиты, град - 16. Пример начала кода на пайтон:
import numpy as np #pip3 install numpy в терминале
from matplotlib import pyplot as plt
import math

mu_m3s2 = 398600.4415e+09

# 1. задать уравнение ММД (правые части СДУ)

def rightsides(t:np.float64, statevec : np.ndarray) -> np.ndarray:

""" Statevec = [x y z Vx Vy Vz] """

global mu_m3s2
r = np.linalg.norm(statevec[0:3])
derivs = np.zeros(shape=[6])
derivs[0:3] = statevec[3:6]
derivs[3:6] = -mu_m3s2 /r**3 * statevec[0:3]
return dervis

r0 = 6371e+03+200e+03
V0_circ = math.sqrt(mu_m3s2/r0)
print(f'V0_circ = {V0_circ}')
Vy0 = V0_circ * math.cos(math.radians(45))
Vz0 = V0_circ * math.sin(math.radians(45))

t0 = 0.0
statevec0 = np.array([
r0, 0.0, 0.0, #x y z
0.0, Vy0, Vz0 # Vx Vy Vz
])


print("rightsides= ", rightsides(t0, statevec0))


# 2. Задать начальные условия

r0 = 6371e+03+200e+03
V0_circ = math.sqrt(mu_m3s2/r0)
print(f'V0_circ = {V0_circ}')
Vy0 = V0_circ * math.cos(math.radians(45))
Vz0 = V0_circ * math.sin(math.radians(45))

t0 = 0.0
statevec0 = np.array([
r0, 0.0, 0.0, #x y z
0.0, Vy0, Vz0 # Vx Vy Vz
])

# 3. Запустить цикл решения СДУ

revs = 0 # виток
while revs < 100:
#TODO расчёт одного шага интегрирования
if(z_old<0 and z>0):
revs += 1
pass

# 4. Визуализация, анализ