>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'])