diff --git a/classes/Rocket.m b/classes/Rocket.m index 0088c2cc31dc65486d83f765dca38ca44d082d77..5cb47ec18224f5503545361a85dabacb95613c1b 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,21 +70,22 @@ 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) obj.length = obj.lengthCenterNoMot + obj.motor.length + obj.parafoil.noseLength; else - obj.length = (obj.absolutePositions("rear") + obj.bays("rear").length) + obj.bays("parafoil").noseLength; + obj.length = (obj.absolutePositions("rear") + obj.bays("rear").length ... + + obj.bays("rear").finsDeltaXFreeChord) + obj.bays("parafoil").noseLength; obj.lengthCenterNoMot = obj.absolutePositions("motor"); end if isempty(obj.crossSection) @@ -92,7 +93,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 +101,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 +115,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 +135,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 +143,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 +164,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 +188,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 +227,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 +240,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 +256,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 +271,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 +284,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 +299,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 +328,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,15 +340,132 @@ 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 img = plot(obj) + arguments + obj Rocket + end + 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; + iMotor = find(strcmp(obj.absolutePositions.keys,"motor")); + xline(0,'r-.', 'Label', 'NOSE', ... + 'LabelOrientation', 'aligned', 'LabelVerticalAlignment', 'bottom', ... + 'LabelHorizontalAlignment', 'right'); + for i = 1:length(lengths) + if i == iMotor + line([lengths(i) lengths(i)], [-1 1], 'LineStyle', '--', 'Color', 'b'); % Draw the line + text(lengths(i), -1 + 0.3, baysBody(i), 'VerticalAlignment', 'top', ... + 'HorizontalAlignment', 'right', 'Rotation', 90, 'Color', 'b'); + else + line([lengths(i) lengths(i)], [-1 1], 'LineStyle', '--', 'Color', 'k'); % Draw the line + text(lengths(i), 1 + 0.1, baysBody(i), 'VerticalAlignment', 'top', ... + 'HorizontalAlignment', 'right', 'Rotation', 90); + end + 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 @@ -355,12 +473,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); @@ -379,7 +497,7 @@ classdef Rocket < Component obj.rear.boatLength; obj.rear.boatFinalDiameter; ], 3); - + geometryTest = round([ obj.coefficients.geometry.diameter; obj.coefficients.geometry.lNose; @@ -392,38 +510,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 new file mode 100644 index 0000000000000000000000000000000000000000..402cf94494293001bebc0111eec2b1917dd0ca31 --- /dev/null +++ b/mainrocketplot.m @@ -0,0 +1,53 @@ +close all +clear +clc + +%% +mission = Mission('design'); +mission2 = Mission('design'); + +rocket = Rocket(mission); +rocket2 = Rocket(mission2); + +r = rocketPlot(mission, rocket); +r2 = rocketPlot(mission2, rocket2); + +%% Merging Figures + +% Get axes from both figures +first_ax = findobj(r, 'type', 'axes'); +second_ax = findobj(r2, 'type', 'axes'); + + + +% Create a new figure for the merged content +merged_fig = figure('Name', 'Merged Figure'); +merged_ax = axes(merged_fig); % Create new axes in the merged figure + +% Copy the children from the first figure's axes +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(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, 'Length [m]'); +title('Rocket Comparison'); +axis equal + +disp('Figures merged successfully.'); diff --git a/provarocketplot.m b/provarocketplot.m new file mode 100644 index 0000000000000000000000000000000000000000..6d9ca1dfc1c9d25e8375586340ec77f33860e83a --- /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 diff --git a/rocketPlot.m b/rocketPlot.m new file mode 100644 index 0000000000000000000000000000000000000000..c95019cb7b293ee336a3c63ced49e62231ae0d53 --- /dev/null +++ b/rocketPlot.m @@ -0,0 +1,109 @@ +function roc = rocketPlot(mission, rocket) +arguments + mission Mission + rocket Rocket +end + deltaXLE = rocket.rear.finsDeltaXFreeChord; + r = rocket.diameter/2; + height = rocket.rear.finsHeight; + C1 = rocket.rear.finsRootChord; + C2 = rocket.rear.finsFreeChord; + + Lcent = rocket.lengthCenter; + Lnos = rocket.parafoil.noseLength; + + Xle1 = Lcent + Lnos - rocket.rear.finsAxialDistance - C1; + + Daft = rocket.rear.boatFinalDiameter; + Laft = rocket.rear.boatLength; + + xNos = Lnos - Lnos*cos(linspace(0, pi/2, 50)); + if strcmp(rocket.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(rocket.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(rocket.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(rocket.parafoil.noseType, 'POWER') + power = rocket.parafoil.nosePower; + Power = @(x) r * (x/Lnos).^(power); + yNos = Power(xNos); + elseif strcmp(rocket.parafoil.noseType, 'MHAACK') + p = rocket.parafoil.nosePMod; + c = rocket.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(rocket.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 + roc = 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 = rocket.absolutePositions.values + 2*Lnos; + bays = rocket.absolutePositions.keys; + xline(0,'r-.', 'Label', 'NOSE', ... + 'LabelOrientation', 'aligned', 'LabelVerticalAlignment', 'bottom', ... + 'LabelHorizontalAlignment', 'right'); + for i = 1:length(lengths) + xline(lengths(i), '--', 'Label', bays(i), ... + 'LabelOrientation', 'aligned', 'LabelVerticalAlignment', 'bottom', ... + 'LabelHorizontalAlignment', 'right'); + end + boatPlotLength = rocket.length - rocket.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(mission.name, '_', ' '))); +end