From 29631aa6f0125727e3388e0266ad1a09d3aa6bf9 Mon Sep 17 00:00:00 2001
From: LolloBici <lorenzo.amici@skywarder.eu>
Date: Sat, 15 Feb 2025 23:09:42 +0100
Subject: [PATCH] [stability-bugs] Fixing plots and commenting openRocket
---
common | 2 +-
stabilityAnalysis/mainStabilityAnalysis.m | 2 +-
stabilityAnalysis/src/getXCP.m | 117 ----------------------
stabilityAnalysis/src/stabilityPlots.m | 48 ++++-----
4 files changed, 26 insertions(+), 143 deletions(-)
delete mode 100644 stabilityAnalysis/src/getXCP.m
diff --git a/common b/common
index 01a8a641..7589b269 160000
--- a/common
+++ b/common
@@ -1 +1 @@
-Subproject commit 01a8a6416273f3ad59a1420fc8a1602640c67b4c
+Subproject commit 7589b2690bf66066247dd1236f27eb2ee7c52cbd
diff --git a/stabilityAnalysis/mainStabilityAnalysis.m b/stabilityAnalysis/mainStabilityAnalysis.m
index 8727aa7a..3164c2bb 100644
--- a/stabilityAnalysis/mainStabilityAnalysis.m
+++ b/stabilityAnalysis/mainStabilityAnalysis.m
@@ -47,7 +47,7 @@ if isempty(environment), environment = Environment(mission, rocket.motor); end
if isempty(wind), wind = WindCustom(mission); end
if isempty(settings), settings = Settings('stabilityAnalysis', 'ode'); end
-Settings.read(settings, options, 'stability');
+% Settings.read(settings, options, 'stability');
settings.stability.dissile.alt = environment.z0;
settings.stability.dissile.xcg = rocket.xCg(1);
diff --git a/stabilityAnalysis/src/getXCP.m b/stabilityAnalysis/src/getXCP.m
deleted file mode 100644
index b914aa4c..00000000
--- a/stabilityAnalysis/src/getXCP.m
+++ /dev/null
@@ -1,117 +0,0 @@
-function [XCPvect, qualityIndex, stabilityCoeffs] = getXCP(coeff, iAlpha, iBeta, input, vars)
-% OUTPUT:
-% XCPvect, double [1, 7]:
-% (1) XCPlon
-% (2) XCPlat
-% (3) XCPtot - std
-% (4) XCP tot - MOD
-% (5) XCPlon - aTot
-% (6) XCPlat - aTot
-% (7) XCPor
-% qualityIndex, double[1, 5]:
-% (1) Angle between Maero and Faero [rad]
-% (2) Angle between Fyz and Myz [rad]
-% (3) Angle between Fyz and vRyz [rad]
-% (4) Error index T
-% (5) Magnitude of Fpar
-% stabilityCoeff, double[1, 10]:
-% (1) Cma (2) CNA (3) Cm (4) CN (5) Cnb
-% (6) CYB (7) Cn (8) CY (9) Cl (10) CA
-% (11) CYrot (12) CNrot (13) CmRot (14) CnRot
-
-XCPvect = nan(7, 1);
-qualityIndex = zeros(4, 1);
-
-%% COEFFICIENTS
-if vars.flagXCP
- Cma = coeff( 6, iAlpha, 1, iBeta); % My D
- CNA = -coeff( 5, iAlpha, 1, iBeta); % Fz D
- Cm = coeff( 2, iAlpha, 1, iBeta); % My
- CN = -coeff( 1, iAlpha, 1, iBeta); % Fz
- Cnb = nan; % Mz D
- CYB = nan; % Fy D
- Cn = coeff( 4, iAlpha, 1, iBeta); % Mz
- CY = coeff( 3, iAlpha, 1, iBeta); % Fy
- Cl = nan;
- CA = nan;
-
-else
- Cma = coeff( 8, iAlpha, 1, iBeta); % My D
- CNA = -coeff( 4, iAlpha, 1, iBeta); % Fz D
- Cm = coeff( 9, iAlpha, 1, iBeta); % My
- CN = -coeff( 5, iAlpha, 1, iBeta); % Fz
- Cnb = coeff(12, iAlpha, 1, iBeta); % Mz D
- CYB = coeff( 2, iAlpha, 1, iBeta); % Fy D
- Cn = coeff(13, iAlpha, 1, iBeta); % Mz
- CY = coeff( 3, iAlpha, 1, iBeta); % Fy
- CA = -coeff(1, iAlpha, 1, iBeta);
- Cl = coeff(6, iAlpha, 1, iBeta);
-end
-
-Fyz = [0, CY, CN];
-Myz = [0, Cm, Cn];
-
-Faero = [CA, CY, CN];
-Maero = [Cl, Cm, Cn];
-vRyz = -1.*[0, tand(input.fltcon.BETA(iBeta)), tand(input.fltcon.ALPHA(iAlpha))];
-
-%% COMPUTE XCPs
-
-%%% LONGITUDINAL
-if abs(input.fltcon.ALPHA(iAlpha)) < 1
- XCPvect(1) = Cma/CNA;
-else
- XCPvect(1) = Cm/CN;
-end
-
-%%% LATERAL
-if input.fltcon.BETA(iBeta) == 0 && not(vars.useAlphaTot)
- XCPvect(2) = -Cnb/CYB;
-else
- XCPvect(2) = -Cn/CY;
-end
-
-%%% TOTAL FORMULATION
-
-Mdir = Myz./norm(Myz);
-
-Fpar = dot(Fyz, Mdir).*Mdir;
-Fnorm = Fyz - Fpar;
-
-%%% TOTAL - STD
-sgn = sign(dot(cross(Myz, Fyz), [1 0 0]));
-XCPvect(3) = sgn * norm(Myz)/norm(Fyz);
-
-%%% TOTAL - MOD
-sgn = sign(dot(cross(Myz, Fnorm), [1 0 0]));
-XCPvect(4) = sgn * norm(Myz)/norm(Fnorm);
-
-%%% ALPHA TOT - REFERENCE FRAME
-[~, phi] = getAlphaPhi(input.fltcon.ALPHA(iAlpha)*pi/180, input.fltcon.BETA(iBeta)*pi/180);
-R = [cos(phi) -sin(phi); sin(phi) cos(phi)];
-Frot = R*Fyz(2:3)';
-Mrot = R*Myz(2:3)';
-
-if norm(Frot) == 0
- XCPvect(5) = XCPvect(1); % Longitudinal
- XCPvect(6) = XCPvect(2); % Lateral
-else
- XCPvect(5) = Mrot(1)/Frot(2); % Longitudinal
- XCPvect(6) = -Mrot(2)/Frot(1); % Lateral
-end
-
-%%% OPEN ROCKET
-if vars.checkOpenRocket
- angle = getAlphaPhi(input.fltcon.ALPHA(iAlpha)*pi/180, input.fltcon.BETA(iBeta)*pi/180);
- XCPvect(7) = aerodynamics(input, angle, vars.modHaack);
-end
-%% SAVE COEFFICIENTS
-stabilityCoeffs = [Cma CNA Cm CN Cnb CYB Cn CY Cl CA Frot(1) Frot(2) Mrot(1) Mrot(2)]';
-
-%% COMPUTE QUALITY INDEXES
-qualityIndex(1) = acos( dot(Faero, Maero)/( norm(Faero)*norm(Maero) ) );
-qualityIndex(2) = acos( dot(Fyz, Myz)/ (norm(Fyz) * norm(Myz)) );
-qualityIndex(3) = acos( dot(vRyz, Fyz)/(norm(vRyz) * norm(Fyz)) );
-qualityIndex(4) = 100 - abs( (90 - acosd(dot(Myz, vRyz)/(norm(Myz) * norm(vRyz))) )/90 );
-qualityIndex(5) = norm(Fpar);
-end
\ No newline at end of file
diff --git a/stabilityAnalysis/src/stabilityPlots.m b/stabilityAnalysis/src/stabilityPlots.m
index 41c8962b..2b3be312 100644
--- a/stabilityAnalysis/src/stabilityPlots.m
+++ b/stabilityAnalysis/src/stabilityPlots.m
@@ -104,31 +104,31 @@ if settings.plots
else
warning('plotXCPlonATOT set to false')
end
- end
-
- if settings.plotQINDEX
- names = {'Angle between Faero and Maero', 'Angle between Fyz and Myz', 'Angle between Fyz and vRyz', 'T - index', 'Norm of Fpar'};
- for i = 1:5
- tit = strcat('Quality index num:', num2str(i));
- figure('Name',tit);
- qInd = squeeze(stabilityData.qIndexMat(i, :, :));
- if i <= 3
- qInd = qInd.*180/pi;
- end
- surf(settings.windMag, settings.windAz*180/pi, qInd', 'FaceColor', 'Interp');
- colormap jet;
- cb = colorbar;
- xlabel('Wind Magnitude (m/s)');
- ylabel('Wind Azimuth (deg)');
-
- if i <= 4
- cb.Label.String = 'qIndex (deg)';
- zlabel('qIndex (deg)');
- else
- cb.Label.String = 'qIndex (-)';
- zlabel('qIndex (-)');
+
+ if settings.plotQINDEX
+ names = {'Angle between Faero and Maero', 'Angle between Fyz and Myz', 'Angle between Fyz and vRyz', 'T - index', 'Norm of Fpar'};
+ for i = 1:5
+ tit = strcat('Quality index num:', num2str(i));
+ figure('Name',tit);
+ qInd = squeeze(stabilityData.qIndexMat(i, :, :));
+ if i <= 3
+ qInd = qInd.*180/pi;
+ end
+ surf(settings.windMag, settings.windAz*180/pi, qInd', 'FaceColor', 'Interp');
+ colormap jet;
+ cb = colorbar;
+ xlabel('Wind Magnitude (m/s)');
+ ylabel('Wind Azimuth (deg)');
+
+ if i <= 4
+ cb.Label.String = 'qIndex (deg)';
+ zlabel('qIndex (deg)');
+ else
+ cb.Label.String = 'qIndex (-)';
+ zlabel('qIndex (-)');
+ end
+ title(names{i})
end
- title(names{i})
end
end
end
--
GitLab