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

CINEMÁTICA

COORDENADAS GENERALIZADAS E QUASI-VELOCIDADES

𝕢[t_] = {x[t], y[t], ψ[t], ϕ[t], θ[t]};
𝕧[t_] = {(Subscript[ω, 1])[t], (Subscript[ω, 2])[t], (Subscript[ω, 3])[t]};

MATRIZES DE ROTAÇÃO (N ↔ E)

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

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"]) = 𝕧[t]

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ÇÕES CINEMÁTICAS

ce = Solve[{(Subsuperscript[v, "G", "N"]) == (Subsuperscript[ℝ, "E", "N"]).(Subsuperscript[v, "G", "E"]), (Superscript[ω, "E"]) == (Subsuperscript[ω, "a", "E"]) + (Subsuperscript[ω, "r" , "E"])}, 𝕢'[t]] // First // FullSimplify

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"]));
D[(Subsuperscript[H, "C", "E"]), t] + (Subsuperscript[ω, "a", "E"])\[Cross](Subsuperscript[H, "C", "E"]) /. ce // FullSimplify
(Subscript[M, "C"]) = (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"])) - (Subscript[M, "C"])) /.ce // FullSimplify;
TQMA /. notation // TableForm

EQUAÇÕES DIFERENCIAIS DE MOVIMENTO

{m𝔽, 𝕄} = Normal @ CoefficientArrays[TQMA, 𝕧'[t]];
𝕄 // MatrixForm

𝔽 = -m𝔽 // FullSimplify;
𝔽 // MatrixForm

EOMa[t_] = Flatten @ Solve[
	(# == 0)& /@ TQMA,
	𝕧'[t]
	] // Collect[#, 𝕧[t], FullSimplify]&;
EOMa[t] // TableForm

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

FORMA DE ESPAÇO DE ESTADOS

𝕩[t_] = {ϕ[t], (Subscript[ω, 1])[t], (Subscript[ω, 2])[t], (Subscript[ω, 3])[t]};
𝕗 = (𝕩'[t]/.ce/.EOMa[t]);
𝕗 // TableForm

Eigenvalues[(({{0, 1, 0, 0}, {(Subscript[A, 21]), 0, (Subscript[A, 23]), (Subscript[A, 24])}, {0, (Subscript[A, 34]), 0, 0}, {0, (Subscript[A, 42]), 0, 0}}))]

CharacteristicPolynomial[(({{0, 1, 0, 0}, {(Subscript[A, 21]), 0, (Subscript[A, 23]), (Subscript[A, 24])}, {0, (Subscript[A, 34]), 0, 0}, {0, (Subscript[A, 42]), 0, 0}})), s]//FullSimplify

SOLUÇÕES EM REGIME PERMANENTE

PAR = {i -> 1/4, j-> 1/2, (Subscript[Ω, "p" ])-> (Sqrt[33])};
RPcond = Join[{(Subscript[ω, 1])'[t] -> 0, (Subscript[ω, 2])'[t] -> 0, (Subscript[ω, 3])'[t] -> 0}, {(Subscript[ω, 1])[t] -> 0, (Subscript[ω, 2])[t] -> (Subscript[Ω, 2]), (Subscript[ω, 3])[t] -> (Subscript[Ω, 3]), ϕ[t] -> Φ, g -> r ((Subscript[Ω, "p"])^2)}]
TQMA/. RPcond

TIPO I:  (Subscript[Ω, 3]) ≠0

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

LINEARIZAÇÃO

TQMAL1 = Linearize[TQMA, RPcond] /. RP1 //Simplify;

EOML1[t_] = Flatten @ Solve[
	(# == 0)& /@ TQMAL1,
	𝕧'[t]
	] // Collect[#, 𝕧[t], FullSimplify]&;
EOML1[t] // TableForm

MATRIZ DE ESTADOS 𝔸

𝔸1 = CoefficientArrays[𝕩'[t]/. ce /.EOML1[t], 𝕩[t]][[2]] // Normal;
𝔸1 //FullSimplify // MatrixForm

ANÁLISE DE ESTABILIDADE DE SOLUÇÕES

𝔸1 /. PAR //MatrixForm

Eigenvalues[𝔸1 /. PAR]//TableForm

EIG1[deg_]:= {{Show[
	Plot[
		Sort@(Re/@Eigenvalues[𝔸1 /. PAR /.{Φ-> deg π/180, (Subscript[Ω, 3])-> yr}]), {yr,0,10}, PlotRange->{-20,+20}, 
		PlotStyle->Cool[[1]], GridLines->Automatic, Frame->True, FrameLabel->{"(Subscript[Ω, 3]) (rad/s)",None},
		PlotLabel->"Eigenvalues of disc model (Φ = "<>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[𝔸1 /. PAR /.{Φ-> deg π/180, (Subscript[Ω, 3])-> yr}]),{yr,0,10}, 
		PlotStyle->{Cool[[3]],Dashed},
		PlotLegends-> Placed[LineLegend[{"Im (rad/s)"},LegendMarkerSize->30,LabelStyle->{GrayLevel[.4],FontFamily->"Roboto",FontSize->16}],Before]
		]
	]},
	{Plot[
		(Subscript[Ω, 2 ])/.RP1 //.PAR /.{Φ-> deg π/180, (Subscript[Ω, 3])-> yr}, {yr,0.1,10}, PlotRange->All, 
		PlotStyle->Cool[[5]], GridLines->Automatic, Frame->True, FrameLabel->{"(Subscript[Ω, 3]) (rad/s)",None},
		PlotLabel->"Steady state value of (Subscript[Ω, 2]) (Φ = "<>ToString[deg]<>"º)",
		PlotLegends->Placed[LineLegend[{"(Subscript[Ω, 2]) (rad/s)"},LegendMarkerSize->30,LabelStyle->{GrayLevel[.4],FontFamily->"Roboto",FontSize->16}], Before],
		ImageSize->{450,450},
		AspectRatio->1
		]}}//TableForm

Manipulate[EIG1[ϕ],{ϕ,Range[0,50,2]},ControlType-> Setter]

TIPO II: (Subscript[Ω, 3]) =0, Φ = 0

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

RP2 = {Φ->0, (Subscript[Ω, 3])->0};

LINEARIZAÇÃO

TQMAL2 = Linearize[TQMA, RPcond] /. RP2  //Simplify

EOML2[t_] = Flatten @ Solve[
	(# == 0)& /@ TQMAL2,
	𝕧'[t]
	] // Collect[#, 𝕧[t], FullSimplify]&;
EOML2[t] // TableForm

MATRIZ DE ESTADOS 𝔸

𝔸2 = CoefficientArrays[𝕩'[t]/. ce /.EOML2[t], 𝕩[t]][[2]] // Normal;
𝔸2 // MatrixForm

ANÁLISE DE ESTABILIDADE DE SOLUÇÕES

𝔸2 /. PAR //MatrixForm

Eigenvalues[𝔸2 /. PAR]//TableForm

EIG2 = Show[
	Plot[
		Sort@(Re/@Eigenvalues[𝔸2 /. PAR /.{(Subscript[Ω, 2])-> yr}]), {yr,0,50}, PlotRange->{-80,+80}, 
		PlotStyle->Cool[[1]], GridLines->Automatic, Frame->True, FrameLabel->{"(Subscript[Ω, 2]) (rad/s)",None},
		PlotLabel->"Eigenvalues of disc model (Φ = 0º, (Subscript[Ω, 3]) = 0)",
		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[𝔸2 /. PAR /.{(Subscript[Ω, 2])-> yr}]),{yr,0,50}, 
		PlotStyle->{Cool[[3]],Dashed},
		PlotLegends-> Placed[LineLegend[{"Im (rad/s)"},LegendMarkerSize->30,LabelStyle->{GrayLevel[.4],FontFamily->"Roboto",FontSize->16}],Before]
		]
	]

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]/.ce/.EOMa[t](*Join[𝕣𝕢[t], EOM2[t]]*)/.PAR)], -1]));
WriteString[fileD, "\n"];
Close[fileD];