From 83e33f51a27794e0973b0c1e08b6072c5d9fa264 Mon Sep 17 00:00:00 2001
From: Mauco03 <marco.gaibotti@skywarder.eu>
Date: Thu, 23 May 2024 12:21:51 +0200
Subject: [PATCH] [refactoring-ode][functions] Stramlined descent functions
---
functions/eventFunctions/eventPad.m | 53 +++++------
functions/simulations/stdAscent.m | 9 --
functions/simulations/stdDescentBal.m | 15 +++-
functions/simulations/stdStability.m | 124 +++++---------------------
settings/odeConfig.m | 4 +-
5 files changed, 61 insertions(+), 144 deletions(-)
diff --git a/functions/eventFunctions/eventPad.m b/functions/eventFunctions/eventPad.m
index 6245c7b..9048dc8 100644
--- a/functions/eventFunctions/eventPad.m
+++ b/functions/eventFunctions/eventPad.m
@@ -1,36 +1,31 @@
-function [value, isterminal, direction] = eventPad(~, Y, ~, environment, ~, ~)
-%{
-eventPad - Event function to stop simulation at the launchpad exit
-
-INPUTS:
-- t, double [1, 1], integration time [s];
-- Y, double [4, 1], integration state, check launchPadFreeDyn for explanation;
-- settings, struct (motor, CoeffsE, CoeffsF, para, ode, stoch, prob, wind),
- simulation data;
-- varargin, cell {3, 1}, for additional inputs of the ode function, un-necessary here.
-
-OUTPUTS:
-- value, double [1, 1], see eventFunction explantion in matlab website
-- isterminal, double [1, 1], see eventFunction explantion in matlab website
-- direction, double [1, 1], see eventFunction explantion in matlab website
-
-CALLED FUNCTIONS: /
-
-REVISIONS:
-- 0 21/10/20, release Adriano Filippo Inno
-Copyright © 2021, Skyward Experimental Rocketry, AFD department
-All rights reserved
-
-SPDX-License-Identifier: GPL-3.0-or-later
-
-%}
+function [value, isterminal, direction] = eventPad(~, Y, ~, environment, varargin)
+% eventPad - Event function to stop simulation at the launchpad exit
+%
+% INPUTS:
+% - t, double [1, 1], integration time [s];
+% - Y, double [4, 1], integration state, check launchPadFreeDyn for explanation;
+% - settings, struct (motor, CoeffsE, CoeffsF, para, ode, stoch, prob, wind),
+% simulation data;
+% - varargin, cell {3, 1}, for additional inputs of the ode function, un-necessary here.
+%
+% OUTPUTS:
+% - value, double [1, 1], see eventFunction explantion in matlab website
+% - isterminal, double [1, 1], see eventFunction explantion in matlab website
+% - direction, double [1, 1], see eventFunction explantion in matlab website
+%
+% CALLED FUNCTIONS: /
+%
+% REVISIONS:
+% - 0 21/10/20, release Adriano Filippo Inno
+% Copyright © 2021, Skyward Experimental Rocketry, AFD department
+% All rights reserved
+%
+% SPDX-License-Identifier: GPL-3.0-or-later
altitude = -Y(3);
-value = altitude - enviroment.rampLength*sin(environment.omega);
+value = altitude - environment.effectiveRampAltitude;
isterminal = 1;
direction = 1;
-
-
end
diff --git a/functions/simulations/stdAscent.m b/functions/simulations/stdAscent.m
index 0035cae..520964b 100644
--- a/functions/simulations/stdAscent.m
+++ b/functions/simulations/stdAscent.m
@@ -34,15 +34,6 @@ end
%
% 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
-
%% DATA PREPARATION
options = odeset(settings.ode.optAscent, 'OutputFcn', @storeData);
tf = settings.ode.finalTime;
diff --git a/functions/simulations/stdDescentBal.m b/functions/simulations/stdDescentBal.m
index 9655db1..d55fcac 100644
--- a/functions/simulations/stdDescentBal.m
+++ b/functions/simulations/stdDescentBal.m
@@ -36,11 +36,18 @@ SPDX-License-Identifier: GPL-3.0-or-later
%% DATA PREPARATION
options = odeset(settings.ode.optDescent, 'OutputFcn', @storeData);
-%% DESCEND
-solution = odextend(ascentSolution, [], settings.ode.finalTime, [], options);
+fun = ascentSolution.extdata.odefun;
+t0 = ascentSolution.x(end);
+tf = settings.ode.finalTime;
+
+Y0 = ascentSolution.y(:, end);
-apogeeIdx = find(solution.x == solution.xe(1), 1, "first");
-descent = getOdeFcn(solution, apogeeIdx);
+wrapper = ascentSolution.extdata.varargin{end};
+wrapper.reset();
+
+%% DESCEND
+solution = ode113(fun, [t0, tf], Y0, options, ascentSolution.extdata.varargin{:});
+descent = getOdeFcn(solution);
% if not(isfield(settings, 'unitTest')) || settings.unitTest == false
% save('descent_plot.mat', 'descent');
diff --git a/functions/simulations/stdStability.m b/functions/simulations/stdStability.m
index 7ff4dfb..cd315f4 100644
--- a/functions/simulations/stdStability.m
+++ b/functions/simulations/stdStability.m
@@ -1,4 +1,11 @@
-function stability = stdStability(Y0, rocket, env, wind, settings)
+function stability = stdStability(rocket, env, wind, settings, wrapper)
+arguments
+ rocket Rocket
+ env Environment
+ wind WindCustom {mustBeA(wind, {'WindCustom', 'WindMatlab'})}
+ settings Settings
+ wrapper DataWrapper
+end
% stochStabilityLoop - This function runs a simulation to compute the rocket
% stability inside the parfor loop
%
@@ -21,108 +28,25 @@ function stability = stdStability(Y0, rocket, env, wind, settings)
%
% 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);
+%% DATA PREPARATION
+options = odeset(settings.ode.optLaunchpad, 'OutputFcn', @storeData);
+tf = settings.ode.finalTime;
-%% retrive the correct wind data
-if settings.wind.model
- Hour = settings.stoch.Hour;
- Day = settings.stoch.Day;
- t = Tpad(end);
+%% INITIAL CONDITIONS
+%%% Attitude
+Q0 = angleToQuat(env.phi, env.omega, 180*pi/180)';
- if settings.stoch.N > 1
- [uw, vw, ww] = windMatlabGenerator(settings, z, t, Hour, Day);
- else
- [uw, vw, ww] = windMatlabGenerator(settings, z, t);
- end
+%%% State
+X0 = [0; 0; 0];
+V0 = [0; 0; 0];
+W0 = [0; 0; 0];
-elseif settings.wind.input
- [uw, vw, ww] = windInputGenerator(settings, z, uncert);
-elseif settings.wind.variable
+Y0 = [X0; V0; W0; Q0; 0];
- [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;
+%% LAUNCHPAD
+t0 = 0;
+solution = ode113(@ballistic, [t0, tf], Y0, options, ...
+ rocket, env, wind, [], 0, wrapper);
+stability = getOdeFcn(solution);
end
\ No newline at end of file
diff --git a/settings/odeConfig.m b/settings/odeConfig.m
index 8e1fb78..d1b5251 100644
--- a/settings/odeConfig.m
+++ b/settings/odeConfig.m
@@ -2,11 +2,11 @@
% ODE settings
ode.finalTime = 2000; % [s] Final integration time
-ode.optAscent = odeset('Events', @eventApogee, 'InitialStep', 1); % ODE options for ascend
+ode.optLaunchpad = odeset('Events', @eventPad, 'InitialStep', 1); % ODE options for the launchpad phase
+ode.optAscent = odeset('Events', @eventApogee, 'InitialStep', 1); % ODE options for ascent
ode.optParachute = odeset('Events', @eventParaCut, 'InitialStep', 1); % ODE options for the parachutes
ode.optDescent = odeset('Events', @eventLanding,...
'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
ode.optDrogue6DOF = odeset('Events', @eventParaCut, 'AbsTol', 1e-6,'RelTol', 1e-6); %ODE options for due to cutting of the drogue chute
--
GitLab