Unprotect["Global`*"]
ClearAll["Global`*"]
Remove["Global`*"]
<<Notation`
Notation[(NotationTemplateTag[(Subscript[x_, i_·])]) ⟺ (NotationTemplateTag[x_[[i_]]])]
Notation[(NotationTemplateTag[(Subscript[x_, i_·j_·])]) ⟺ (NotationTemplateTag[x_[[i_,j_]]])]
Notation[(NotationTemplateTag[(Subscript[x_, i_·j_·k_·])]) ⟺ (NotationTemplateTag[x_[[i_,j_,k_]]])]

Notation[(NotationTemplateTag[(Subscript[x_[t_], i_·])]) ⟺ (NotationTemplateTag[x_[t_][[i_]]])]
Notation[(NotationTemplateTag[(Subscript[x_[t_], i_·j_·])]) ⟺ (NotationTemplateTag[x_[t_][[i_,j_]]])]
Notation[(NotationTemplateTag[(Subscript[x_[t_], i_·j_·k_·])]) ⟺ (NotationTemplateTag[x_[t_][[i_,j_,k_]]])]

Symbolize[(NotationTemplateTag[(Subscript[x_, t])])]
Symbolize[(NotationTemplateTag[(Underscript[x_, (UnderBar[_])])])]
Symbolize[(NotationTemplateTag[(Subscript[(Underscript[x_, (UnderBar[_])]), 0])])]
Symbolize[(NotationTemplateTag[((x_)^𝒪)])]
Symbolize[(NotationTemplateTag[((UnderBar[x_])^𝒪)])]
Symbolize[(NotationTemplateTag[(OverHat[x_])])]
Symbolize[(NotationTemplateTag[(Underoverscript[x_, (_), (_)])])]
Symbolize[(NotationTemplateTag[(Underoverscript[x_, (_), (^)])])]
Symbolize[(NotationTemplateTag[(UnderBar[x_])])]
Symbolize[(NotationTemplateTag[(Subscript[(OverHat[x_]), i_])])]
<<HaneeshPackages`
ϵ=LeviCivitaTensor[3]//Normal;
δ=Table[KroneckerDelta[i,j],{i,1,3},{j,1,3}];

ψ[t_]=Ω t;


(Underscript[R, (UnderBar[_])])::usage="Matrix representation of the rotation tensor in the (Subscript[(OverHat[(UnderBar[E])]), i]) basis";
Echo[((Underscript[R, (UnderBar[_])])=RotationMatrix[ψ[t],{0,0,1}].RotationMatrix[θ[t],{0,1,0}])//MatrixForm,"Rotation tensor \n"];



(Underscript[Ω, (UnderBar[_])])::usage="Matrix representation of the angular velocity tensor in the (Subscript[(OverHat[(UnderBar[E])]), i]) basis";
Echo[((Underscript[Ω, (UnderBar[_])])=(D[(Underscript[R, (UnderBar[_])]),t].(Underscript[R, (UnderBar[_])])\[Transpose])//Simplify )//MatrixForm,   "Angular velocity tensor \n"];

(UnderBar[ω])::usage="Matrix representation of the angular velocity vector with respect to the (Subscript[(OverHat[(UnderBar[E])]), i]) basis";
(UnderBar[ω])=-(1/2)Table[Sum[(Subscript[ϵ, i·j·k·])(Subscript[(Underscript[Ω, (UnderBar[_])]), j·k·]),{j,1,3},{k,1,3}],{i,1,3}]; 

((Subscript[(Underscript[J, (UnderBar[_])]), 0]))::usage="Matrix representation of the inertia tensor at time=0 with respect to the (Subscript[(OverHat[(UnderBar[E])]), i]) basis";
Echo[((Subscript[(Underscript[J, (UnderBar[_])]), 0])={{((m (ℓ^2))/12),0,0},{0,((m (ℓ^2))/12),0},{0,0,0}})//MatrixForm,"Inertia tensor (Subscript[(Underscript[J, (UnderBar[_])]), 0]) \n"];

(Underscript[J, (UnderBar[_])])::usage="Matrix representation of the inertia tensor with respect to the (Subscript[(OverHat[(UnderBar[E])]), i]) basis";
Echo[((Underscript[J, (UnderBar[_])])=((Underscript[R, (UnderBar[_])]).((Subscript[(Underscript[J, (UnderBar[_])]), 0])).(Underscript[R, (UnderBar[_])])\[Transpose])//Simplify)//MatrixForm,"Inertia tensor (Underscript[J, (UnderBar[_])]) \n"];
Echo[((UnderBar[H])=(Underscript[J, (UnderBar[_])]).(UnderBar[ω])//Simplify)//MatrixForm,"Angular momentum vector (UnderBar[H]) \n"];





(*Angular velocity vector*)

Protect[m,ℓ, g];
m::usage="mass of the bar, constant";
l::usage="length of the bar, constant";
g::usage="acceleration due to gravity, constant";
Echo[((Underoverscript[X, (_), (_)])={0,0,ℓ/2}//Simplify)//MatrixForm,"Position vector of the center of mass in reference configuration, \!(Underoverscript[X, (_), (_)]) \n"];
Echo[((Underoverscript[x, (_), (_)])=(Underscript[R, (UnderBar[_])]).(Underoverscript[X, (_), (_)]) //Simplify)//MatrixForm,"Position vector of the center of mass in spatial configuration, (Underoverscript[x, (_), (_)]) \n"];
Echo[((Underoverscript[v, (_), (_)])=D[(Underoverscript[x, (_), (_)]),t]//Simplify)//MatrixForm,  "velocity vector of the center of mass, (Underoverscript[v, (_), (_)]) \n"];
Echo[((UnderBar[G])=(m (Underoverscript[v, (_), (_)]) )//Simplify)//MatrixForm,    "linear momentum, (UnderBar[G]) \n"];
Echo[((Underoverscript[F, (_), (_)])={0,0,m g})//MatrixForm,"Total force, (Underoverscript[F, (_), (_)]) \n"];
Echo[(((UnderBar[H])^𝒪)=((UnderBar[H])+(Underoverscript[x, (_), (_)])\[Cross](UnderBar[G]) )//Simplify)//MatrixForm,"Angular momentum about origin, ((UnderBar[H])^𝒪) \n"];
Echo[(((UnderBar[M])^𝒪)=((Underoverscript[x, (_), (_)])\[Cross](Underoverscript[F, (_), (_)]))//Simplify)//MatrixForm,"Net moment about about origin, ((UnderBar[M])^𝒪) \n"];

Res::usage="matrix reprsentation of the residual moment vector, ((OverDot[((UnderBar[H])^𝒪)]))-((UnderBar[M])^𝒪), on the solid in the (Subscript[(OverHat[(UnderBar[E])]), i]) basis";
Echo[(Res=D[((UnderBar[H])^𝒪),t]-((UnderBar[M])^𝒪))//Simplify//MatrixForm,"Balance of angular momementum about origin \n"];

Echo[((Underscript[R, (UnderBar[_])])\[Transpose]).Res//Simplify//MatrixForm,"Matrix reprsentation of the residual moment vector, ((OverDot[((UnderBar[H])^𝒪)]))-((UnderBar[M])^𝒪), on the solid in the (Subscript[(UnderBar[(OverHat[e])]), i]) basis \n"];