From 36a3e03be3cbce512d308fe0f66b750f4c07aa7a Mon Sep 17 00:00:00 2001 From: Mauco03 <marco.gaibotti@skywarder.eu> Date: Fri, 3 May 2024 22:52:56 +0200 Subject: [PATCH] [refactoring][functions] Testing new recallOde --- classes/DataWrapper.m | 94 +++++++++++++++++ functions/eventFunctions/eventApogee.m | 2 +- functions/miscellaneous/recallOdeFcn.m | 139 +++++++------------------ functions/odeFunctions/ballistic.m | 21 ++-- functions/odeOutputFcn/storeData.m | 10 ++ settings/odeConfig.m | 8 +- 6 files changed, 160 insertions(+), 114 deletions(-) create mode 100644 classes/DataWrapper.m create mode 100644 functions/odeOutputFcn/storeData.m diff --git a/classes/DataWrapper.m b/classes/DataWrapper.m new file mode 100644 index 0000000..04c8273 --- /dev/null +++ b/classes/DataWrapper.m @@ -0,0 +1,94 @@ +classdef DataWrapper < handle + % OutputWrapper: A fast way to exchange data by pointer + % + % Constructor: + % - OutputWrapper: Creates an instance of the Wrapper class. + % + % WARNING: Intended for data exchange where output arguments are inaccessible. + + properties(SetAccess = private) + data struct + end + + properties(Access = private) + cache + counter int32 = 0 + dummyStrc struct + chunkSize int32 % Size of preallocated memory + dataLen int32 % Current size of cache + end + + methods + function obj = DataWrapper(dataDummy, chunkSize) + arguments + dataDummy struct + chunkSize {mustBeInteger, mustBePositive} = 5 + end + obj.chunkSize = chunkSize; + + fields = fieldnames(dataDummy)'; fields{2, 1} = []; + obj.dummyStrc = struct(fields{:}); + + obj.cache = obj.dummyStrc; + obj.data = obj.dummyStrc; % Set struct dataType + obj.data(obj.chunkSize) = obj.dummyStrc; % Preallocate struct + + obj.dataLen = obj.chunkSize; + end + + function setCache(obj, data) + obj.cache = data; + end + + function setData(obj) + obj.counter = obj.counter + 1; + if obj.counter > obj.dataLen + len = obj.dataLen + obj.chunkSize; + obj.dataLen = len; + obj.data(len) = obj.dummyStrc; + end + obj.data(obj.counter) = obj.cache; + end + + function out = popData(obj) + out = obj.packageData(obj.data(1:obj.counter)); + obj.reset(); + end + + function out = getData(obj) + out = obj.data(1:obj.counter); + end + + function obj = reset(obj) + obj.cache = obj.dummyStrc; + obj.data = obj.dummyStrc; % Set struct dataType + obj.data(obj.chunkSize) = obj.dummyStrc; % Preallocate struct + + obj.dataLen = obj.chunkSize; + obj.counter = 0; + end + end + + methods(Static) + function out = packageData(data) + % Converts struct array to struct with array fields + % + % Main assumption: fields to concatenate are colums vectors or + % scalars + + fields = fieldnames(data); + + for i = 1:length(fields) + currentField = fields{i}; + currentValue = [data.(currentField)]; + + if isstruct(currentValue) + out.(currentField) = DataWrapper.packageData([data.(currentField)]); + else + out.(currentField) = currentValue; + end + end + end + end +end + diff --git a/functions/eventFunctions/eventApogee.m b/functions/eventFunctions/eventApogee.m index a8a8116..e586a95 100644 --- a/functions/eventFunctions/eventApogee.m +++ b/functions/eventFunctions/eventApogee.m @@ -1,4 +1,4 @@ -function [value, isterminal, direction] = eventApogee(t, Y, rocket, ~, ~, ~) +function [value, isterminal, direction] = eventApogee(t, Y, rocket, varargin) %{ eventApogee - Event function to stop simulation at apogee checking when a value tends to zero; the value taken is to account is the diff --git a/functions/miscellaneous/recallOdeFcn.m b/functions/miscellaneous/recallOdeFcn.m index 1a69a51..6af4a4e 100644 --- a/functions/miscellaneous/recallOdeFcn.m +++ b/functions/miscellaneous/recallOdeFcn.m @@ -1,100 +1,39 @@ -function [allSteps] = recallOdeFcn(fun, T, Y, varargin) -%{ -recallOdeFcn - This function allows to compute some parameters used - inside the ODE integrations. - -INPUTS: - - fun, function, ode function used; - - T, double [n° variation, 1], integration time vector; - - Y, double [n° variation, 16], state matrix, - [ x y z | u v w | p q r | q0 q1 q2 q3 | m | Ixx Iyy Izz ]: - * (x y z), NED{north, east, down} horizontal frame; - * (u v w), body frame velocities; - * (p q r), body frame angular rates; - * m , total mass; - * (Ixx Iyy Izz), Inertias; - * (q0 q1 q2 q3), attitude unit quaternion. - -OUTPUTS: - - allSteps, struct, which contains all the parameters needed. - -CALLED FUNCTIONS: - - -REVISIONS: -- #0 16/11/2018, Release, Adriano Filippo Inno -Copyright © 2021, Skyward Experimental Rocketry, AFD department -All rights reserved - -SPDX-License-Identifier: GPL-3.0-or-later - -%} -[~,firstStep] = fun(T(1),Y(:,1),varargin{:}); - -namesFields = fieldnames(firstStep); -NT = length(T); - -steps = cell(NT,1); -for k = 1:NT - [~,steps{k}] = fun(T(k),Y(:,k),varargin{:}); -end - -for i = 1:numel(namesFields) - currentFieldName = namesFields{i}; - currentField = firstStep.(currentFieldName); - - if isstruct(currentField) - namesArrays = fieldnames(currentField); - - for j = 1:numel(namesArrays) - currentArrayName = namesArrays{j}; - currentArray = currentField.(currentArrayName); - sizeArray = size(currentArray); - sizeArrayCut = sizeArray; - - if any(sizeArray == 1) - sizeArrayCut = max(sizeArray); - end - - allSteps.(currentFieldName).(currentArrayName) = zeros([sizeArrayCut,NT]); - - for k = 1:NT - currentStep = steps{k}; - - if all(sizeArray ~= 1) - allSteps.(currentFieldName).(currentArrayName)(:,:,k) =... - currentStep.(currentFieldName).(currentArrayName); - else - allSteps.(currentFieldName).(currentArrayName)(:,k) =... - currentStep.(currentFieldName).(currentArrayName); - end - end - - end - - else - sizeArray = size(currentField); - - sizeArrayCut = sizeArray; - if any(sizeArray == 1) - sizeArrayCut = max(sizeArray); - end - - allSteps.(currentFieldName) = zeros([sizeArrayCut,NT]); - - for k = 1:NT - currentStep = steps{k}; - - if all(sizeArray ~= 1) - allSteps.(currentFieldName)(:,:,k) =... - currentStep.(currentFieldName); - else - allSteps.(currentFieldName)(:,k) =... - currentStep.(currentFieldName); - end - end - - end - -end -end - +function out = recallOdeFcn(fun, t, Y, varargin) +% recallOdeFcn - This function allows to compute some parameters used +% inside the ODE integrations. +% +% INPUTS: +% - fun, function, ode function used; +% - T, double [n° variation, 1], integration time vector; +% - Y, double [n° variation, 16], state matrix, +% [ x y z | u v w | p q r | q0 q1 q2 q3 | m | Ixx Iyy Izz ]: +% * (x y z), NED{north, east, down} horizontal frame; +% * (u v w), body frame velocities; +% * (p q r), body frame angular rates; +% * m , total mass; +% * (Ixx Iyy Izz), Inertias; +% * (q0 q1 q2 q3), attitude unit quaternion. +% +% OUTPUTS: +% - allSteps, struct, which contains all the parameters needed. +% +% CALLED FUNCTIONS: - +% +% REVISIONS: +% - #0 16/11/2018, Release, Adriano Filippo Inno +% Copyright © 2021, Skyward Experimental Rocketry, AFD department +% All rights reserved +% +% SPDX-License-Identifier: GPL-3.0-or-later + +wrapper = varargin{end}; +dataRaw = wrapper.getData(); + +%% Adding t = 0, Fixing t(end) + +[~, data0] = fun(t(1), Y(:,1), varargin{1:end-1}); +data = [data0, dataRaw]; + +[~, data(end)] = fun(t(end), Y(:,end), varargin{1:end-1}); +out = wrapper.packageData(data); +wrapper.reset(); \ No newline at end of file diff --git a/functions/odeFunctions/ballistic.m b/functions/odeFunctions/ballistic.m index c8e56f1..c4ddb8d 100644 --- a/functions/odeFunctions/ballistic.m +++ b/functions/odeFunctions/ballistic.m @@ -1,4 +1,4 @@ -function [dY, parout] = ballistic(t, Y, rocket, environment, wind, angleRef, timeChangeRef, t0) +function [dY, parout] = ballistic(t, Y, rocket, environment, wind, angleRef, timeChangeRef, t0, wrapper) arguments t Y @@ -8,6 +8,7 @@ arguments angleRef double = [] timeChangeRef double = [] t0 double = 0 + wrapper = [] %DataWrapper = DataWrapper.empty end % ballistic - ode function of the 6DOF Rigid Rocket Model % @@ -319,31 +320,29 @@ dY = dY'; %% SAVING THE QUANTITIES FOR THE PLOTS -if nargout == 2 - %parout.integration.time = theta; - +if nargout == 2 || ~isempty(wrapper) parout.interp.mach = mach; parout.interp.alpha = alphaOut; parout.interp.beta = betaOut; parout.interp.alt = altitude; parout.interp.mass = m; - parout.interp.inertias = [Ixx, Iyy, Izz]; + parout.interp.inertias = [Ixx; Iyy; Izz]; - parout.wind.NED = [uw, vw, ww]; + parout.wind.NED = [uw; vw; ww]; parout.wind.body = windVels; parout.rotations.dcm = dcm; parout.velocities = vels; - parout.forces.aero = [fX, fY, fZ]; + parout.forces.aero = [fX; fY; fZ]; parout.forces.T = T; parout.air.rho = rho; parout.air.P = P; - parout.accelerations.body = [du, dv, dw]; - parout.accelerations.angular = [dp, dq, dr]; + parout.accelerations.body = [du; dv; dw]; + parout.accelerations.angular = [dp; dq; dr]; parout.coeff.CA = CA; parout.coeff.CYB = CYB; @@ -359,5 +358,7 @@ if nargout == 2 parout.coeff.Cnp = Cnp; parout.coeff.XCPlon = XCPlon; parout.coeff.XCPlat = XCPlat; - parout.coeff.XCPtot = XCPtot; + parout.coeff.XCPtot = XCPtot; + + if ~isempty(wrapper), wrapper.setCache(parout); end end \ No newline at end of file diff --git a/functions/odeOutputFcn/storeData.m b/functions/odeOutputFcn/storeData.m new file mode 100644 index 0000000..54d45fb --- /dev/null +++ b/functions/odeOutputFcn/storeData.m @@ -0,0 +1,10 @@ +function s = storeData(~, ~, flag, varargin) +s = 0; +switch flag + case {'init', 'done'} + case '' + wrapper = varargin{end}; + wrapper.setData(); +end +end + diff --git a/settings/odeConfig.m b/settings/odeConfig.m index 3d957cc..15f92ea 100644 --- a/settings/odeConfig.m +++ b/settings/odeConfig.m @@ -2,11 +2,13 @@ % ODE settings ode.finalTime = 2000; % [s] Final integration time -ode.optAscent = odeset('Events', @eventApogee, 'InitialStep', 1); % ODE options for ascend +ode.optAscent = odeset('Events', @eventApogee, 'InitialStep', 1, ... + 'OutputFcn', @storeData); % 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); % ODE options for the parachutes +ode.optParachute = odeset('Events', @eventParaCut, 'InitialStep', 1, .... + 'OutputFcn', @storeData); % ODE options for the parachutes ode.optDescent = odeset('Events', @eventLanding,... -'RelTol', 1e-3, 'AbsTol', 1e-3); % ODE options for ballistic descent +'RelTol', 1e-3, 'AbsTol', 1e-3, 'OutputFcn', @storeData); % ODE options for ballistic descent ode.optLaunchpad = odeset('Events', @eventPad); % ODE options for the launchpad phase % Settings for descent 6dof simulation -- GitLab