# Newton's and Kepler's laws

# Equlibrium fo Two Forces
restart;
R1 := G*m[1]*m[x]/rx^2 = G*m[x]*m[2]/(r - rx)^2;
Sr := solve(R1, rx);
rx := factor(Sr[2]):
rnx := simplify(subs(r = 1, m[2] = N^2*m[1], rx),symbolic):
rnx := subs(N = sqrt(n), rnx);
plot(rnx, n=0..10, labels=[`n`, 'rx']);

# Equlibrium fo Three Forces
r1 := sqrt((x1 - x)^2 + y^2): 
r2 := sqrt((x2 - x)^2 + y^2): 
r3 := sqrt(x^2 + (y3 - y)^2):
F1 := G*m1*m/r1^2: 
F2 := G*m2*m/r2^2: 
F3 := G*m3*m/r3^2:
F1x := F1*(x1 - x)/r1: F1y := -F1*y/r1:
F2x := F2*(x2 - x)/r2: F2y:= -F2*y/r2:
F3x := -F3*x/r3: F3y := F3*(y3 - y)/r3:
Rx := F1x + F2x + F3x=0;  Ry := F1y + F2y + F3y=0;
System := subs(G = 1, m = 1, m1 = 12, x1 = -5, m2 = 9, 
               x2 = 3, m3 = 6, y3 = 8, {Rx, Ry});
with(plots):
implicitplot(System, x = -4..2, y = 0..6, grid=[60,60], 
             labels = [`x`,`y`], color = black);
Digits:=20:
NS1 := fsolve(System, {x, y}, {x = -1..0, y = 0..1});
NS2 := fsolve(System, {x, y}, {x = -1..0, y = 4..5});
evalf(subs(NS1, System)); evalf(subs(NS2, System));

# Equlibrium fo Three Forces Computed from the Potential Energy
Digits:=10:
U := - G*m1/r1 - G*m2/r2 - G*m3/r3:
with(linalg):
g := grad(U,[x,y]);
Un := subs(G = 1, m = 1, m1 = 12, x1 = -5, 
           m2 = 9, x2 = 3, m3 = 6, y3 = 8, U);
gn := grad(Un, [x, y]);
implicitplot({seq(Un = -i/6, i = 24..50)}, x=-6..3, y=-1..9, 
  scaling = constrained, labels = [`x`, `y`], color = black, 
  grid=[60,60]);
plot3d(Un, x = -6..3, y = -1..9, view = -12..-2, 
  orientation = [-100, 75], style = hidden, color = black, 
  grid=[60,60], axes = boxed, labels = [`x`, `y`, 'U']);
NS1 := fsolve({gn[1], gn[2]}, {x, y}, {x = -1..0, y= 0..1});
NS2 := fsolve({gn[1], gn[2]}, {x, y}, {x = -1..0, y= 4..5});
subs(NS1, hessian(Un, [x, y]));
subs(NS2, hessian(Un, [x, y]));

# Gravitation of the Massive Line Segment
restart;
with(plots): with(linalg):
r := sqrt((xi - x)^2 +y^2 +z^2):
sigma := m/L: # linear mass density
U := Int(G*sigma/r, xi = -L/2..L/2);
U := factor(value(U));
Us := subs(G=1,m=1,L=1, z=1/5, U);
plot3d(Us, x = -2..2, y=-2..2, axes = boxed, style = hidden, 
        color = black, orientation = [-80, -130], 
        grid=[35,35], labels = [`x`, `y`, 'U']);
p1 := implicitplot({seq(Us = i/10, i = 5..30)}, x = -1..1, 
        y = -1..1, color = black, scaling = constrained,grid=[40,40]):
p2 := gradplot(Us, x = -1..1, y = -1..1, color = black, 
        scaling = constrained, arrows = SLIM, grid=[15,15]):
display({p1,p2},labels = [`x`,`y`]); 
Uinf := Limit(U, L=0);
Uinf := value(Uinf);
ginf := grad(Uinf, [x, y, z]);
abs(ginf) = simplify(norm(ginf, 2), symbolic);
Us := subs(G = 1, m = 1, L = 1, U):
D2r := [diff(x(t), t, t), diff(y(t), t, t), diff(z(t), t, t)]:
g := subs(x = x(t), y = y(t), z = z(t), grad(Us, [x, y, z])):
IniC:= x(0) = 1, D(x)(0) = 0, y(0) = 0, D(y)(0) = 1, 
       z(0) = 3/4, D(z)(0) = 0:
Ns := dsolve({seq(D2r[i] = g[i], i = 1..3), IniC},
             {x(t), y(t), z(t)}, numeric);
for i from 0 to 500 do;
  T := i/12.5;
  NsT := Ns(T):
  X[i] := subs(NsT,x(t)); Vx[i] := subs(NsT,diff(x(t),t));
  Y[i] := subs(NsT,y(t)); Vy[i] := subs(NsT,diff(y(t),t));
  Z[i] := subs(NsT,z(t)); Vz[i] := subs(NsT,diff(z(t),t));
  KepVec[i] := convert(crossprod([X[i], Y[i] ,Z[i]], 
    [Vx[i], Vy[i], Vz[i]]), list)[2..3];
  KepAbs[i] := norm(KepVec[i], 2);  
od:
spacecurve({[seq([X[i], Y[i], Z[i]], i = 0..500)], 
  [[-1/2, 0, 0], [1/2, 0, 0]]},  labels=[`x`, `y`, `z`],color=black,axes=boxed);
