diff --git a/classes/@Rocket/Rocket.m b/classes/@Rocket/Rocket.m
index 8da0e3c0bbaafc4d47a86fbf3cecd820b2ac5003..5d91e4b504e335ecff22dacd88b6b63d839e856c 100644
--- a/classes/@Rocket/Rocket.m
+++ b/classes/@Rocket/Rocket.m
@@ -18,7 +18,8 @@ classdef Rocket < Config
         cutoffInertia       (3, 1)  double          % [kg*m^2]  Total inertia at motor cutoff
         inertiaDot          (3, :)  double          % [kg*m^2/s]Total inertia time derivative
 
-        stagesMass                  double          % [kg]      Mass of stage 2, 3, .. without chutes
+        stagesMass          (1, :)  double          % [kg]      Mass of stage 2, 3, .. without chutes
+        stagesInertia       (3, :)  double          % [kg*m^2]  Inertia of stage 2, 3, .. without chutes
         crossSection        (1, 1)  double          % [m^2]     Rocket cross sectional area
     end
 
@@ -101,6 +102,7 @@ classdef Rocket < Config
                 options.checkGeometry       logical = true
             end
             obj@Config(mission, varsIn);
+            obj.mission = mission;
             
             %% Initialize base properties
             bayNames = camel2UpperSnake(obj.getProperties(Superclass='Bay'));
diff --git a/classes/@Rocket/update.m b/classes/@Rocket/update.m
index 4400d73971f1171aa8da36f1e1899ddd58a1c526..ae246b7ba4aef47fe361667779210b34848682a5 100644
--- a/classes/@Rocket/update.m
+++ b/classes/@Rocket/update.m
@@ -38,8 +38,21 @@ if ~override('INERTIA')
 end
 
 %% Cutoff
