>restart;`read`(orbit.sav)
> dsnumsort := proc (numpr::list, Coor::list)
 local i, j, n; global C1, C2, C3, V1, V2, V3;
 n := nops(Coor); print("Order of the variables:");
 for i from 2 to 2*n+1 do 
  for j to n do
  if [numpr[i]] = select(has, numpr, diff(Coor[j](t), t)) then C[j] := i-1; V[j] := i;
 print(Coor[j], C[j], " ", diff(Coor[j](t), t), V[j]) 
fi;
od;
od;
end:
 >dx := diff(x(t), t, t) = -G*Mz*x(t)/(x(t)^2+y(t)^2)^(3/2);
 >dy := diff(y(t), t, t) = -G*Mz*y(t)/(x(t)^2+y(t)^2)^(3/2); 
 >IniC := x(0) = 7*10^6, (D(x))(0) = 1, y(0) = 0, (D(y))(0) = 9*10^3;
 >Ns := dsolve({dx, dy, IniC}, {x(t), y(t)}, numeric);
 >dsnumsort(Ns(0), [x, y]):
 >for i from 0 to 328 do T := 40*i; NsT := Ns(T); 
 >X[i], Y[i], Vx[i], Vy[i] := (eval([x(t), y(t), diff(x(t), t), diff(y(t), t)], NsT))[]
 >od:
 >i := 'i'; j := i+1;
 >with(plots);
 > plot([seq([X[i]-XS[j], Y[i]-YS[j]], i = 0 .. 328)], labels = ['dx', 'dy']);
 >plot([seq([Vx[i]-VxS[j], Vy[i]-VyS[j]], i = 0 .. 328)], labels = ['dx', 'dy'])
 > R := sqrt(XS[j]^2+YS[j]^2);
 >sin := YS[j]/R; cos := X[j]/R;
 > GCX := (X[i]-XS[j])*cos-(Y[i]-YS[i])*sin;
 > GCY := (X[i]-XS[j])*sin+(Y[i]-YS[i])*cos;
 > GCVx := (Vx[i]-VxS[j])*cos-(Vy[i]-VyS[i])*sin;
 > GCVy := (Vx[i]-VxS[j])*sin+(Vy[i]-VyS[i])*cos;
 > plot([seq([GCX, GCY], i = 0 .. 328)], labels = ['dx', 'dy']);
 >plot([seq([GCVx, GCVy], i = 0 .. 327)], labels = ['dVx', 'dVy'])