Я пытался написать очень простую программу для расчета времени жизни наноспутника на орбите. Я получил плотность атмосферы как функцию высоты с этого сайта.

ALT(km)         DENSITY(kg/m^3)
0               1.17E+00
20              9.49E-02
40              4.07E-03
60              3.31E-04
80              1.68E-05
100             5.08E-07
120             1.80E-08
140             3.26E-09
160             1.18E-09
180             5.51E-10
200             2.91E-10
220             1.66E-10
240             9.91E-11
260             6.16E-11
280             3.94E-11
300             2.58E-11
320             1.72E-11
340             1.16E-11
360             7.99E-12
380             5.55E-12
400             3.89E-12
420             2.75E-12
440             1.96E-12
460             1.40E-12
480             1.01E-12
500             7.30E-13
520             5.31E-13
540             3.88E-13
560             2.85E-13
580             2.11E-13
600             1.56E-13
620             1.17E-13
640             8.79E-14
660             6.65E-14
680             5.08E-14
700             3.91E-14

• программа использует простой интегратор rk4 с землей как сферой
• программа находит перетаскивание =$1/2 C_d\rho A v^2$, с$C_d$как 2.2
• программа предполагает среднюю солнечную активность
• вектор состояния инициализируется для круговой орбиты
• плотность атмосферы берется из приведенной выше ссылки и линейно интерполируется для промежуточных высот.

$C_d$согласно этой статье, которую я только что просмотрел, говорится, что можно взять между 2.0-2.2. Часть программы, которая вычисляет ускорение, выглядит следующим образом:

#define REARTH (6400.0)       // Radius of earth in Km
#define AREA   (0.1 * 1.0E-6) // 0.1 m^2
#define MASS   (14)           // 14 Kg
#define CD     (2.2)
// mu = 3.986004418E5 km^3/s^2
  pos = subm(sv, range(0, 2), range(0, 0));
  vel = subm(sv, range(3, 5), range(0, 0));

  r = length(pos);
  g = - mu / (r*r*r) * pos;
  alt = r - REARTH;
  // gd() returns density in kg/km^3 
  pho = gd(alt);

  drag = -  1.0 * CD * pho * AREA / MASS * vel * length(vel) / 2.0;

  dydx(0) = vel(0);
  dydx(1) = vel(1);
  dydx(2) = vel(2);

  gpdrag = g + drag;

  dydx(3) = gpdrag(0);
  dydx(4) = gpdrag(1);
  dydx(5) = gpdrag(2);

когда я запускаю эту программу, в течение 70 дней орбита уменьшается до менее 100,0 км на 350,0 км. Кажется, это слишком мало времени. Я предполагаю, что здесь что-то не так или совсем не так?

Я предположил, что земля сферическая, разница между полюсами и экватором всего 22 км. Итак, имеет ли значение для этого вывода значение SSPO или экваториальной орбиты?

Я добавил графики для различных высот круговой орбиты.