+cutoffXcg = interp1(obj.time, obj.xcg, obj.cutoffTime, "linear", 'extrap');
 obj.cutoffMass = interp1(obj.time, obj.mass, obj.cutoffTime, "linear", 'extrap');
 obj.cutoffInertia = interp1(obj.time, obj.inertia', obj.cutoffTime, "linear", 'extrap')';
 
-end
+%% Stages
+% TODO! generalize for n stages, variable input
+obj.stagesMass = zeros(1, 2);
+obj.stagesMass(1) = obj.cutoffMass - obj.PARAFOIL.mass; % MAIN mass
+obj.stagesMass(2) = obj.PARAFOIL.mass;
+
+IParafoil = obj.PARAFOIL.inertia + ...
+    (cutoffXcg - (obj.PARAFOIL.xcg + obj.PARAFOIL.position))^2 * [0, 1, 1] * ...
+    obj.PARAFOIL.mass;
 
+obj.stagesInertia = zeros(3, 2);
+obj.stagesInertia(:, 1) = obj.cutoffInertia - IParafoil'; % MAIN mass
+obj.stagesInertia(:, 2) = obj.PARAFOIL.inertia';
+end
\ No newline at end of file
diff --git a/functions/eventFunctions/eventParaCut.m b/functions/eventFunctions/eventParaCut.m
index 7e7d3be656aeab0fea6d4ffaad2d0d1872040bba..211813d5a3a9516709bd727bfe751bd7c8393aa9 100644
--- a/functions/eventFunctions/eventParaCut.m
+++ b/functions/eventFunctions/eventParaCut.m
@@ -40,7 +40,7 @@ altitude = -Y(3);
 para = descentData.para;
 stage = descentData.stage;
 
-value = altitude - rocket.parachutes(para, stage).finalAltitude;
+value = altitude - rocket.parachutes{para, stage}.finalAltitude;
 
 isterminal = 1;
 direction = 0;
diff --git a/functions/odeFunctions/descentParachute.m b/functions/odeFunctions/descentParachute.m
index bbca96a724686eee8529a72ff6163447a23bea1c..47beab92eac2b5ffad9759890298be2472815807 100644
--- a/functions/odeFunctions/descentParachute.m
+++ b/functions/odeFunctions/descentParachute.m
@@ -47,9 +47,9 @@ dY = zeros(6, 1);
 para = descentData.para; % defined in stdRun {para = i; settings.paraNumber = para;}
 stage = descentData.stage;
 
-S = rocket.parachutes(para, stage).surface;          % [m^2]   Surface
-CD = rocket.parachutes(para, stage).cd;        % [/] Parachute Drag Coefficient
-CL = rocket.parachutes(para, stage).cl;        % [/] Parachute Lift Coefficient
+S = rocket.parachutes{para, stage}.surface;          % [m^2]   Surface
+CD = rocket.parachutes{para, stage}.cd;        % [/] Parachute Drag Coefficient
+CL = rocket.parachutes{para, stage}.cl;        % [/] Parachute Lift Coefficient
 
 m = rocket.stagesMass(stage);
 
diff --git a/functions/odeFunctions/descentParafoil.m b/functions/odeFunctions/descentParafoil.m
index ec73f328509b77df5b25b3c5dae785e0eedfe77a..3b40dac04bdfb91727c16a8e3893c9a02f58fd36 100644
--- a/functions/odeFunctions/descentParafoil.m
+++ b/functions/odeFunctions/descentParafoil.m
@@ -50,7 +50,7 @@ qz = Y(13);
 deltaA = Y(14);
 
 dY = zeros(13, 1);
-parafoil = rocket.parachutes(descentData.para, descentData.stage);
+parafoil = rocket.parachutes{descentData.para, descentData.stage};
 
 %% GETTING DATA
 % environment
diff --git a/functions/simulations/stdAscent.m b/functions/simulations/stdAscent.m
index f0ef8d5130612b979bae233f9e4d58f79a44ea01..60f9812a90d023f7253bcf17a52f34b434340bfa 100644
--- a/functions/simulations/stdAscent.m
+++ b/functions/simulations/stdAscent.m
@@ -54,8 +54,8 @@ t0 = 0;
 solution = ode113(@ballistic, [t0, tf], Y0, options, ...
     rocket, env, wind, [], 0, wrapper);
 
-if settings.simulator.parachute && ~isempty(rocket.parachutes(1, 1).openingTime) && rocket.parachutes(1, 1).openingTime
-    solution = odextend(solution, [], solution.x(end) + rocket.parachutes(1, 1).openingTime);
+if settings.simulator.parachute && ~isempty(rocket.parachutes{1, 1}.openingTime) && rocket.parachutes{1, 1}.openingTime
+    solution = odextend(solution, [], solution.x(end) + rocket.parachutes{1, 1}.openingTime);
 end
 
 if nargout == 2
diff --git a/functions/simulations/stdDescent.m b/functions/simulations/stdDescent.m
index 10076a50420c1eebfeb8d08f78897831964e4192..b4b077e94c552e33f26340d42972d130946592ab 100644
--- a/functions/simulations/stdDescent.m
+++ b/functions/simulations/stdDescent.m
@@ -107,13 +107,13 @@ for i = 1:nParachutes
     descentData.para = i;
     for j = 1:nStages
         descentData.stage = j;        
-        if isempty(rocket.parachutes(i, j).name), continue; end
+        if isempty(rocket.parachutes{i, j}.name), continue; end
         if i ~= 1 && ~isempty(descent(i-1, j).state)
             prev = descent(i-1, j); 
         end
         wrapper.reset();
 
-        switch class(rocket.parachutes(i, j))
+        switch class(rocket.parachutes{i, j})
             case 'Parachute'
                 Y0 = [prev.state.Y(1:3, end);
                     prev.state.Y(4:6, end)];
@@ -144,14 +144,14 @@ for i = 1:nParachutes
                     rocket, environment, wind, descentData, [], 0, wrapper);
             otherwise
                 error('Unrecognized para type for descent: %s', ...
-                    class(rocket.parachutes(i, j)))
+                    class(rocket.parachutes{i, j}))
         end
 
-        if i < nParachutes && rocket.parachutes(i+1, j).openingTime
-            solution = odextend(solution, [], solution.x(end) + rocket.parachutes(1, 1).openingTime);
+        if i < nParachutes && rocket.parachutes{i+1, j}.openingTime
+            solution = odextend(solution, [], solution.x(end) + rocket.parachutes{1, 1}.openingTime);
         end
 
-        if i < nParachutes && -solution.y(3, end) < rocket.parachutes(i+1, j).finalAltitude
+        if i < nParachutes && -solution.y(3, end) < rocket.parachutes{i+1, j}.finalAltitude
             error("Chute (%d, %d) has arrived at a lower altitude than Chute (%d, %d) target altitude." + ...
                 "   Please lower Chute (%d, %d) opening delay or target altitude.", i, j, i+1, j);
         end