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