example.m

function example
% Программа, демонстрирующая движение нескольких искусственных спутников Земли

  R = 6371; % Радиус земли в километрах

% Перечисляются элементы орбиты каждого из спутников (в примере три спутника):
  count = 3;                     % количество спутников
  a     = [4*R  3*R     2*R];    % Большие полуоси
  e     = [0.3  0.05    0.15];   % Эксцентриситет
  i     = [pi/9 pi/4    pi/3];   % Наклонение
  Omega = [pi/6 pi/3    pi];     % Долгота восходящего узла
  omega = [pi/3 2*pi/3  5*pi/3]; % Аргумент перицентра
  T     = [180  90      60];     % Период обращения

% Цвета для каждого спутника:
  color_track = [   'r',  'b',    'g'];   % Цвет линии орбиты
  color_bind  = { '--r',  '--b',  '--g'}; % Цвет линии, соединяющей спутник и центр земли

% Предварительный расчет координат для каждого из спутников:
  len = max(T); % len - сколько итераций рисовать (выбрана длина максимального периода)
  sp_x = sp_y = sp_z = zeros(count, len); % векторы для декартовых координат спутников
  max_value = 0;  % в процессе расчетов будем искать максимальное значение среди всех координат чтобы пронормировать
  for id = 1:count % в цикле по всем функциям
    % вызов функции, расчитывающей координаты одного спутника:
    [sp_x(id, :), sp_y(id, :), sp_z(id, :)] = one_sputnik(a(id), e(id), i(id), Omega(id), omega(id), T(id), len);
    % поиск и обновление максимума среди координат:
    tm1 = max(abs(sp_x(id, :)));
    tm2 = max(abs(sp_y(id, :)));
    tm3 = max(abs(sp_z(id, :)));
    max_value = max([max_value, tm1, tm2, tm3]);
  end

% Нормировка радиуса Земли:
  R /= max_value;

% Рисование сферы:
  [X, Y, Z] = sphere(25);
  mesh(R*X, R*Y, R*Z);

% Настраивание осей:
  axis ([-1 1 -1 1 -1 1]);
  axis equal;
  hold on;
  xlabel('x');
  ylabel('y');
  zlabel('z');

% Нормировка координат спутников и рисование орбит:
  for id = 1:count
    sp_x(id, :) /= max_value;
    sp_y(id, :) /= max_value;
    sp_z(id, :) /= max_value;
    % рисуем орбиту спутника с номером id:
    plot3(sp_x(id, [1:T(id), 1]), sp_y(id, [1:T(id), 1]), sp_z(id, [1:T(id), 1]), color_track(id));
  end

  r = R / 25; % радиус сферы для отображения спутника
  X *= r;     % сжатие сферы
  Y *= r;
  Z *= r;


% Рисуем каждый спутник и линию, соединяющую спутник и центр Земли в начальный момент времени:
  h = sputnik = zeros(1, count); % вектора, в которые мы будем сохранять ссылки на графические объекты для их удаления
  for id = 1:count
    % рисуем сферу, представляющую спутник:
    sputnik(id) = surf(X+sp_x(id, 1), Y+sp_y(id, 1), Z+sp_z(id, 1));
    % рисуем линию, соединяющую спутник и центр Земли:
    h(id) = plot3([0 sp_x(id, 1)], [0 sp_y(id, 1)], [0 sp_z(id, 1)], color_bind{id});
  end

% счетчик моментов времени для каждого спутника:
  k = ones(1, count);

% основная анимация:
  for t = 2:len % текущий момент времени
    for id = 1:count % цикл по всем спутникам
      k(id)++; % увеличиваем момент времени каждого спутника
      if (k(id) > T(id)) % если превысили период обращения
        k(id) = 1; % то возвращаемся в начало
      end

      % удаление предыдущих и рисование текущих сфер и линий:
      delete(sputnik(id));
      delete(h(id));
      sputnik(id) = surf(X+sp_x(id, k(id)), Y+sp_y(id, k(id)), Z+sp_z(id, k(id)));
      h(id) = plot3([0 sp_x(id, k(id))], [0 sp_y(id, k(id))], [0 sp_z(id, k(id))], color_bind{id});
    end
    pause(0.001); % величина паузы
  end
end

function [x, y, z] = one_sputnik(a, e, i, Omega, omega, T, len)
% функция, расчитывающая координаты одного спутника в каждый из моментов времени
% a - большая полуось орбиты
% e - эксцентриситет
% i - наклонение
% Omega - долгота восходящего узла
% omega - аргумент перицентра
% T - Период обращения
% len - длина массива, которую необходимо выделить
% введена, так как периоды у всех спутников разные, а длина ассивов должна быть одинаковая

% i и Omega не изменяются, поэтому заранее вычисляем их косинусы и синусы:
  cos_Omega = cos(Omega);
  sin_Omega = sin(Omega);

  cos_i = cos(i);
  sin_i = sin(i);

% массив для моментов времени:
  dt = 2*pi/T;
  t = [0:dt:2*pi-dt, zeros(1, len-T)];

% массив координат:
  x = y = z = zeros(1, length(t));

% вычисляем координаты для каждого момента времени:
  for k = 1:length(t)
    r = a*(1-e^2)/(1+e*cos(t(k)));
    sin_omega = sin(t(k)+omega);
    cos_omega = cos(t(k)+omega);
    x(k) = r*(cos_Omega*cos_omega-sin_Omega*sin_omega*cos_i);
    y(k) = r*(sin_Omega*cos_omega+cos_Omega*sin_omega*cos_i);
    z(k) = r*sin_omega*sin_i;
  end
end