Euler’s Equations

Consider an object with general (Subscript[I, 1]), (Subscript[I, 2 ])and (Subscript[I, 3]) tumbling in outer space with no external torque.
In the lab frame, the whole (OverVector[L]) is conserved, while in the frame of the object, only the
magnitude is constant as the direction of (OverVector[L]) wanders around.
At any time, the three principal directions form the rotation matrix.
Let’s agree to use Euler angles to parameterize these.  Then here is the DiffEQ we’ll have to solve to 
translate knowledge of ω-vector in to a rotation matrix.

$Post=MatrixForm
euler[ϕ_,θ_,ψ_] = RotationMatrix[ϕ,{0,0,1}].RotationMatrix[θ,{1,0,0}].RotationMatrix[ψ,{0,1,0}];
rot[t_] := euler[ϕ[t],θ[t],ψ[t]]
Ωmat= Transpose[ rot[t]].D[rot[t],t]//Simplify (* form the antisymmetric Ω matrix *)
(eQ1:=Ωmat[[1,2]]==-ωvec[t][[3]]; (* map the matrix elements to components of ω vector *)
eQ2 :=Ωmat[[1,3]] == ωvec[t][[2]];
eQ3 := Ωmat[[2,3]] ==-ωvec[t][[1]];)

This time consider a point and stick, as in Morin 9.39

params={m->1,M->2,ell->1.4,h->1.4};
(I1=m (h M/(m+M))^2+M (h m/(m+M))^2/.params;I2=(1/12) M ell^2/.params; I3=I1+I2;
ω1i = 1; ω2i=2; ω3i=0;
Tf=15.3;
eqn1 = I1 ω1'[t] == (I2-I3) ω2[t]*ω3[t];
eqn2 = I2 ω2'[t] == (I3-I1) ω3[t]*ω1[t];
eqn3 = I3 ω3'[t] == (I1-I2) ω1[t]*ω2[t];
soln=NDSolve[{eqn1,eqn2,eqn3,ω1[0]==ω1i, ω2[0]==ω2i, ω3[0]==ω3i},{ω1[t],ω2[t],ω3[t]},{t,0,Tf}];
ωvec[t_] = {ω1[t],ω2[t],ω3[t]}/.soln[[1]];
Plot[{Evaluate[ωvec[t]]},{t,0,Tf},PlotStyle->{Red,Green,Blue}])
(srot=NDSolve[{eQ1,eQ2,eQ3,
ϕ[0]==.01,θ[0]==0.01,ψ[0]==.01},{ϕ[t],θ[t],ψ[t]},{t,0,Tf}][[1]];
transform[t_]=rot[t]/.srot;
myobj[t_]:=If[t>0,{GeometricTransformation[myobj0,transform[t]],{Green,Ball[{ell/2,-h m/(M+m),t},.03],Arrow[{{0,0,0},transform[t].ωvec[t]/5}]}},{myobj0,{Green,Ball[{ell/2,-h m/(M+m),-2t},.03]}}]/.params;
myobj0=With[{y1=h M/(M+m)/.params,y2=h m/(m+M)/.params,x2=ell/2/.params},
(pt1={0,y1,0};pt2a={-x2,-y2,0};pt2b={x2,-y2,0};{{Red,Ball[pt1,.1]},{Blue,Cylinder[{pt2a,pt2b},.07]},{Purple,Ball[pt2a,.07],LightBlue,Ball[pt2b,.07]},{Gray,Ball[{0,0,0},.02]},{Thick,Line[{pt2b,pt1,pt2a}]}})/.params];
Manipulate[Graphics3D[myobj[t],Boxed->False,PlotRange->{{-1,1},{-1,1},{-1,1}}],{t,-3/2,Tf}])

Tides

To make a vector field plot of the familiar 1/r^2 gravity, we can define the -1/r potential ϕ and take (minus) its gradient.
Here we make a basic VectorPlot, and then add a black dot at the center.

ϕ[x_,y_] = -1/Sqrt[x^2+y^2]
g[x_,y_] = - D[ϕ[x,y],{{x,y}}]
VectorPlot[g[x,y],{x,.4,1.4},{y,-.5,.5}]
Graphics[(vfield={{Black,Disk[{0,0},.08]},%[[1]]})]

Now add a small planet of radius ϵ centered at the point (1,0), and zoom in to look
at the local gravitational field, which is pretty close to constant on this scale.

ϵ=.02;
planet={Blue,Opacity[.4],Disk[{1,0},ϵ]};
Graphics[{vfield,planet}]
Graphics[{VectorPlot[g[x,y],{x,1-ϵ,1+ϵ},{y,-ϵ,ϵ}][[1]],planet}]

The tidal force is the difference between this (almost constant) force and the (actually constant) force at the middle:

VectorPlot[g[x,y]-g[1,0],{x,1-ϵ,1+ϵ},{y,-ϵ,ϵ}]
Graphics[{%[[1]],planet}]

Note this force field comes from a potential too:

ϕtide[x_,y_] = ϕ[x,y] + g[1,0].{ x,0}
gtide[x_,y_] = -D[ϕtide[x,y],{{x,y}}]
field=VectorPlot[gtide[x,y],{x,1-ϵ,1+ϵ},{y,-ϵ,ϵ}];
Show[field,Graphics[planet]]

Now add in the local gravity from the we planet itself, and we can plot equipotential shapes.
We choose the gravity to be small, which exaggerates the tidal bulge.

ϕtot[x_,y_] = ϕtide[x,y]- (ϵ^2/4)/Sqrt[(x-1)^2+y^2]; (*play with the size of the ϵ^2 term*)
q=ContourPlot[ϕtot[x,y]==ϕtot[1,ϵ*.95],{x,1-1.2ϵ,1+1.2ϵ},{y,-1.2ϵ,1.2ϵ},ContourStyle->Red];
Show[q,field,Graphics[planet]]

BTW if you put in the dimensionful parameters we have set to 1, you can get a prediction for the difference
between high and low tide.  Adding the contributions from the Moon and Sun, the result is about 1 meter.
In the real world, tide charts are several times larger.  Why?!