From a7d6defb05a3c1eb31b37e260fc6a22c5835d487 Mon Sep 17 00:00:00 2001
From: Mauco03 <marco.gaibotti@skywarder.eu>
Date: Sat, 4 May 2024 19:06:10 +0200
Subject: [PATCH] [refactoring-test-recall][functions] Done testing
---
classes/DataWrapper.m | 43 +--
classes/components/Environment.m | 13 +-
functions/miscellaneous/recallOdeFcn.m | 22 +-
functions/odeFunctions/ascentMultipleAB.m | 336 ----------------------
functions/odeFunctions/ballistic.m | 22 +-
functions/odeFunctions/descentBal.m | 249 ----------------
functions/odeFunctions/descentParachute.m | 58 ++--
settings/odeConfig.m | 6 +-
8 files changed, 93 insertions(+), 656 deletions(-)
delete mode 100644 functions/odeFunctions/ascentMultipleAB.m
delete mode 100644 functions/odeFunctions/descentBal.m
diff --git a/classes/DataWrapper.m b/classes/DataWrapper.m
index 04c8273..d708671 100644
--- a/classes/DataWrapper.m
+++ b/classes/DataWrapper.m
@@ -11,43 +11,40 @@ classdef DataWrapper < handle
end
properties(Access = private)
- cache
- counter int32 = 0
- dummyStrc struct
- chunkSize int32 % Size of preallocated memory
- dataLen int32 % Current size of cache
+ counter int32
+ shift int32
+ dummyStrc struct
+ chunkSize int32 % Size of preallocated memory
+ dataLen int32 % Current size of cache
end
methods
- function obj = DataWrapper(dataDummy, chunkSize)
+ function obj = DataWrapper(dataDummy, chunkSize, shift)
arguments
dataDummy struct
- chunkSize {mustBeInteger, mustBePositive} = 5
+ chunkSize {mustBeInteger, mustBePositive} = 1000
+ shift = 0
end
obj.chunkSize = chunkSize;
+ obj.shift = shift; % Useful to preallocate initial value
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;
+ obj.reset();
end
function setCache(obj, data)
- obj.cache = data;
+ obj.data(obj.counter + 1) = data;
end
function setData(obj)
obj.counter = obj.counter + 1;
- if obj.counter > obj.dataLen
+ 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)
@@ -55,17 +52,20 @@ classdef DataWrapper < handle
obj.reset();
end
- function out = getData(obj)
- out = obj.data(1:obj.counter);
+ function out = getData(obj, startIdx)
+ arguments
+ obj
+ startIdx = 1
+ end
+ out = obj.data(startIdx: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;
+ obj.counter = obj.shift;
end
end
@@ -85,6 +85,11 @@ classdef DataWrapper < handle
if isstruct(currentValue)
out.(currentField) = DataWrapper.packageData([data.(currentField)]);
else
+ sz = size(data(1).(currentField));
+ if all(sz > 1)
+ szIn = num2cell(sz); szIn{end+1} = []; %#ok<AGROW>
+ currentValue = reshape(currentValue, szIn{:});
+ end
out.(currentField) = currentValue;
end
end
diff --git a/classes/components/Environment.m b/classes/components/Environment.m
index bacae3c..cbdc0e7 100644
--- a/classes/components/Environment.m
+++ b/classes/components/Environment.m
@@ -29,8 +29,10 @@ classdef Environment < Component
properties(SetAccess = private)
g0 double % [-] Gravity costant at launch latitude and altitude
+ local double % [-] Vector conatining inputs for atmosphereData
pinDistance double % [m] Distance of the upper pin from the rail base (upper pin-boat + boat-rail base)
effectiveRampLength double % [m] Total launchpad length
+ effectiveRampAltitude double % [m] Projection of effectiveRampLength on z axis
end
properties(Access = protected)
@@ -56,7 +58,8 @@ classdef Environment < Component
obj.updateG0;
obj.updatePinDistance;
- obj.updateEffectiveRampLength;
+ obj.updateRamp;
+ obj.updateLocal;
end
end
@@ -71,8 +74,14 @@ classdef Environment < Component
+ obj.motor.tankLength;
end
- function obj = updateEffectiveRampLength(obj)
+ function obj = updateRamp(obj)
obj.effectiveRampLength = obj.rampLength - obj.pinDistance;
+ obj.effectiveRampAltitude = obj.effectiveRampLength*sin(obj.omega);
+ end
+
+ function obj = updateLocal(obj)
+ obj.local = [obj.z0, obj.temperature, ... % vector containing inputs for atmosphereData
+ obj.pressure, obj.rho];
end
end
end
\ No newline at end of file
diff --git a/functions/miscellaneous/recallOdeFcn.m b/functions/miscellaneous/recallOdeFcn.m
index 6af4a4e..68ad613 100644
--- a/functions/miscellaneous/recallOdeFcn.m
+++ b/functions/miscellaneous/recallOdeFcn.m
@@ -1,4 +1,8 @@
-function out = recallOdeFcn(fun, t, Y, varargin)
+function out = recallOdeFcn(solution, startIdx)
+arguments
+ solution
+ startIdx = 1
+end
% recallOdeFcn - This function allows to compute some parameters used
% inside the ODE integrations.
%
@@ -26,14 +30,14 @@ function out = recallOdeFcn(fun, t, Y, varargin)
%
% SPDX-License-Identifier: GPL-3.0-or-later
-wrapper = varargin{end};
-dataRaw = wrapper.getData();
+args = solution.extdata.varargin;
+fun = solution.extdata.odefun;
-%% Adding t = 0, Fixing t(end)
+wrapper = args{end};
+data = wrapper.getData(startIdx);
-[~, data0] = fun(t(1), Y(:,1), varargin{1:end-1});
-data = [data0, dataRaw];
+%% Fixing points at t = 0, t = end
-[~, data(end)] = fun(t(end), Y(:,end), varargin{1:end-1});
-out = wrapper.packageData(data);
-wrapper.reset();
\ No newline at end of file
+[~, data(1)] = fun(solution.x(startIdx), solution.y(:,startIdx), args{1:end-1});
+[~, data(end)] = fun(solution.x(end), solution.y(:,end), args{1:end-1});
+out = wrapper.packageData(data);
\ No newline at end of file
diff --git a/functions/odeFunctions/ascentMultipleAB.m b/functions/odeFunctions/ascentMultipleAB.m
deleted file mode 100644
index 052faaf..0000000
--- a/functions/odeFunctions/ascentMultipleAB.m
+++ /dev/null
@@ -1,336 +0,0 @@
-function [dY, parout] = ascentMultipleAB(t, Y, t0, iAB, settings)
-%{
-ascentMultipleAB - ode function of the 6DOF Rigid Rocket Model with
- multiple and smooth airbrakes opening configurations
-
-INPUTS:
-- t, double [1, 1] integration time [s];
-- Y, double [13, 1] state vector [ x y z | u v w | p q r | q0 q1 q2 q3]:
-
- * (x y z), NED{north, east, down} horizontal frame;
- * (u v w), body frame velocities;
- * (p q r), body frame angular rates;
- * (q0 q1 q2 q3), attitude unit quaternion;
-- t0, double [1, 1] starting time for airbrakes change of configuration;
-- iAB, double [1, 1] i-th airbrake configuration;
-- settings, struct(motor, CoeffsE, CoeffsF, para, ode, stoch, prob, wind), rocket data structure;
-
-OUTPUTS:
-- dY, double [16, 1] state derivatives
-- parout, struct, interesting fligth quantities structure (aerodyn coefficients, forces and so on..)
-
-
-CALLED FUNCTIONS: windMatlabGenerator, windInputGenerator, quatToDcm, interpCoeffs
-
-NOTE: To get the NED velocities the body-frame must be multiplied by the
-conjugated of the current attitude quaternion
-E.G. quatrotate(quatconj(Y(:,10:13)),Y(:,4:6))
-
-REVISIONS:
--#0 02/10/2021, Release, Davide Rosato
-
-Copyright © 2021, Skyward Experimental Rocketry, AFD department
-All rights reserved
-
-SPDX-License-Identifier: GPL-3.0-or-later
-
-%}
-
-
-% recalling the states
-% x = Y(1);
-% y = Y(2);
-z = Y(3);
-u = Y(4);
-v = Y(5);
-w = Y(6);
-p = Y(7);
-q = Y(8);
-r = Y(9);
-q0 = Y(10);
-q1 = Y(11);
-q2 = Y(12);
-q3 = Y(13);
-
-
-
-%% CONSTANTS
-S = settings.S; % [m^2] cross surface
-C = settings.C; % [m] caliber
-g = settings.g0/(1 + (-z*1e-3/6371))^2; % [N/kg] module of gravitational field
-tb = settings.tb; % [s] Burning Time
-local = settings.Local; % vector containing inputs for atmosphereData
-
-OMEGA = settings.OMEGA;
-%% INERTIAS
-if t < tb
- if t < settings.timeEngineCut
- I = interpLinear(settings.motor.expTime, settings.I, t);
- Idot = interpLinear(settings.motor.expTime, settings.Idot, t);
- else
- I = settings.IengineCut;
- Idot = zeros(3, 1);
- end
-else
- if settings.timeEngineCut < tb
- I = settings.IengineCut;
- else
- I = settings.I(:, end);
- end
- Idot = zeros(3, 1);
-end
-Ixx = I(1); Iyy = I(2); Izz = I(3);
-Ixxdot = Idot(1); Iyydot = Idot(2); Izzdot = Idot(3);
-
-%% QUATERION ATTITUDE
-Q = [q0 q1 q2 q3];
-Q = Q/norm(Q);
-
-%% ADDING WIND (supposed to be added in NED axes);
-if settings.wind.model
- [uw, vw, ww] = windMatlabGenerator(settings, z, t);
-else
- [uw, vw, ww] = computeWindVels(settings.wind, -z);
-end
-
-dcm = quatToDcm(Q);
-wind = dcm*[uw; vw; ww];
-
-% Relative velocities (plus wind);
-ur = u - wind(1);
-vr = v - wind(2);
-wr = w - wind(3);
-
-% Body to Inertial velocities
-Vels = dcm'*[u; v; w];
-V_norm = norm([ur vr wr]);
-
-%% ATMOSPHERE DATA
-if -z < 0 % z is directed as the gravity vector
- z = 0;
-end
-
-absoluteAltitude = -z + settings.z0;
-[~, a, P, rho] = atmosphereData(absoluteAltitude, g, local);
-
-M = V_norm/a;
-M_value = M;
-
-%% TIME-DEPENDENTS VARIABLES
-if t < tb
- if t < settings.timeEngineCut
- m = interpLinear(settings.motor.expTime, settings.mTotalTime, t);
- T = interpLinear(settings.motor.expTime, settings.motor.expThrust, t);
- Pe = interpLinear(settings.motor.expTime, settings.motor.Pe, t);
- T = T + settings.motor.Ae*(Pe - P);
- else
- m = settings.expMengineCut + settings.ms;
- T = 0;
- end
-
-else % for t >= tb the fligth condition is the empty one(no interpolation needed)
- if settings.timeEngineCut < tb
- m = settings.ms + settings.expMengineCut;
- else
- m = settings.ms + settings.motor.expM(end);
- end
- T = 0;
-end
-
-%% AERODYNAMICS ANGLES
-if not(abs(ur) < 1e-9 || V_norm < 1e-9)
- alpha = atan(wr/ur);
- beta = atan(vr/ur); % beta = asin(vr/V_norm) is the classical notation, Datcom uses this one though.
- % alpha_tot = atan(sqrt(wr^2 + vr^2)/ur); % datcom 97' definition
-else
- alpha = 0;
- beta = 0;
-end
-
-alpha_value = alpha;
-beta_value = beta;
-
-%% CHOSING THE EMPTY CONDITION VALUE
-% interpolation of the coefficients with the value in the nearest condition of the Coeffs matrix
-if t0 == 0
- c = settings.control(iAB);
-else
- c = settings.control(iAB);
- cPrev = settings.control(iAB - 1);
- delay = settings.delayControl(iAB - 1);
-end
-
-%% INTERPOLATE AERODYNAMIC COEFFICIENTS:
-[coeffsValues, angle0] = interpCoeffs(t, alpha, M, beta, absoluteAltitude, c, settings);
-
-if t0 ~= 0 && t >= t0 && t < t0 + delay
- [coeffsValuesPrev, ~] = interpCoeffs(t, alpha, M, beta, absoluteAltitude, cPrev, settings);
- coeffsValues = coeffsValuesPrev + (coeffsValues - coeffsValuesPrev) * (t - t0)/delay;
-end
-
-% Retrieve Coefficients
-CA = coeffsValues(1); CYB = coeffsValues(2); CY0 = coeffsValues(3);
-CNA = coeffsValues(4); CN0 = coeffsValues(5); Cl = coeffsValues(6);
-Clp = coeffsValues(7); Cma = coeffsValues(8); Cm0 = coeffsValues(9);
-Cmad = coeffsValues(10); Cmq = coeffsValues(11); Cnb = coeffsValues(12);
-Cn0 = coeffsValues(13); Cnr = coeffsValues(14); Cnp = coeffsValues(15);
-% XCP_value = coeffsValues(16);
-
-
-% compute CN,CY,Cm,Cn (linearized with respect to alpha and beta):
-alpha0 = angle0(1); beta0 = angle0(2);
-
-CN = (CN0 + CNA*(alpha - alpha0));
-CY = (CY0 + CYB*(beta - beta0));
-Cm = (Cm0 + Cma*(alpha - alpha0));
-Cn = (Cn0 + Cnb*(beta - beta0));
-
-% XCPlon = Cm/CN;
-
-if abs(alpha) <= pi/180
- XCPlon = Cma/CNA;
-else
- XCPlon = Cm/CN;
-end
-
-% XCPlat = Cn/CY;
-
-
-if abs(beta) <= pi/180
- XCPlat = Cnb/CYB;
-else
- XCPlat = Cn/CY;
-end
-
-% if Cn == 0 && CY == 0
-% XCPlat = -5;
-% end
-
-%%
-if -z < settings.lrampa*sin(OMEGA) % No torque on the launchpad
-
- Fg = m*g*sin(OMEGA); % [N] force due to the gravity
- X = 0.5*rho*V_norm^2*S*CA;
- F = -Fg +T -X;
- du = F/m;
-
- dv = 0;
- dw = 0;
- dp = 0;
- dq = 0;
- dr = 0;
-
- alpha_value = NaN;
- beta_value = NaN;
- Y = 0;
- Z = 0;
- XCPlon = NaN;
- XCPlat = NaN;
-
- if T < Fg % No velocity untill T = Fg
- du = 0;
- end
-
- XCPtot = NaN;
-else
- %% FORCES
- % first computed in the body-frame reference system
- qdyn = 0.5*rho*V_norm^2; % [Pa] dynamics pressure
- qdynL_V = 0.5*rho*V_norm*S*C;
-
- X = qdyn*S*CA; % [N] x-body component of the aerodynamics force
- Y = qdyn*S*CY; % [N] y-body component of the aerodynamics force
- Z = qdyn*S*CN; % [N] z-body component of the aerodynamics force
- Fg = dcm*[0; 0; m*g]; % [N] force due to the gravity in body frame
-
- F = Fg + [-X+T, Y, -Z]'; % [N] total forces vector
-
- %% STATE DERIVATIVES
- % velocity
- du = F(1)/m - q*w + r*v;
- dv = F(2)/m - r*u + p*w;
- dw = F(3)/m - p*v + q*u;
-
- % Rotation
- dp = (Iyy - Izz)/Ixx*q*r + qdynL_V/Ixx*(V_norm*Cl+Clp*p*C/2) - Ixxdot*p/Ixx;
- dq = (Izz - Ixx)/Iyy*p*r + qdynL_V/Iyy*(V_norm*Cm + (Cmad+Cmq)*q*C/2)...
- - Iyydot*q/Iyy;
- dr = (Ixx - Iyy)/Izz*p*q + qdynL_V/Izz*(V_norm*Cn + (Cnr*r+Cnp*p)*C/2)...
- - Izzdot*r/Izz;
-
- % Compute the aerodynamici roll angle
- [~, phi] = getAlphaPhi(alpha, beta);
-
- % Aerodynamic-force coefficient in the alpha-total plane
- CFaTot = sin(phi)*CY + cos(phi)*(-CN);
- % Aerodynanic-moment coefficient in the alpha-total plane
- CMaTot = cos(phi)*Cm - sin(phi)*Cn;
-
- if CFaTot ~= 0
- XCPtot = CMaTot/CFaTot;
- else
- XCPtot = XCPlon;
- end
-
-end
-% Quaternions
-OM = [ 0 -p -q -r ;
- p 0 r -q ;
- q -r 0 p ;
- r q -p 0 ];
-
-dQQ = 1/2*OM*Q';
-
-%% FINAL DERIVATIVE STATE ASSEMBLING
-dY(1:3) = Vels;
-dY(4) = du;
-dY(5) = dv;
-dY(6) = dw;
-dY(7) = dp;
-dY(8) = dq;
-dY(9) = dr;
-dY(10:13) = dQQ;
-dY = dY';
-
-%% SAVING THE QUANTITIES FOR THE PLOTS
-parout.integration.t = t;
-
-parout.interp.M = M_value;
-parout.interp.alpha = alpha_value;
-parout.interp.beta = beta_value;
-parout.interp.alt = -z;
-parout.interp.mass = m;
-parout.interp.inertias = [Ixx, Iyy, Izz];
-
-parout.wind.NED_wind = [uw, vw, ww];
-parout.wind.body_wind = wind;
-
-parout.rotations.dcm = dcm;
-
-parout.velocities = Vels;
-
-parout.forces.AeroDyn_Forces = [X, Y, Z];
-parout.forces.T = T;
-
-parout.air.rho = rho;
-parout.air.P = P;
-
-parout.accelerations.body_acc = [du, dv, dw];
-parout.accelerations.ang_acc = [dp, dq, dr];
-
-parout.coeff.CA = CA;
-parout.coeff.CYB = CYB;
-parout.coeff.CNA = CNA;
-parout.coeff.Cl = Cl;
-parout.coeff.Clp = Clp;
-parout.coeff.Cma = Cma;
-parout.coeff.Cmad = Cmad;
-parout.coeff.Cmq = Cmq;
-parout.coeff.Cnb = Cnb;
-parout.coeff.Cnr = Cnr;
-parout.coeff.Cnp = Cnp;
-parout.coeff.XCPlon = XCPlon;
-parout.coeff.XCPlat = XCPlat;
-
-parout.coeff.XCPtot = XCPtot;
-end
\ No newline at end of file
diff --git a/functions/odeFunctions/ballistic.m b/functions/odeFunctions/ballistic.m
index c4ddb8d..9ed6be9 100644
--- a/functions/odeFunctions/ballistic.m
+++ b/functions/odeFunctions/ballistic.m
@@ -61,16 +61,16 @@ q2 = Y(12);
q3 = Y(13);
angle = Y(14);
+dY = zeros(14, 1);
% if altitude < 0, altitude = 0; end
-%% CONSTANTS
+%% GETTING DATA
S = rocket.crossSection; % [m^2] cross surface
C = rocket.diameter; % [m] caliber
g = environment.g0/(1 + (altitude*1e-3/6371))^2; % [N/kg] module of gravitational field
tb = rocket.motor.cutoffTime; % [s] Burning Time
theta = t-t0; % [s] Time shift (needed for control)
-local = [environment.z0, environment.temperature, ... % vector containing inputs for atmosphereData
- environment.pressure, environment.rho];
+local = environment.local; % vector containing inputs for atmosphereData
omega = environment.omega;
isControlActive = false;
@@ -81,10 +81,11 @@ Q = [q0 q1 q2 q3];
Q = Q/norm(Q);
%% ADDING WIND (supposed to be added in NED axes);
-if isa(wind, 'WindMatlab')
- [uw, vw, ww] = wind.getVels(altitude, theta);
-else
- [uw, vw, ww] = wind.getVels(altitude);
+switch class(wind)
+ case 'WindMatlab'
+ [uw, vw, ww] = wind.getVels(altitude, theta);
+ case 'WindCustom'
+ [uw, vw, ww] = wind.getVels(altitude);
end
dcm = quatToDcm(Q);
@@ -206,7 +207,7 @@ end
% end
%% RAMP / FREE FLIGHT
-if altitude < environment.effectiveRampLength*sin(omega) % No torque on the launchpad
+if altitude < environment.effectiveRampAltitude % No torque on the launchpad
Fg = m*g*sin(omega); % [N] force due to the gravity
fX = 0.5*rho*velsNorm^2*S*CA;
F = -Fg +T -fX;
@@ -316,7 +317,6 @@ dY(4:6) = [du; dv; dw];
dY(7:9) = [dp; dq; dr];
dY(10:13) = dQQ;
dY(14) = dAngle;
-dY = dY';
%% SAVING THE QUANTITIES FOR THE PLOTS
@@ -360,5 +360,7 @@ if nargout == 2 || ~isempty(wrapper)
parout.coeff.XCPlat = XCPlat;
parout.coeff.XCPtot = XCPtot;
- if ~isempty(wrapper), wrapper.setCache(parout); end
+ if ~isempty(wrapper)
+ wrapper.setCache(parout);
+ end
end
\ No newline at end of file
diff --git a/functions/odeFunctions/descentBal.m b/functions/odeFunctions/descentBal.m
deleted file mode 100644
index 3e2bb42..0000000
--- a/functions/odeFunctions/descentBal.m
+++ /dev/null
@@ -1,249 +0,0 @@
-function [dY, parout] = descentBal(t, Y, settings)
-
-%{
-descentBal - ode function of the ballistic descent
-
-INPUTS:
-- t, double [1, 1], integration time [s];
-- Y, double [13, 1], state vector [ x y z | u v w | p q r | q0 q1 q2 q3]:
-
- * (x y z), NED{north, east, down} horizontal frame;
- * (u v w), body frame velocities;
- * (p q r), body frame angular rates;
- * (q0 q1 q2 q3), attitude unit quaternion;
-- settings, struct(motor, CoeffsE, CoeffsF, para, ode, stoch, prob, wind), rocket data structure;
-
-OUTPUTS:
-- dY , double [13, 1] state derivatives
-- parout, struct, interesting fligth quantities structure (aerodyn coefficients, forces and so on..)
-
-
-CALLED FUNCTIONS: windMatlabGenerator, windInputGenerator, quatToDcm, interpCoeffs
-
-NOTE: To get the NED velocities the body-frame must be multiplied by the
-conjugated of the current attitude quaternion
-E.G. quatrotate(quatconj(Y(:,10:13)),Y(:,4:6))
-
-REVISIONS:
--#0 31/12/2014, Release, Ruben Di Battista
-
--#1 16/04/2016, Second version, Francesco Colombi
-
--#2 13/01/2018, Third version, Adriano Filippo Inno
-
--#3 21/07/2023, update, Riccardo Cadamuro, Maria Teresa Cazzola
- Bugfix when the rocket begin descent with ur < 0:
- included larger aerodynamic matrix and linear
- coefficients interpolation;
- alpha and beta computed with atan2
-
-Copyright © 2021, Skyward Experimental Rocketry, AFD department
-All rights reserved
-
-SPDX-License-Identifier: GPL-3.0-or-later
-
-%}
-
-% recalling the state
-% x = Y(1);
-% y = Y(2);
-z = Y(3);
-u = Y(4);
-v = Y(5);
-w = Y(6);
-p = Y(7);
-q = Y(8);
-r = Y(9);
-q0 = Y(10);
-q1 = Y(11);
-q2 = Y(12);
-q3 = Y(13);
-m = settings.mEndAscent;
-Ixx = settings.IxxEndAscent;
-Iyy = settings.IyyEndAscent;
-Izz = settings.IzzEndAscent;
-
-persistent flippedFlag;
-
-if isempty(flippedFlag)
- flippedFlag = false;
-end
-
-%% CONSTANTS
-% Everything related to empty condition (descent-fase)
-S = settings.S; % [m^2] cross surface
-C = settings.C; % [m] caliber
-g = settings.g0/(1 + (-z*1e-3/6371))^2; % [N/kg] module of gravitational field
-T = 0;
-local = settings.Local; % vector containing inputs for atmosphereData
-
-%% QUATERION ATTITUDE
-Q = [q0 q1 q2 q3];
-Q = Q/norm(Q);
-
-%% ADDING WIND (supposed to be added in NED axes);
-if settings.wind.model
- [uw, vw, ww] = windMatlabGenerator(settings, z, t);
-else
- [uw, vw, ww] = computeWindVels(settings.wind, -z);
-end
-
-dcm = quatToDcm(Q);
-wind = dcm*[uw; vw; ww];
-
-% Relative velocities (plus wind);
-ur = u - wind(1);
-vr = v - wind(2);
-wr = w - wind(3);
-
-% Body to Inertial velocities
-Vels = dcm'*[u; v; w];
-V_norm = norm([ur vr wr]);
-
-%% ATMOSPHERE DATA
-if -z < 0 % z is directed as the gravity vector
- z = 0;
-end
-
-absoluteAltitude = -z + settings.z0;
-[~, a, P, rho] = atmosphereData(absoluteAltitude, g, local);
-
-M = V_norm/a;
-M_value = M;
-
-%% AERODYNAMICS ANGLES
-if not(abs(ur) < 1e-9 || V_norm < 1e-9)
- alpha = atan2(wr,ur);
- beta = atan2(vr,ur); % beta = asin(vr/V_norm); is the classical notation, Datcom uses this one though.
-else
- alpha = 0;
- beta = 0;
-end
-
-alpha_value = alpha;
-beta_value = beta;
-
-%% CHOSING THE CONDITION VALUE
-% interpolation of the coefficients with the value in the nearest condition of the Coeffs matrix
-
-c = 1; % descent with no aerobrakes
-
-%% INTERPOLATE AERODYNAMIC COEFFICIENTS:
-
-if abs(alpha)>25*pi/180 || abs(beta)>25*pi/180
- coeffsValues = interpN( settings.highAOA.Coeffs,...
- {settings.highAOA.State.Alphas, settings.highAOA.State.Machs, settings.highAOA.State.Betas, settings.highAOA.State.Altitudes}, ...
- [alpha, M, beta, absoluteAltitude]);
- angle0 = [alpha beta];
- flippedFlag = true;
-else
- if flippedFlag
- coeffsValues = interpN( settings.Coeffs(:, :, :, :, :, c, end),...
- {settings.State.Alphas, settings.State.Machs, settings.State.Betas, settings.State.Altitudes}, ...
- [alpha, M, beta, absoluteAltitude]);
- angle0 = [alpha beta];
- else
- [coeffsValues, angle0] = interpCoeffs(t, alpha, M, beta, absoluteAltitude, c, settings);
- end
-end
-
-% Retrieve Coefficients
-CA = coeffsValues(1); CYB = coeffsValues(2); CY0 = coeffsValues(3);
-CNA = coeffsValues(4); CN0 = coeffsValues(5); Cl = coeffsValues(6);
-Clp = coeffsValues(7); Cma = coeffsValues(8); Cm0 = coeffsValues(9);
-Cmad = coeffsValues(10); Cmq = coeffsValues(11); Cnb = coeffsValues(12);
-Cn0 = coeffsValues(13); Cnr = coeffsValues(14); Cnp = coeffsValues(15);
-
-% compute CN,CY,Cm,Cn (linearized with respect to alpha and beta):
-alpha0 = angle0(1); beta0 = angle0(2);
-
-CN = (CN0 + CNA*(alpha - alpha0));
-CY = (CY0 + CYB*(beta - beta0));
-Cm = (Cm0 + Cma*(alpha - alpha0));
-Cn = (Cn0 + Cnb*(beta - beta0));
-
-%% FORCES
-% first computed in the body-frame reference system
-qdyn = 0.5*rho*V_norm^2; % [Pa] dynamics pressure
-qdynL_V = 0.5*rho*V_norm*S*C; %
-
-fX = qdyn*S*CA; % [N] x-body component of the aerodynamics force
-fY = qdyn*S*CY; % [N] y-body component of the aerodynamics force
-fZ = qdyn*S*CN; % [N] z-body component of the aerodynamics force
-Fg = dcm*[0; 0; m*g]; % [N] force due to the gravity
-
-F = Fg +[-fX, +fY, -fZ]'; % [N] total forces vector
-
-%% STATE DERIVATIVES
-% velocity
-du = F(1)/m - q*w + r*v;
-dv = F(2)/m - r*u + p*w;
-dw = F(3)/m - p*v + q*u;
-
-% Rotation
-dp = (Iyy - Izz)/Ixx*q*r + qdynL_V/Ixx*(V_norm*Cl + Clp*p*C/2);
-dq = (Izz - Ixx)/Iyy*p*r + qdynL_V/Iyy*(V_norm*Cm + (Cmad + Cmq)*q*C/2);
-dr = (Ixx - Iyy)/Izz*p*q + qdynL_V/Izz*(V_norm*Cn + (Cnr*r + Cnp*p)*C/2);
-
-% Quaternion
-OM = 1/2* [ 0 -p -q -r ;
- p 0 r -q ;
- q -r 0 p ;
- r q -p 0 ];
-
-dQQ = OM*Q';
-
-%% FINAL DERIVATIVE STATE ASSEMBLING
-dY(1:3) = Vels;
-dY(4) = du;
-dY(5) = dv;
-dY(6) = dw;
-dY(7) = dp;
-dY(8) = dq;
-dY(9) = dr;
-dY(10:13) = dQQ;
-dY = dY';
-
-%% SAVING THE QUANTITIES FOR THE PLOTS
-parout.integration.t = t;
-
-parout.interp.M = M_value;
-parout.interp.alpha = alpha_value;
-parout.interp.beta = beta_value;
-parout.interp.alt = -z;
-
-parout.wind.NED_wind = [uw, vw, ww];
-parout.wind.body_wind = wind;
-
-parout.rotations.dcm = dcm;
-
-parout.velocities = Vels;
-
-parout.forces.AeroDyn_Forces = [fX, fY, fZ];
-parout.forces.T = T;
-
-parout.air.rho = rho;
-parout.air.P = P;
-
-parout.accelerations.body_acc = [du, dv, dw];
-parout.accelerations.ang_acc = [dp, dq, dr];
-
-parout.forces.AeroDyn_Forces = [fX, fY, fZ];
-parout.forces.T = T;
-parout.coeff.CA = CA;
-parout.coeff.CYB = CYB;
-parout.coeff.CNA = CNA;
-parout.coeff.Cl = Cl;
-parout.coeff.Clp = Clp;
-parout.coeff.Cma = Cma;
-parout.coeff.Cmad = Cmad;
-parout.coeff.Cmq = Cmq;
-parout.coeff.Cnb = Cnb;
-parout.coeff.Cnr = Cnr;
-parout.coeff.Cnp = Cnp;
-
-if eventLanding(t, Y, settings)<0
- clear flippedFlag;
-end
-
-end
\ No newline at end of file
diff --git a/functions/odeFunctions/descentParachute.m b/functions/odeFunctions/descentParachute.m
index 8d5b06f..9f1fcdd 100644
--- a/functions/odeFunctions/descentParachute.m
+++ b/functions/odeFunctions/descentParachute.m
@@ -1,12 +1,13 @@
-function [dY, parout] = descentParachute(t, Y, rocket, environment, wind, settings, descentData)
+function [dY, parout] = descentParachute(~, Y, rocket, environment, wind, settings, descentData, wrapper)
arguments
- t
+ ~
Y
rocket Rocket
environment Environment
- wind WindCustom {mustBeA(wind, {'WindCustom', 'WindMatlab'})}
+ wind % WindCustom {mustBeA(wind, {'WindCustom', 'WindMatlab'})}
settings Settings
descentData struct
+ wrapper = [] %DataWrapper = DataWrapper.empty
end
%{
descentParachute - ode function of the descent with parachute
@@ -44,7 +45,9 @@ u = Y(4);
v = Y(5);
w = Y(6);
-%% CONSTANTS
+dY = zeros(6, 1);
+
+%% GETTING DATA
if isfield(settings.simulator, 'stochNumber')
para = settings.stoch.para; %!!!!!!!!!!
@@ -64,10 +67,11 @@ local = [environment.z0, environment.temperature, ... % vector containing inpu
environment.pressure, environment.rho];
%% ADDING WIND (supposed to be added in NED axes);
-if isa(wind, 'WindMatlab')
- [uw, vw, ww] = wind.getVels(altitude, t);
-else
- [uw, vw, ww] = wind.getVels(altitude);
+switch class(wind)
+ case 'WindMatlab'
+ [uw, vw, ww] = wind.getVels(altitude, theta);
+ case 'WindCustom'
+ [uw, vw, ww] = wind.getVels(altitude);
end
windVels = [uw vw ww];
@@ -124,27 +128,27 @@ dv = F(2)/m;
dw = F(3)/m;
%% FINAL DERIVATIVE STATE ASSEMBLING
-dY(1:3) = [u v w];
+dY(1:3) = [u; v; w];
dY(4) = du;
dY(5) = dv;
dY(6) = dw;
-dY = dY';
-
%% SAVING THE QUANTITIES FOR THE PLOTS
-%if settings.plots
-parout.integration.t = t;
-parout.interp.alt = altitude;
-parout.interp.mass = m;
-parout.wind.body_wind = [uw, vw, ww];
-parout.wind.NED_wind = [uw, vw, ww];
-
-parout.air.rho = rho;
-parout.air.P = P;
-
-parout.accelerations.body_acc = [du, dv, dw];
-
-parout.velocities = [u, v, w];
-parout.accelerometer = [];
-
-%end
+if nargout == 2 || ~isempty(wrapper)
+ parout.interp.alt = altitude;
+ parout.interp.mass = m;
+ parout.wind.body = [uw, vw, ww];
+ parout.wind.NED = [uw, vw, ww];
+
+ parout.air.rho = rho;
+ parout.air.P = P;
+
+ parout.accelerations.body = [du, dv, dw];
+
+ parout.velocities = [u, v, w];
+ parout.accelerometer = [];
+
+ if ~isempty(wrapper)
+ wrapper.setCache(parout);
+ end
+end
\ No newline at end of file
diff --git a/settings/odeConfig.m b/settings/odeConfig.m
index 15f92ea..031e351 100644
--- a/settings/odeConfig.m
+++ b/settings/odeConfig.m
@@ -2,11 +2,9 @@
% ODE settings
ode.finalTime = 2000; % [s] Final integration time
-ode.optAscent = odeset('Events', @eventApogee, 'InitialStep', 1, ...
- 'OutputFcn', @storeData); % ODE options for ascend
+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, ....
- 'OutputFcn', @storeData); % ODE options for the parachutes
+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
ode.optLaunchpad = odeset('Events', @eventPad); % ODE options for the launchpad phase
--
GitLab