example.m

function example
  R = 6371;

  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 = [90 60 180];

  color_track = ['r', 'b' 'g'];
  color_bind = {'--b', '--g', '--r'};

  len = max(T);

  count = length(a);
  sp_x = sp_y = sp_z = zeros(count, len);
  max_value = 0;
  for id = 1:length(a)
    [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

  for id = 1:length(a)
    sp_x(id, :) /= max_value;
    sp_y(id, :) /= max_value;
    sp_z(id, :) /= max_value;
  end

  R /= max_value;

  [X, Y, Z] = sphere(25);
  mesh(R*X, R*Y, R*Z);

  axis ([-1 1 -1 1 -1 1]);
  hold on;
  xlabel('x');
  ylabel('y');
  zlabel('z');


  count = length(a);
  k = 1;
  r = R / 10;
  x0 = y0 = z0 = h = sputnik = zeros(1, count);
  for id = 1:length(a)
    sputnik(id) = surf(X*r+sp_x(id, k), Y*r+sp_y(id, k), Z*r+sp_z(id, k));
    h(id) = plot3([0 sp_x(id, k)], [0 sp_y(id, k)], [0 sp_z(id, k)], color_bind{id});
    x0(id) = sp_x(id, k);
    y0(id) = sp_y(id, k);
    z0(id) = sp_z(id, k);
  end

  k = ones(1, count);
  for t = 2:len
    for id = 1:count
      k(id)++;
      if (k(id) > T(id))
        k(id) = 1;
      end
      plot3([x0(id), sp_x(id, k(id))], [y0(id), sp_y(id, k(id))], [z0(id), sp_z(id, k(id))], color_track(id));
%      disp(size(x0));
%      disp(size(sp_x(id, k)));
      x0(id) = sp_x(id, k(id));
      y0(id) = sp_y(id, k(id));
      z0(id) = sp_z(id, k(id));
      delete(sputnik(id));
      delete(h(id));
      sputnik(id) = surf(X*r+sp_x(id, k(id)), Y*r+sp_y(id, k(id)), Z*r+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)
  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