diff --git a/classes/DataWrapper.m b/classes/DataWrapper.m
index 62981eb589d4bc3a4cf7b0a49d416abd813dddfa..3fa54c678173b1ecbbce4cdee89e04f995674ec1 100644
--- a/classes/DataWrapper.m
+++ b/classes/DataWrapper.m
@@ -15,12 +15,12 @@ classdef DataWrapper < handle
% improve code maintainability
properties(SetAccess = private)
- data struct % Cached data
- dummy struct % Empty struct for field reference
+ data struct % Cached data
+ dummy struct % Empty struct for field reference
+ counter int32 % Current size of data stored
end
properties(Access = private)
- counter int32 % Current size of data stored
shift int32 % Number of fields to skip in cache, useful to preallocate initial values / starting conditions
chunkSize int32 % Size of preallocated memory
dataLen int32 % Current preallocated size of cache
diff --git a/functions/odeFunctions/ballistic.m b/functions/odeFunctions/ballistic.m
index 5992de85a73f1091cc44fbd481fec2190b6745e0..6c4a1f8857f7987097cedee58f0be5a18259ec03 100644
--- a/functions/odeFunctions/ballistic.m
+++ b/functions/odeFunctions/ballistic.m
@@ -350,6 +350,7 @@ if nargout == 2 || ~isempty(wrapper)
parout.coeff.XCPlat = XCPlat;
parout.coeff.XCPtot = XCPtot;
+ parout.uncertanty = [];
if ~isempty(wrapper)
wrapper.setCache(parout);
end
diff --git a/functions/odeFunctions/descentParachute.m b/functions/odeFunctions/descentParachute.m
index 6ec95378cec9be33c9d4dfb51b50b31bf024b771..556fbbdf136a2735edae5668f01ca0bcea716055 100644
--- a/functions/odeFunctions/descentParachute.m
+++ b/functions/odeFunctions/descentParachute.m
@@ -150,6 +150,7 @@ if nargout == 2 || ~isempty(wrapper)
parout.accelerations.body = [du; dv; dw];
parout.coeff = [];
+ parout.uncertanty = [];
if ~isempty(wrapper)
wrapper.setCache(parout);
diff --git a/functions/odeFunctions/descentParafoil.m b/functions/odeFunctions/descentParafoil.m
index 66ceca02f3c4575132657618bcdb30fa604ce276..279a78574691ba2555eaba5debe418f2cfd445ee 100644
--- a/functions/odeFunctions/descentParafoil.m
+++ b/functions/odeFunctions/descentParafoil.m
@@ -227,6 +227,7 @@ if nargout == 2 || ~isempty(wrapper)
parout.accelerations.body = bodyAcc;
parout.coeff = [];
+ parout.uncertanty = [];
if ~isempty(wrapper)
wrapper.setCache(parout);
diff --git a/functions/simulations/quickApogeeOnly.m b/functions/simulations/OLD/quickApogeeOnly.m
similarity index 100%
rename from functions/simulations/quickApogeeOnly.m
rename to functions/simulations/OLD/quickApogeeOnly.m
diff --git a/functions/simulations/stochAscentLoop.m b/functions/simulations/OLD/stochAscentLoop.m
similarity index 100%
rename from functions/simulations/stochAscentLoop.m
rename to functions/simulations/OLD/stochAscentLoop.m
diff --git a/functions/simulations/stochDescentBalLoop.m b/functions/simulations/OLD/stochDescentBalLoop.m
similarity index 100%
rename from functions/simulations/stochDescentBalLoop.m
rename to functions/simulations/OLD/stochDescentBalLoop.m
diff --git a/functions/simulations/stochDescentLoop.m b/functions/simulations/OLD/stochDescentLoop.m
similarity index 100%
rename from functions/simulations/stochDescentLoop.m
rename to functions/simulations/OLD/stochDescentLoop.m
diff --git a/functions/simulations/stochPostProcess.m b/functions/simulations/OLD/stochPostProcess.m
similarity index 100%
rename from functions/simulations/stochPostProcess.m
rename to functions/simulations/OLD/stochPostProcess.m
diff --git a/functions/simulations/stochStabilityLoop.m b/functions/simulations/OLD/stochStabilityLoop.m
similarity index 100%
rename from functions/simulations/stochStabilityLoop.m
rename to functions/simulations/OLD/stochStabilityLoop.m
diff --git a/functions/simulations/stochWindData.m b/functions/simulations/OLD/stochWindData.m
similarity index 100%
rename from functions/simulations/stochWindData.m
rename to functions/simulations/OLD/stochWindData.m
diff --git a/functions/simulations/stdStability.m b/functions/simulations/stdStability.m
new file mode 100644
index 0000000000000000000000000000000000000000..7ff4dfb59d74b77a02194258d43615b19fd6c82b
--- /dev/null
+++ b/functions/simulations/stdStability.m
@@ -0,0 +1,128 @@
+function stability = stdStability(Y0, rocket, env, wind, settings)
+% stochStabilityLoop - This function runs a simulation to compute the rocket
+% stability inside the parfor loop
+%
+% INPUTS:
+% - settings, struct: stores data of the rocket and of the simulation.
+% - Y0, double [13, 1]: initial state vector of the integration
+% [x, y, z, u, v, w, p, q, r, q0, q1, q2, q3]
+%
+% OUTPUTS:
+% - data_stability,cell [1, k]: cell containg the struct with the rocket
+% stability data at the launchpad exit
+%
+% CALLED FUNCTIONS:
+%
+% REVISIONS:
+% - 0 14/06/2023, Release, Riccardo Cadamuro
+%
+% Copyright © 2023, Skyward Experimental Rocketry, AFD department
+% All rights reserved
+%
+% SPDX-License-Identifier: GPL-3.0-or-later
+
+if rocket.airbrakes.enabled && length(rocket.airbrakes.extension) > 1 && ...
+ length(rocket.airbrakes.deltaTime) < length(rocket.airbrakes.extension) - 1
+ error('In airbrakes smooth opening simulations, AB configuration usage time vector must be at least of length length(airbrakes.extension)-1, check rocketConfig.m')
+end
+
+if any(rocket.airbrakes.extension > 3) || any(rocket.airbrakes.extension < 1)
+ error('At least a value in airbrakes.extension is set incorrectly in rocketConfig, they must be between 1 and 3')
+end
+
+%% ASCENT
+t0 = 0;
+tf = 10;
+solution = ode113(@ballistic, [t0, tf], Y0, options, ...
+ rocket, env, wind);
+
+if nargout == 2
+ t = solution.x;
+ Y = solution.y;
+
+ [~, stability] = ballistic(t(end), Y(:, end), rocket, env, wind);
+
+ ascent = getOdeFcn(solution);
+ ascent.state.Y = Y;
+ ascent.state.time = t;
+end
+
+%%
+
+%% comupute launchpad run
+[Tpad, Ypad] = ode113(@ascent, [0, 10], Y0, settings.ode.optionspad, settings);
+
+z = -Ypad(end, 3);
+[~, a, ~, ~] = atmosisa(settings.z0 + z);
+
+%% retrive the correct wind data
+if settings.wind.model
+ Hour = settings.stoch.Hour;
+ Day = settings.stoch.Day;
+ t = Tpad(end);
+
+ if settings.stoch.N > 1
+ [uw, vw, ww] = windMatlabGenerator(settings, z, t, Hour, Day);
+ else
+ [uw, vw, ww] = windMatlabGenerator(settings, z, t);
+ end
+
+elseif settings.wind.input
+ [uw, vw, ww] = windInputGenerator(settings, z, uncert);
+elseif settings.wind.variable
+
+ [uw, vw, ww] = windVariableGenerator(z, settings.stoch);
+else
+ uw = settings.stoch.uw; vw = settings.stoch.vw; ww = settings.stoch.ww;
+end
+inertialWind = [uw, vw, ww];
+
+%% compute the exit conditions in the aerodynamic frame
+vExit = norm(Ypad(end, 4));
+machExit = vExit/a;
+
+bodyWind = quatrotate(Y0(10:13)', inertialWind);
+bodyVelocity = [Ypad(end, 4), 0, 0];
+Vr = bodyVelocity - bodyWind;
+ur = Vr(1); vr = Vr(2); wr = Vr(3);
+
+if ur > 1e-9
+ alpha = atan(wr/ur)*180/pi;
+ beta = atan(vr/ur)*180/pi;
+else
+ alpha = 0;
+ beta = 0;
+end
+
+%% create dissileMatcom input data struct
+dissileVars.Alpha = sort([-alpha, alpha, 0, -1, 1]);
+dissileVars.Beta = beta;
+dissileVars.Alt = settings.z0 + z;
+dissileVars.Mach = machExit;
+dissileVars.xcg = interpLinear(settings.motor.expTime, settings.xcg, Tpad(end));
+dissileVars.hprot = [];
+input = createDissileInput(settings, dissileVars);
+
+input.flagXCP = true;
+coeff = dissileMatcom(input);
+
+vars.useAlphaTot = false;
+vars.checkOpenRocket = false;
+vars.modHaack = true;
+
+if dissileVars.Alpha(1) == alpha
+ iAlpha = 1;
+else
+ iAlpha = 5;
+end
+
+%% compute the stability margin and coresponding quality indices
+[XCPvect, qualityIndex, stabilityCoeffs] = getXCP(coeff, iAlpha, 1, input, vars);
+stability.time = Tpad;
+stability.state = Ypad;
+stability.wind = [uw, vw, ww];
+stability.XCPvect = XCPvect;
+stability.qualityIndex = qualityIndex;
+stability.stabilityCoeffs = stabilityCoeffs;
+
+end
\ No newline at end of file
diff --git a/settings/odeConfig.m b/settings/odeConfig.m
index 031e3512587a7466740006fef8ba1858497723d3..8e1fb785d62abce6dd04d17bf068ddece5c78dc7 100644
--- a/settings/odeConfig.m
+++ b/settings/odeConfig.m
@@ -3,10 +3,9 @@
% ODE settings
ode.finalTime = 2000; % [s] Final integration time
ode.optAscent = odeset('Events', @eventApogee, 'InitialStep', 1); % ODE options for ascend
-ode.optAscentDelayPara= odeset('InitialStep', 1); % ODE options for due to the opening delay of the parachute
ode.optParachute = odeset('Events', @eventParaCut, 'InitialStep', 1); % ODE options for the parachutes
ode.optDescent = odeset('Events', @eventLanding,...
-'RelTol', 1e-3, 'AbsTol', 1e-3, 'OutputFcn', @storeData); % ODE options for ballistic descent
+'RelTol', 1e-3, 'AbsTol', 1e-3); % ODE options for ballistic descent
ode.optLaunchpad = odeset('Events', @eventPad); % ODE options for the launchpad phase
% Settings for descent 6dof simulation