SetDirectory @ NotebookDirectory[];
<<Rotation.m
<<Linearize.m
<<PlotCustom.m
<<NotationRules.m

CINEMÁTICA

MATRIZES DE ROTAÇÃO (N ↔ E)

(Subsuperscript[ℝ, "E", "N"]) = Rotation["z", "x"][ψ[t], ϕ[t]];
(Subsuperscript[ℝ, "N", "E"]) = Transpose[(Subsuperscript[ℝ, "E", "N"])];

(Subsuperscript[ℝ, "E", "N"])//TraditionalForm

POSIÇÃO DO CENTRO DE MASSA

(Subsuperscript[p, "C/O", "N"]) = {x[t], y[t], 0}; 
(Subsuperscript[p, "G/C", "E"]) = {0, 0, -r};
(Subsuperscript[p, "G/O", "N"]) = (Subsuperscript[ℝ, "E", "N"]).(Subsuperscript[p, "G/C", "E"]) + (Subsuperscript[p, "C/O", "N"])

VELOCIDADES ANGULARES

(Subsuperscript[ω, "a", "E"]) = (Subsuperscript[ℝ, "N", "E"]).{0, 0, ψ'[t]} + {ϕ'[t], 0, 0}
(Subsuperscript[ω, "r" , "E"])= {0, θ'[t], 0}
(Superscript[ω, "E"]) =(Subsuperscript[ ω, "a", "E"]) + (Subsuperscript[ω, "r", "E"])

VELOCIDADE DO CENTRO DE MASSA

(Subsuperscript[v, "G", "N"]) = D[(Subsuperscript[p, "G/O", "N"]),t] // FullSimplify
(Subsuperscript[v, "G", "E"]) = (Superscript[ω, "E"])\[Cross]((Subsuperscript[p, "G/C", "E"])) // FullSimplify

EQUAÇÃO DE VÍNCULO

ce = Solve[(Subsuperscript[v, "G", "N"]) == (Subsuperscript[ℝ, "E", "N"]).(Subsuperscript[v, "G", "E"]), {x'[t],y'[t]}] //First

VELOCIDADE DO PONTO C

(Subsuperscript[v, "C", "N"]) = {x'[t], y'[t], 0} /. ce
(Subsuperscript[v, "C", "E"]) = (Subsuperscript[ℝ, "N", "E"]) . (Subsuperscript[v, "C", "N"]) // FullSimplify

DINÂMICA

TQMA – POLO C

(Subsuperscript[𝕀, "C", "E"]) = DiagonalMatrix[{i+1, j+1, i}] m (r^2);
(Subsuperscript[H, "C", "E"]) = ((Subsuperscript[𝕀, "C", "E"]).(Superscript[ω, "E"]))
(Subsuperscript[M, "C", "E"]) = (Subsuperscript[p, "G/C", "E"]) \[Cross]((Subsuperscript[ℝ, "N", "E"]).{0, 0, m g})

TQMA = (D[(Subsuperscript[H, "C", "E"]), t] + (Subsuperscript[ω, "a", "E"])\[Cross](Subsuperscript[H, "C", "E"]) + (Subsuperscript[v, "C", "E"])\[Cross](m (Subsuperscript[v, "G", "E"])) - (Subsuperscript[M, "C", "E"])) // FullSimplify;
TQMA /. notation // TableForm

EQUAÇÕES DIFERENCIAIS DE MOVIMENTO

EOMa[t_] = Flatten @ Solve[
	(# == 0)& /@ TQMA,
	{ψ''[t], ϕ''[t], θ''[t]}
	] // Collect[#, {ψ'[t], ϕ'[t], θ'[t]}, FullSimplify]&;
EOMa[t] // TableForm

EOMa[t] /. notation // FullSimplify // TableForm

SOLUÇÕES EM REGIME PERMANENTE

RPcond = {ψ''[t] -> 0, ϕ''[t] -> 0, θ''[t] -> 0, ψ'[t] -> (Subscript[Ω, "y"]), ϕ'[t] -> 0, θ'[t] -> (Subscript[Ω, "s"]), ϕ[t] -> (SuperStar[ϕ]), g -> r ((Subscript[Ω, "p"])^2)};
TQMA/. RPcond

RPsol = Solve[(# == 0)& /@ (TQMA/. RPcond), {(Subscript[Ω, "s"])}] //First //FullSimplify

PAR = {i -> 1/4, j-> 1/2, (Subscript[Ω, "p" ])-> 1};

ContourPlot[(Subscript[Ω, "s"])/.RPsol/.PAR,{(Subscript[Ω, "y"]),0.1, 10},{(SuperStar[ϕ]),0,0.95 π/2},
PlotTheme->"Detailed",ColorFunction->"DarkRainbow",GridLines->{Range[0,10,0.5],Range[0,1.50,0.05]}, 
PlotLegends->BarLegend[Automatic,LabelStyle->{FontSize->20}],FrameLabel->{"(Subscript[Ω, y]) [rad/s]","(SuperStar[ϕ]) [rad]"}, PlotLabel->"(Subscript[Ω, s])"]

TQMA/. RPcond /.{(Subscript[Ω, "y"])->0}

TQMA/. RPcond/. RPsol //Simplify

(Subsuperscript[v, "G", "E"]) /. RPcond // FullSimplify

LINEARIZAÇÃO

TQMAL = Linearize[TQMA, RPcond] /.RPsol //Simplify;

EOMl[t_] = Flatten @ Solve[
	(# == 0)& /@ TQMAL,
	{ψ''[t], ϕ''[t], θ''[t]}
	] // Collect[#, {ψ'[t], ϕ'[t], θ'[t]}, FullSimplify]&;
EOMl[t] // TableForm

MODELO MATEMÁTICO - FORMA DE ESPAÇO DE ESTADOS

MATRIZ DE ESTADOS 𝔸

𝕩[t_] = {ϕ[t], ψ'[t], ϕ'[t], θ'[t]};
𝔸 = CoefficientArrays[𝕩'[t]/.EOMl[t], 𝕩[t]][[2]] // Normal;
𝔸 // MatrixForm

ANÁLISE DE ESTABILIDADE

𝔸 /. PAR //MatrixForm

Eigenvalues[𝔸 /. PAR]//TableForm

EIG[deg_]:= Show[
	Plot[
		Sort@(Re/@Eigenvalues[𝔸 /. PAR /.{(SuperStar[ϕ])-> deg π/180, (Subscript[Ω, "y"])-> yr}]), {yr,0,2}, PlotRange->{-2,+2}, 
		PlotStyle->Cool[[1]], GridLines->Automatic, Frame->True, FrameLabel->{"(Subscript[Ω, y]) (rad/s)",None},
		PlotLabel->"Eigenvalues of disc model ((SuperStar[ϕ]) = "<>ToString[deg]<>"º)",
		PlotLegends->Placed[LineLegend[{"Re (rad/s)"},LegendMarkerSize->30,LabelStyle->{GrayLevel[.4],FontFamily->"Roboto",FontSize->16}], Before],
		ImageSize->{450,450},
		AspectRatio->1
		],
	Plot[
		Sort@(*RedundantElim@*)Chop@(Im/@Eigenvalues[𝔸 /. PAR /.{(SuperStar[ϕ])-> deg π/180, (Subscript[Ω, "y"])-> yr}]),{yr,0,2}, 
		PlotStyle->{Cool[[3]],Dashed},
		PlotLegends-> Placed[LineLegend[{"Im (rad/s)"},LegendMarkerSize->30,LabelStyle->{GrayLevel[.4],FontFamily->"Roboto",FontSize->16}],Before]
		]
	]

Manipulate[EIG[ϕ],{ϕ,Range[0,30,2]},ControlType-> Setter]

EXPORT TO OCTAVE

RULES

Unprotect[Power];
Format[Power[a_,n_Integer],CForm]:=Distribute[ConstantArray[Hold[a],n],Hold,List,HoldForm,Times]
Protect[Power];

crule = MapThread[#1->#2&, {#, x/@((Range[Length@#]))}& @ 𝕩[t]];
csorule = {"Cos"->"cos", "Sin"->"sin", "Tan"->"tan", "Sec"->"sec"};
csqrule = {};

EXPORT FILE

fileD = OpenWrite["export/disc_OCT.m"];
(*WriteString[fileD, (#/.(({ii_}->vv_):>(StringReplace[StringReplace[ToString @ StringForm["dy(``) = ``;", ii, (vv /. crule //CForm)], csqrule], csorule]<>"\n")))]&/@ (Normal@(Select[#,Not @(NumericQ[#])&]&@Association@Delete[ArrayRules@SparseArray[(𝕪'[t]/.EOM[t])], -1]));*)
WriteString[fileD, (#/.(({ii_}->vv_):>(StringReplace[StringReplace[ToString @ StringForm["    ``;", (vv /. crule //CForm)], csqrule], csorule]<>"\n")))]&/@ (Normal@(Select[#,Not @(NumericQ[#])&]&@Association@Delete[ArrayRules@SparseArray[(𝕩'[t]/.EOMa[t](*Join[𝕣𝕢[t], EOM2[t]]*)/.PAR)], -1]));
WriteString[fileD, "\n"];
Close[fileD];