example.m

function example
  e = 0.7;
  i = pi/9;
  Omega = pi/6;
  omega = pi/3;

  cos_Omega = cos(Omega);
  sin_Omega = sin(Omega);

  sin_omega = sin(omega);
  cos_omega = cos(omega);

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


  e_x = e*(cos_Omega*cos_omega-sin_Omega*sin_omega*cos_i);
  e_y = e*(sin_Omega*cos_omega+cos_Omega*sin_omega*cos_i);
  e_z = e*(sin_omega*sin_i);

  a = 0;
  b = 2*pi;
  t = [a:(b-a)/180:b];



  x = y = z = zeros(1, length(t));

  R = 6371;
  a = 5*R;

  for k = 1:length(t)
    r = a*(1-e^2)/(1+e*cos(t(k)));
    x(k) = r*(cos_Omega*cos(t(k)+omega)-sin_Omega*sin(t(k)+omega)*cos_i);
    y(k) = r*(sin_Omega*cos(t(k)+omega)+cos_Omega*sin(t(k)+omega)*cos_i);
    z(k) = r*sin(t(k)+omega)*sin_i;
  end

  total_max = max([abs(x), abs(y), abs(z)]);

  x /= total_max;
  y /= total_max;
  z /= total_max;
  R /= total_max;

  [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');

  [X, Y, Z] = sphere(25);
  k = 1;
  r = R / 10;
  sputnik = surf(X*r+x(k), Y*r+y(k), Z*r+z(k));
  h = plot3([0 x(k)], [0 y(k)], [0 z(k)], '-b');
  x0 = x(k);
  y0 = y(k);
  z0 = z(k);
  for k = 2:length(t)
    plot3([x0, x(k)], [y0, y(k)], [z0, z(k)], 'r');
    x0 = x(k);
    y0 = y(k);
    z0 = z(k);
    delete(sputnik);
    delete(h);
    sputnik = surf(X*r+x(k), Y*r+y(k), Z*r+z(k));
    h = plot3([0 x(k)], [0 y(k)], [0 z(k)], '-b');
    pause(0.01);
  end
end