plot([seq(KepVec[i], i = 0..500)], labels=[`y`, `z`],color=black);
plot([seq([i/25, KepAbs[i]], i = 0..500)],
  labels = ['t', 'MofI']);
IniC := x(0) = 1, D(x)(0) = 0, y(0) = 0, D(y)(0) = 1, 
        z(0) = 1/2, D(z)(0) = 1/2:

# The Earth Satellite
restart;
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):
G:=6.67*10^(-11): Mz:=6*10^24:
Inic := x(0) = 7*10^(6), D(x)(0)=0, y(0)=0, D(y)(0)=9*10^3:
Digits := 15:
Ns := dsolve({dx, dy, Inic}, {x(t),y(t)}, numeric):
for i from 0 to 400 do; 
  T := i*40;
  NsT := Ns(T): 
  X[i] := subs(NsT,x(t)); Vx[i] := subs(NsT,diff(x(t),t));
  Y[i] := subs(NsT,y(t)); Vy[i] := subs(NsT,diff(y(t),t));
  MofI[i] := X[i]*Vy[i]-Y[i]*Vx[i]; 
od: 
with(plots):   
p1 := polarplot(6378*10^3, phi = 0..2*Pi):
p2 := plot([seq([X[i], Y[i]], i = 0..327)], thickness=2):
display({p1, p2}, labels = ['x', 'y'], scaling = constrained);
plot([seq([i*40,MofI[i] - 0.63*10^11], i = 0..400)], 
  labels=['t','Wp']);
tx := evalf((t1+t2)/2): 
t1 := 326*40: t2 := 327*40:
y1 := subs(Ns(t1),y(t));  
y2 := subs(Ns(t2),y(t));
while (t1 < tx) and (tx < t2) do:
  yx := subs(Ns(tx),y(t)):
  if yx > 0 then
    y2 := yx: t2 := tx:
   else 
    y1 := yx: t1 := tx: 
  fi;
od:
Tx = evalf(tx);
Hx := floor(tx/3600):
Mx := floor((tx - Hx*3600)/60):
Sx := tx - Hx*3600 - Mx*60:
Hx, Mx, Sx;
XS := [seq(X[i], i = 0..328)]:   YS := [seq(Y[i], i = 0..328)]:
VxS := [seq(Vx[i], i = 0..328)]: VyS := [seq(Vy[i], i = 0..328)]: 
save(G, Mz, XS, YS, VxS, VyS,`orbit.sav`);

# The Earth Satellite, Second Solution
restart;with(plots):
read(`orbit.sav`):
ax := -G*Mz*x/(x^2 + y^2)^(3/2):
ay := -G*Mz*y/(x^2 + y^2)^(3/2):
i := 'i': j := i + 1: 
for k from 0 to 3 do:
  x := 7*10^6: Vx := 0: 
  y := 0: Vy := 9000: 
  dt := evalf(1/2^k); 
  for i from 0 to 328 do:
     X[i] := evalf(x); Y[i] := evalf(y): 
     for n from 1 to 40*2^k do:
        x := evalf(ax*dt^2/2 + Vx*dt + x);
        y := evalf(ay*dt^2/2 + Vy*dt + y);
        Vx := evalf(ax*dt + Vx); 
        Vy := evalf(ay*dt + Vy);
     od:
     if i mod 41 = 0 then  
        dX[k,i]:=X[i]-XS[j]; dY[k,i]:=Y[i]-YS[j];
     fi;
  od:
  p[k] := plot([seq([(X[i] - XS[j])/1000, 
        (Y[i] - YS[j])/1000], i = 0..328)], color = black):
od:
p1 := display({seq(p[k], k = 0..3)}, thickness = 3):p1;
SI := [seq(i*41, i = 0..8)]:
p2 := plot({seq([seq([dX[k, i]/1000, dY[k, i]/1000], 
       k = 0..1), [0, 0]], i = SI)}, color = black): 
display({p1, p2}, scaling = constrained, labels = ['dx', 'dy']);
display({p1, p2}, view = [-0.1..0.5, -0.4..0.2], 
       scaling = constrained, labels = ['dx', 'dy']);

# The Lost Screw
restart;
read(`orbit.sav`):
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): 
for i from 0 to 328 do; 
   T := i*40;  NsT := Ns(T); 
   X[i] := subs(NsT,x(t)); Vx[i] := subs(NsT,diff(x(t),t)); 
   Y[i] := subs(NsT,y(t)); Vy[i] := subs(NsT,diff(y(t),t));
od:
i := 'i': j := i + 1:
plot([seq([X[i] - XS[j], Y[i] - YS[j]], i = 0..328)], 
     labels = [`d x`, `d y`]);
plot([seq([Vx[i] - VxS[j], Vy[i] - VyS[j]], i = 0..328)], 
     labels = [`d Vx`, `d Vy`]); 
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[j])*Sin:
GCY := (X[i] - XS[j])*Sin + (Y[i] - YS[j])*Cos:
GCVx := (Vx[i] - VxS[j])*Cos - (Vy[i] - VyS[j])*Sin:
GCVy := (Vx[i] - VxS[j])*Sin + (Vy[i] - VyS[j])*Cos:
plot([seq([GCX, GCY], i = 0..328)], labels = [`d x`, `d y`]);
plot([seq([GCVx, GCVy], i = 0..327)], labels = [`d Vx`,`d Vy`]);