From cfaf29f1fc6b62834e810c20f2c9f2daaaab9aac Mon Sep 17 00:00:00 2001
From: giuliaghirardini <giuliaghirardini2001@gmail.com>
Date: Fri, 29 Nov 2024 16:14:20 +0100
Subject: [PATCH] [rocket-plot][classes] Implemented plot function in Rocket
class
---
classes/Rocket.m | 295 ++++++++++++++++++++++++++++++----------------
mainrocketplot.m | 21 ++--
provarocketplot.m | 7 ++
3 files changed, 214 insertions(+), 109 deletions(-)
create mode 100644 provarocketplot.m
diff --git a/classes/Rocket.m b/classes/Rocket.m
index 8a78fd6..8ae797c 100644
--- a/classes/Rocket.m
+++ b/classes/Rocket.m
@@ -9,18 +9,18 @@ classdef Rocket < Component
% - mission: Mission, mission object
% - varIn: (optional) config source. Alternative to config.m
% file
-
+
properties
- parafoil Parafoil % [-] Parafoil bay
+ parafoil Parafoil % [-] Parafoil bay
recovery Recovery % [-] Recovery bay
electronics Electronics % [-] Electronics bay
airbrakes Airbrakes % [-] Airbrakes bay
motor Motor % [-] Motor bay
rear Rear % [-] Rear bay
-
+
pitot Pitot % [-] Pitot component
parachutes Para % [-] (nParachutes, nStages) Parachutes onboard
-
+
length double % [m] Total length
lengthCenter double % [m] Center length - no nose, boat
lengthCenterNoMot double % [m] Center length - no nose, no motor
@@ -31,13 +31,13 @@ classdef Rocket < Component
inertiaNoMotor double % [kg*m^2] Inertia without motor
xCg double % [m] Total xCg
xCgNoMotor double % [m] xCg without motor
-
+
stagesMass double % [kg] Mass of stage 2, 3, .. without chutes
-
+
coefficients % [-] Aerodynamic coefficients
coefficientsHighAOA % [-] Aerodynamic coefficients for high angle of attacks
end
-
+
properties(SetAccess = protected)
crossSection % [m^2] Rocket cross sectional area
inertiaDot double % [kg*m^2/s]Total inertia time derivative
@@ -48,18 +48,18 @@ classdef Rocket < Component
bays dictionary ... % [Bay] All bays
= dictionary()
end
-
+
properties(Access = protected, Hidden)
configName = 'rocketConfig.m'
variableName = ''
mission Mission
end
-
+
methods % Updaters
function updateAbsolutePositions(obj)
bDictionay = obj.bays;
b = values(bDictionay);
-
+
if ~isempty(obj.lengthCenterNoMot)
obj.absolutePositions = dictionary(["motor", "rear"], ...
[obj.lengthCenterNoMot, obj.lengthCenterNoMot + obj.rear.position]);
@@ -70,15 +70,15 @@ classdef Rocket < Component
shift = zeros(size(b));
overrides = ~cellfun(@isempty, {b.position});
overridesShift = logical([overrides(2:end), 0]);
-
+
shift(overridesShift) = b(overrides).position;
shift(~overridesShift) = [b(~overridesShift).length];
-
+
absPositions = cumsum([-shift(1); shift(1:end-1)]);
absPositions = absPositions + ...
(obj.parafoil.length - obj.parafoil.noseLength);
obj.absolutePositions = dictionary(keys(bDictionay), absPositions);
-
+
end
function updateGeometry(obj)
if ~isempty(obj.lengthCenterNoMot)
@@ -92,7 +92,7 @@ classdef Rocket < Component
end
obj.lengthCenter = obj.length - obj.parafoil.noseLength - obj.rear.boatLength;
end
-
+
function updateMassNoMotor(obj)
if ~isempty(obj.massNoMotor)
return;
@@ -100,11 +100,11 @@ classdef Rocket < Component
bNoMotor = values(remove(obj.bays, 'motor'));
obj.massNoMotor = sum([bNoMotor.mass]);
end
-
+
function updateMass(obj)
obj.mass = obj.massNoMotor + obj.motor.mass;
end
-
+
function updateXCgNoMotor(obj)
if ~isempty(obj.xCgNoMotor)
return;
@@ -114,13 +114,13 @@ classdef Rocket < Component
obj.xCgNoMotor = ((aPosNoMotor + [bNoMotor.xCg])* ...
[bNoMotor.mass]')/obj.massNoMotor;
end
-
+
function updateXCg(obj)
obj.xCg = (obj.xCgNoMotor*obj.massNoMotor + ...
obj.motor.mass.*(obj.motor.xCg+obj.lengthCenterNoMot))./ ...
obj.mass;
end
-
+
function updateInertiaNoMotor(obj)
if ~isempty(obj.inertiaNoMotor)
return;
@@ -134,7 +134,7 @@ classdef Rocket < Component
(obj.xCgNoMotor' - (aPosNoMotor + [bNoMotor.xCg])).^2 .* [bNoMotor.mass];
obj.inertiaNoMotor = [sum(inertias(1, :)); sum(temp, 2)];
end
-
+
function updateInertia(obj)
% [3xN] Ix, Iy, Iz
% Assumption: xCgs are close to rocket axis
@@ -142,19 +142,19 @@ classdef Rocket < Component
(obj.xCg - (obj.absolutePositions('motor') + obj.motor.xCg)).^2 .* obj.motor.mass;
baysInertia = obj.inertiaNoMotor(2:3, :) + ...
(obj.xCg - obj.xCgNoMotor).^2 .* obj.massNoMotor;
-
+
obj.inertia = [obj.motor.inertia(1, :) + obj.inertiaNoMotor(1, :);
motorInertia + baysInertia];
-
+
obj.inertiaDot = diff(obj.inertia, 1, 2)./diff(obj.motor.time);
obj.inertiaDot(:, end + 1) = obj.inertiaDot(:, end);
end
-
+
function updateCutoff(obj)
obj.cutoffInertia = interpLinear(obj.motor.time, obj.inertia, obj.motor.cutoffTime);
obj.cutoffMass = interpLinear(obj.motor.time, obj.mass, obj.motor.cutoffTime);
end
-
+
function updateStagesMass(obj)
stage1 = obj.cutoffMass - (obj.parachutes(1,2).mass + obj.parachutes(2,2).mass + obj.parafoil.mass);
% Everything at cut off without parafoil, parafoil drogue and
@@ -163,20 +163,20 @@ classdef Rocket < Component
% only parafoil: parafoil, parafoil drogue and parafoil airfoil
obj.stagesMass = [stage1 stage2];
end
-
+
function updateAll(obj)
% Note: properties without motor must be updated first
% if obj.motor.isHRE
- obj.updateAbsolutePositions;
- obj.updateGeometry;
- obj.updateMassNoMotor;
- obj.updateMass;
- obj.updateXCgNoMotor;
- obj.updateXCg;
- obj.updateInertiaNoMotor;
- obj.updateInertia;
- obj.updateCutoff;
- obj.updateStagesMass;
+ obj.updateAbsolutePositions;
+ obj.updateGeometry;
+ obj.updateMassNoMotor;
+ obj.updateMass;
+ obj.updateXCgNoMotor;
+ obj.updateXCg;
+ obj.updateInertiaNoMotor;
+ obj.updateInertia;
+ obj.updateCutoff;
+ obj.updateStagesMass;
% end
% obj.updateAbsolutePositions;
% obj.updateGeometry;
@@ -187,30 +187,30 @@ classdef Rocket < Component
% obj.updateStagesMass;
end
end
-
+
methods(Access = protected) % Loaders
function obj = loadData(obj)
- if isempty(obj.motor)
+ if isempty(obj.motor)
return;
end
varNames = {'total', 'geometry', 'state', 'finsCN'};
motorName = obj.motor.name;
aeroPath = fullfile(obj.mission.dataPath, 'aeroCoefficients.mat');
aeroHighAOAPath = fullfile(obj.mission.dataPath, 'aeroCoefficientsHighAOA.mat');
-
+
if ~(exist(aeroPath, 'file') && exist(aeroHighAOAPath,'file'))
return;
end
dataCoeffs = load(aeroPath);
dataCoeffsHighAOA = load(aeroHighAOAPath);
-
+
if isfield(dataCoeffs, motorName) && isfield(dataCoeffsHighAOA, motorName)
dataCoeffs = dataCoeffs.(motorName);
dataCoeffsHighAOA = dataCoeffsHighAOA.(motorName);
else
return;
end
-
+
if isfield(dataCoeffs, "finsCN") && ...
isfield(dataCoeffsHighAOA, "finsCN")
obj.coefficients.finsCN = dataCoeffs.finsCN;
@@ -226,10 +226,10 @@ classdef Rocket < Component
obj.coefficientsHighAOA.geometry = dataCoeffsHighAOA.geometry;
obj.coefficients.state = dataCoeffs.state;
obj.coefficientsHighAOA.state = dataCoeffsHighAOA.state;
- end
+ end
end
end
-
+
methods
function obj = Rocket(mission, varIn, options)
arguments
@@ -239,12 +239,12 @@ classdef Rocket < Component
options.checkGeometry logical = true
end
obj@Component(mission, varIn);
-
+
%% Loading data
if isempty(mission.name)
return;
end
-
+
vars = obj.getConfig(mission); % Load config once and apply to other bays
obj.parafoil = Parafoil(mission, vars);
obj.recovery = Recovery(mission, vars);
@@ -255,11 +255,11 @@ classdef Rocket < Component
obj.pitot = Pitot(mission, vars);
paras = Para(mission);
obj.parachutes = paras(:, :, 1);
-
+
[bayNames, bayValues] = obj.getProperties(Superclass='Bay', Heterogeneous=1);
obj.bays = dictionary(bayNames, bayValues);
obj.updateAll();
-
+
if options.loadCoefficients
obj.loadData();
answer = '';
@@ -270,7 +270,7 @@ classdef Rocket < Component
answer = questdlg(['Coefficient matrices differ from rocket geometry. ' ...
'Do you want to create new matrices?']);
end
-
+
switch answer
case 'Yes'
parserPath = fullfile(mission.msaPath, 'autoMatricesProtub');
@@ -283,7 +283,7 @@ classdef Rocket < Component
end
end
end
-
+
function [coeffsValues, angle0] = interpCoeffs(obj, t, alpha, mach, beta, alt, c)
% interpCoeffs - interpolation of aerodynamic coefficients.
%
@@ -298,20 +298,20 @@ classdef Rocket < Component
% OUTPUTS:
% - coeffsValues, array [16,1], aerodynamic coefficients;
% - angle0, array [2,1], array of the reference aerodynamic angles.
-
+
if mach > obj.airbrakes.maxMach
c = 0;
end
h = c*obj.coefficients.state.hprot(end);
-
+
[nearAlpha, indexAlpha] = nearestValSorted(obj.coefficients.state.alphas*pi/180, alpha);
[~, indexMach] = nearestValSorted(obj.coefficients.state.machs, mach);
[nearBeta, indexBeta] = nearestValSorted(obj.coefficients.state.betas*pi/180, beta);
[~, indexAlt] = nearestValSorted(obj.coefficients.state.altitudes, alt);
[~, indexControl] = nearestValSorted(obj.coefficients.state.hprot, h);
-
+
angle0 = [nearAlpha; nearBeta];
-
+
if t >= obj.motor.cutoffTime
% Interpolate on airbrake control
% Uses end xCg, even if tCutOff < tBurnTime
@@ -327,10 +327,10 @@ classdef Rocket < Component
index1 = index1 - 1;
end
end
-
+
VF = obj.coefficients.total(:, indexAlpha, indexMach, indexBeta, indexAlt, index1, end);
VE = obj.coefficients.total(:, indexAlpha, indexMach, indexBeta, indexAlt, index2, end);
-
+
deltaH = obj.coefficients.state.hprot(index2) - obj.coefficients.state.hprot(index1);
coeffsValues = ( (h - obj.coefficients.state.hprot(index1))/deltaH)*(VE - VF) + VF;
else
@@ -339,25 +339,122 @@ classdef Rocket < Component
indexXcg = t >= obj.coefficients.state.xcgTime(1:end-1) & t < obj.coefficients.state.xcgTime(2:end);
index1 = num(indexXcg);
index2 = index1 + 1;
-
+
VF = obj.coefficients.total(:, indexAlpha, indexMach, indexBeta, indexAlt, indexControl, index1);
VE = obj.coefficients.total(:, indexAlpha, indexMach, indexBeta, indexAlt, indexControl, index2);
-
+
deltaT = obj.coefficients.state.xcgTime(index2) - obj.coefficients.state.xcgTime(index1);
coeffsValues = ( (t - obj.coefficients.state.xcgTime(index1))/deltaT)*(VE - VF) + VF;
end
end
-
- % function rocketPlot(obj)
- % deltaXLE = obj.rear.finsDeltaXFreeChord;
- % r = obj.diameter/2;
- % height = obj.rear.finsHeight;
- % C1 = obj.rear.finsRootChord;
- % C2 = obj.rear.finsFreeChord;
- %
- % %%% still developing. this is just a try
- % end
-
+
+ function img = rocketPlot(obj, mission)
+ deltaXLE = obj.rear.finsDeltaXFreeChord;
+ r = obj.diameter/2;
+ height = obj.rear.finsHeight;
+ C1 = obj.rear.finsRootChord;
+ C2 = obj.rear.finsFreeChord;
+
+ Lcent = obj.lengthCenter;
+ Lnos = obj.parafoil.noseLength;
+
+ Xle1 = Lcent + Lnos - obj.rear.finsAxialDistance - C1;
+
+ Daft = obj.rear.boatFinalDiameter;
+ Laft = obj.rear.boatLength;
+
+ xNos = Lnos - Lnos*cos(linspace(0, pi/2, 50));
+ if strcmp(obj.parafoil.noseType, 'KARMAN')
+ theta = @(x) acos( 1 - ( (2*x)./Lnos ) );
+ Karman = @(x) ( r/sqrt(pi) ) * sqrt( theta(x) - ( sin(2*theta(x))./ 2 ) );
+ yNos = Karman(xNos);
+ elseif strcmp(obj.parafoil.noseType, 'HAACK')
+ theta = @(x) acos( 1 - ( (2*x)./Lnos ) );
+ Haack = @(x) ( r/sqrt(pi) ) * sqrt( theta(x) - ( sin(2*theta(x))./ 2 ) + (1/3)*sin(theta(x)).^3 );
+ yNos = Haack(xNos);
+ elseif strcmp(obj.parafoil.noseType, 'OGIVE')
+ rho = (r^2 + Lnos^2)/(2*r);
+ Ogive = @(x) sqrt(rho^2 - (Lnos - x).^2) + r - rho;
+ yNos = Ogive(xNos);
+ elseif strcmp(obj.parafoil.noseType, 'POWER')
+ power = obj.parafoil.nosePower;
+ Power = @(x) r * (x/Lnos).^(power);
+ yNos = Power(xNos);
+ elseif strcmp(obj.parafoil.noseType, 'MHAACK')
+ p = obj.parafoil.nosePMod;
+ c = obj.parafoil.noseCMod;
+ xMod = @(x, p) (x/Lnos).^(p)*Lnos;
+ thetaMod = @(x, p) acos( 1 - ( (2*xMod(x, p))./Lnos ) );
+ haackSeriesMod = @(x, p, C) ( r/sqrt(pi) ) * sqrt( thetaMod(x, p) ...
+ - ( sin(2*thetaMod(x, p))./ 2 ) + C*sin(thetaMod(x, p)).^3 );
+ yNos = haackSeriesMod(xNos, p, c);
+ end
+
+ if strcmp(obj.rear.boatType, 'OGIVE') % [-] Boat-tail shape. 0: Cone, 1: Tangent Ogive
+ [xBoat, yBoat] = computeTangentBoatPoints(2*r, Daft, Laft);
+ else
+ xBoat = [0 Laft];
+ yBoat = [r Daft/2];
+ end
+
+ %%% figure begin
+ img = figure();
+ hold on
+
+ %%% NOSECONE
+ plot(xNos, yNos, 'k');
+ plot(xNos, -yNos, 'k');
+ plot([Lnos Lnos], [-r r], 'k');
+
+ %%% CENTERBODY
+ plot([Lnos Lnos+Lcent], [r r], 'k');
+ plot([Lnos Lnos+Lcent], [-r -r], 'k');
+ plot([Lnos+Lcent Lnos+Lcent], [-r r], 'k');
+
+ %%% BOAT-TAIL PLOT
+ plot(xBoat+Lnos+Lcent, yBoat, 'k');
+ plot(xBoat+Lnos+Lcent, -yBoat, 'k');
+ plot([Lnos+Lcent+Laft Lnos+Lcent+Laft], [-Daft/2 Daft/2], 'k');
+
+ %%% FINS PLOT
+ % top
+ plot([Xle1 Xle1+deltaXLE], [r r+height],'k');
+ plot([Xle1+C1 Xle1+deltaXLE+C2], [r r+height],'k');
+ plot([Xle1+deltaXLE Xle1+deltaXLE+C2], [r+height r+height],'k');
+ plot([Xle1 Xle1+C1], [r r], 'k');
+ % bottom
+ plot([Xle1 Xle1+deltaXLE], [-r -r-height],'k');
+ plot([Xle1+C1 Xle1+deltaXLE+C2], [-r -r-height],'k');
+ plot([Xle1+deltaXLE Xle1+deltaXLE+C2], [-r-height, -r-height],'k');
+ plot([Xle1 Xle1+C1], [-r -r], 'k');
+
+ % plot([Xle1 Xle1+deltaXLE], [r+height r+height], 'k--')
+ % plot([Xle1+deltaXLE Xle1+deltaXLE], [r+height r],'k--')
+
+ %%% BAYS
+ lengths = obj.absolutePositions.values + 2*Lnos;
+ baysBody = obj.absolutePositions.keys;
+ xline(0,'r-.', 'Label', 'NOSE', ...
+ 'LabelOrientation', 'aligned', 'LabelVerticalAlignment', 'bottom', ...
+ 'LabelHorizontalAlignment', 'right');
+ for i = 1:length(lengths)
+ xline(lengths(i), '--', 'Label', baysBody(i), ...
+ 'LabelOrientation', 'aligned', 'LabelVerticalAlignment', 'bottom', ...
+ 'LabelHorizontalAlignment', 'right');
+ end
+ boatPlotLength = obj.length - obj.rear.boatLength;
+ xline(boatPlotLength,'r-.', 'Label', 'BOAT-TAIL', ...
+ 'LabelOrientation', 'aligned', 'LabelVerticalAlignment', 'bottom', ...
+ 'LabelHorizontalAlignment', 'right');
+
+ axis equal
+ % set(gca, 'xlim', [-100, Lnos+Lcent+Laft+100]);
+
+ %%% TITLE
+ title(string(strrep(obj.mission.name, '_', ' ')));
+ end
+
+
function checks = checkGeometry(obj)
% checkGeometry - This methods checks if the rocket geometry
% is consistent with the geometry of the
@@ -365,12 +462,12 @@ classdef Rocket < Component
%
% OUTPUTS:
% - checks (n fields of geometry, 1): boolean value of the geometry checks
-
+
xCgRocket = round([ ...
obj.xCg(1); ...
obj.xCg(end) ...
], 3);
-
+
xCgTest = round([
obj.coefficients.geometry.xcg(1);
obj.coefficients.geometry.xcg(end);
@@ -389,7 +486,7 @@ classdef Rocket < Component
obj.rear.boatLength;
obj.rear.boatFinalDiameter;
], 3);
-
+
geometryTest = round([
obj.coefficients.geometry.diameter;
obj.coefficients.geometry.lNose;
@@ -402,38 +499,38 @@ classdef Rocket < Component
obj.coefficients.geometry.boatL;
obj.coefficients.geometry.boatD;
], 3);
- else % MHAAK ogive case, modified p and c coefficients
+ else % MHAAK ogive case, modified p and c coefficients
geometryRocket = round([
- obj.diameter;
- obj.parafoil.noseLength;
- obj.lengthCenter;
- obj.parafoil.noseCMod;
- obj.parafoil.nosePMod;
- obj.rear.finsRootChord;
- obj.rear.finsFreeChord;
- obj.rear.finsHeight;
- obj.rear.finsDeltaXFreeChord;
- obj.rear.nPanel;
- obj.rear.boatLength;
- obj.rear.boatFinalDiameter;
- ], 3);
-
- geometryTest = round([
- obj.coefficients.geometry.diameter;
- obj.coefficients.geometry.lNose;
- obj.coefficients.geometry.lCenter;
- obj.coefficients.geometry.cMod;
- obj.coefficients.geometry.pMod;
- obj.coefficients.geometry.chord1;
- obj.coefficients.geometry.chord2;
- obj.coefficients.geometry.height;
- obj.coefficients.geometry.deltaXLE;
- obj.coefficients.geometry.nPanel;
- obj.coefficients.geometry.boatL;
- obj.coefficients.geometry.boatD;
- ], 3);
+ obj.diameter;
+ obj.parafoil.noseLength;
+ obj.lengthCenter;
+ obj.parafoil.noseCMod;
+ obj.parafoil.nosePMod;
+ obj.rear.finsRootChord;
+ obj.rear.finsFreeChord;
+ obj.rear.finsHeight;
+ obj.rear.finsDeltaXFreeChord;
+ obj.rear.nPanel;
+ obj.rear.boatLength;
+ obj.rear.boatFinalDiameter;
+ ], 3);
+
+ geometryTest = round([
+ obj.coefficients.geometry.diameter;
+ obj.coefficients.geometry.lNose;
+ obj.coefficients.geometry.lCenter;
+ obj.coefficients.geometry.cMod;
+ obj.coefficients.geometry.pMod;
+ obj.coefficients.geometry.chord1;
+ obj.coefficients.geometry.chord2;
+ obj.coefficients.geometry.height;
+ obj.coefficients.geometry.deltaXLE;
+ obj.coefficients.geometry.nPanel;
+ obj.coefficients.geometry.boatL;
+ obj.coefficients.geometry.boatD;
+ ], 3);
end
-
+
checks = [
xCgRocket == xCgTest;
strcmp(obj.coefficients.geometry.ogType, obj.parafoil.noseType);
diff --git a/mainrocketplot.m b/mainrocketplot.m
index b3b6816..402cf94 100644
--- a/mainrocketplot.m
+++ b/mainrocketplot.m
@@ -18,13 +18,7 @@ r2 = rocketPlot(mission2, rocket2);
first_ax = findobj(r, 'type', 'axes');
second_ax = findobj(r2, 'type', 'axes');
-% Modify the appearance of the second rocket's plot
-ch2 = get(second_ax, 'children'); % Get children of second figure's axes
-for i = 1:length(ch2)
- if isprop(ch2(i), 'Color') % Check if the child has a 'Color' property
- set(ch2(i), 'Color', 'blue', 'LineStyle', '--'); % Set color to blue, dashed line
- end
-end
+
% Create a new figure for the merged content
merged_fig = figure('Name', 'Merged Figure');
@@ -35,17 +29,24 @@ ch1 = get(first_ax, 'children'); % Direct children only
copyobj(ch1, merged_ax); % Copy to merged axes
% Copy the modified children from the second figure's axes
+% Modify the appearance of the second rocket's plot
+ch2 = get(second_ax, 'children'); % Get children of second figure's axes
+
+for i = 1:length(ch2)
+ if isprop(ch2(i), 'Color') && ~isprop(ch2(i), 'Label') % Check if the child has a 'Color' property
+ set(ch2(i), 'Color', 'blue', 'LineStyle', '--', 'LineWidth', 1.2); % Set color to blue, dashed line
+ end
+end
copyobj(ch2, merged_ax); % Copy to merged axes
% Set up legend with correct colors
-legend_labels = [ch1(1), ch2(1)]; % Use the first plot from each figure for the legend
+legend_labels = [ch1(9), ch2(end)]; % Use the first plot from each figure for the legend
legend(merged_ax, legend_labels, ...
{string(strrep(mission.name, '_', ' ')), string(strrep(mission2.name, '_', ' '))}, ...
'Location', 'best');
% Adjust labels and appearance
-xlabel(merged_ax, 'Merged X-axis Label');
-ylabel(merged_ax, 'Merged Y-axis Label');
+xlabel(merged_ax, 'Length [m]');
title('Rocket Comparison');
axis equal
diff --git a/provarocketplot.m b/provarocketplot.m
new file mode 100644
index 0000000..6d9ca1d
--- /dev/null
+++ b/provarocketplot.m
@@ -0,0 +1,7 @@
+clear
+clc
+
+mission = Mission('design');
+rocket = Rocket(mission);
+
+gpx = rocket.rocketPlot(mission);
\ No newline at end of file
--
GitLab