From 531cb7c386ee21c973e56627ee798dc6debd57e2 Mon Sep 17 00:00:00 2001
From: Mauco03 <marco.gaibotti@skywarder.eu>
Date: Thu, 29 Feb 2024 11:34:16 +0100
Subject: [PATCH] [code-refactoring][tests] WIP
---
functions/eventFunctions/eventApogee.m | 3 +-
functions/miscellaneous/interpCoeffs.m | 28 +-
.../odeFunctions/{ascent.m => OLDascent.m} | 0
.../config/environmentConfig.m | 13 +-
.../config/inertiaConfig.m | 9 +-
missions/classes/Aerodynamics.m | 44 +-
missions/classes/Environment.m | 2 +-
missions/classes/Inertia.asv | 49 ---
missions/classes/Inertia.m | 13 +-
missions/classes/Motor.asv | 90 ----
tests/Test.m | 14 +-
tests/TestOverride.m | 22 +
tests/argValidation.m | 7 +-
tests/ascent.asv | 390 ------------------
tests/ascent.m | 46 ++-
tests/matlab.mat | 3 +
tests/speedBenchmark.m | 48 ++-
tests/testAccess.m | 10 +
18 files changed, 173 insertions(+), 618 deletions(-)
rename functions/odeFunctions/{ascent.m => OLDascent.m} (100%)
delete mode 100644 missions/classes/Inertia.asv
delete mode 100644 missions/classes/Motor.asv
create mode 100644 tests/TestOverride.m
delete mode 100644 tests/ascent.asv
create mode 100644 tests/matlab.mat
create mode 100644 tests/testAccess.m
diff --git a/functions/eventFunctions/eventApogee.m b/functions/eventFunctions/eventApogee.m
index 27c3ae3..b18fcec 100644
--- a/functions/eventFunctions/eventApogee.m
+++ b/functions/eventFunctions/eventApogee.m
@@ -1,4 +1,4 @@
-function [value, isterminal, direction] = eventApogee(t, Y, settings, varargin)
+function [value, isterminal, direction] = eventApogee(t, Y, varargin)
%{
eventApogee - Event function to stop simulation at apogee checking when a
value tends to zero; the value taken is to account is the
@@ -30,6 +30,7 @@ SPDX-License-Identifier: GPL-3.0-or-later
%}
+settings = varargin{end};
Q = Y(10:13)';
% Inertial Frame velocities
diff --git a/functions/miscellaneous/interpCoeffs.m b/functions/miscellaneous/interpCoeffs.m
index d6d867a..f44568f 100644
--- a/functions/miscellaneous/interpCoeffs.m
+++ b/functions/miscellaneous/interpCoeffs.m
@@ -1,4 +1,4 @@
-function [coeffsValues, angle0] = interpCoeffs(t, alpha, M, beta, alt, c, aerodynamics, settings)
+function [coeffsValues, angle0] = interpCoeffs(t, alpha, M, beta, alt, c, aero, settings)
%{
interpCoeffs - interpolation of aerodynamic coefficients.
@@ -31,28 +31,28 @@ SPDX-License-Identifier: GPL-3.0-or-later
%}
-[nearAlpha, indexAlpha] = nearestValSorted(settings.State.Alphas*pi/180, alpha);
-[~, indexMach] = nearestValSorted(settings.State.Machs, M);
-[nearBeta, indexBeta] = nearestValSorted(settings.State.Betas*pi/180, beta);
-[~, indexAlt] = nearestValSorted(settings.State.Altitudes, alt);
+[nearAlpha, indexAlpha] = nearestValSorted(aero.state.Alphas*pi/180, alpha);
+[~, indexMach] = nearestValSorted(aero.state.Machs, M);
+[nearBeta, indexBeta] = nearestValSorted(aero.state.Betas*pi/180, beta);
+[~, indexAlt] = nearestValSorted(aero.state.Altitudes, alt);
-if t > settings.timeEngineCut && t < settings.tb
+if t > settings.timeEngineCut && t < settings.tControl
t = settings.timeEngineCut;
end
-if t >= settings.tb
- coeffsValues = settings.Coeffs(:, indexAlpha, indexMach, indexBeta, indexAlt, c, end);
+if t >= settings.tControl
+ coeffsValues = aero.coeffs(:, indexAlpha, indexMach, indexBeta, indexAlt, c, end);
else
- num = 1:1:length(settings.State.xcgTime)-1;
- indexXcg = t >= settings.State.xcgTime(1:end-1) & t < settings.State.xcgTime(2:end);
+ num = 1:1:length(aero.state.xcgTime)-1;
+ indexXcg = t >= aero.state.xcgTime(1:end-1) & t < aero.state.xcgTime(2:end);
index1 = num(indexXcg);
index2 = index1 + 1;
- VF = settings.Coeffs(:, indexAlpha, indexMach, indexBeta, indexAlt, c, index1);
- VE = settings.Coeffs(:, indexAlpha, indexMach, indexBeta, indexAlt, c, index2);
+ VF = aero.coeffs(:, indexAlpha, indexMach, indexBeta, indexAlt, c, index1);
+ VE = aero.coeffs(:, indexAlpha, indexMach, indexBeta, indexAlt, c, index2);
- deltaT = settings.State.xcgTime(index2) - settings.State.xcgTime(index1);
- coeffsValues = ( (t - settings.State.xcgTime(index1))/deltaT)*(VE - VF) + VF;
+ deltaT = aero.state.xcgTime(index2) - aero.state.xcgTime(index1);
+ coeffsValues = ( (t - aero.state.xcgTime(index1))/deltaT)*(VE - VF) + VF;
end
angle0 = [nearAlpha; nearBeta];
\ No newline at end of file
diff --git a/functions/odeFunctions/ascent.m b/functions/odeFunctions/OLDascent.m
similarity index 100%
rename from functions/odeFunctions/ascent.m
rename to functions/odeFunctions/OLDascent.m
diff --git a/missions/2024_Lyra_Roccaraso_September/config/environmentConfig.m b/missions/2024_Lyra_Roccaraso_September/config/environmentConfig.m
index 8b6efbf..37a0f6d 100644
--- a/missions/2024_Lyra_Roccaraso_September/config/environmentConfig.m
+++ b/missions/2024_Lyra_Roccaraso_September/config/environmentConfig.m
@@ -2,18 +2,19 @@
%
% Use this file to write launch data, independent from the rocket itself
-Environment.omega = 85*pi/180; % [rad] Minimum Elevation Angle, user input in degrees (ex. 80)
+environment = Environment();
+
+environment.omega = 85*pi/180; % [rad] Minimum Elevation Angle, user input in degrees (ex. 80)
environment.phi = 0*pi/180; % [rad] Minimum Azimuth Angle from North Direction, user input in degrees (ex. 90)
environment.lat0 = 41.8084579; % [deg] Launchpad latitude
environment.lon0 = 14.0546408; % [deg] Launchpad longitude
environment.z0 = 1414; % [m] Launchpad Altitude
-environment.g0 = gravitywgs84(environment.z0, environment.lat0); % Gravity costant at launch latitude and altitude
environment.pin1Length = 000; % [m] Distance from the upper pin to the upper tank cap
environment.pin2Length = 000; % [m] Distance from the lower pin to the lower tank cap
-environment.rampLength = 000; % [m] Total launchpad length
+environment.rampLength = 10; % [m] Total launchpad length
% environment.altitude ?
-environment.temperature = NaN;
-environment.pressure = NaN;
-environment.rho = NaN;
\ No newline at end of file
+environment.temperature = 0;
+environment.pressure = 0;
+environment.rho = 0;
\ No newline at end of file
diff --git a/missions/2024_Lyra_Roccaraso_September/config/inertiaConfig.m b/missions/2024_Lyra_Roccaraso_September/config/inertiaConfig.m
index e601a96..d2bbe70 100644
--- a/missions/2024_Lyra_Roccaraso_September/config/inertiaConfig.m
+++ b/missions/2024_Lyra_Roccaraso_September/config/inertiaConfig.m
@@ -6,10 +6,9 @@
% y-axis: right wing
% z-axis: downward
-[geometry, mass, motor] = ...
- loadConfig('geometryConfig.m', 'massConfig.m', 'motorConfig.m');
+inertia = Inertia();
%%% Inertias for fuselage only (no engine)
-inertia.rocket.Ixx = 0.06535397; % [kg*m^2] Inertia to x-axis
-inertia.rocket.Iyy = 12.07664659; % [kg*m^2] Inertia to y-axis
-inertia.rocket.Izz = 12.07701314; % [kg*m^2] Inertia to z-axis
\ No newline at end of file
+inertia.rocketIxx = 0.06535397; % [kg*m^2] Inertia to x-axis
+inertia.rocketIyy = 12.07664659; % [kg*m^2] Inertia to y-axis
+inertia.rocketIzz = 12.07701314; % [kg*m^2] Inertia to z-axis
\ No newline at end of file
diff --git a/missions/classes/Aerodynamics.m b/missions/classes/Aerodynamics.m
index ea36f02..5b574a1 100644
--- a/missions/classes/Aerodynamics.m
+++ b/missions/classes/Aerodynamics.m
@@ -3,7 +3,9 @@ classdef Aerodynamics < Config
% Detailed explanation goes here
properties
-
+ coeffs
+ geometry
+ state
end
properties(Access = protected)
@@ -12,41 +14,23 @@ classdef Aerodynamics < Config
end
methods
- function obj = Aerodynamics(inputArg1,inputArg2)
- %AERODYNAMICS Construct an instance of this class
- % Detailed explanation goes here
- obj.Property1 = inputArg1 + inputArg2;
- end
-
- function outputArg = method1(obj,inputArg)
- %METHOD1 Summary of this method goes here
- % Detailed explanation goes here
- outputArg = obj.Property1 + inputArg;
+ function obj = Aerodynamics(mission)
+ arguments (Input)
+ mission Mission = Mission()
+ end
+ obj.mission = mission;
+ obj.loadData();
end
end
methods(Access = private)
function obj = loadData(obj)
- if isempty(obj.mission.name) || isempty(obj.name), return; end
- load(fullfile(obj.mission.dataPath, 'aeroCoefficients.mat'), 'CoeffsTot', 'Geometry', 'State');
- chosenMotor = motors(strcmp({motors.MotorName}, obj.name));
- if isempty(chosenMotor), error(strcat('Unable to find engine: ', obj.name)); end
+ if isempty(obj.mission.name), return; end
+ load(fullfile(obj.mission.dataPath, 'aerodynamics.mat'), 'CoeffsTot', 'Geometry', 'State');
- obj.motorLength = chosenMotor.L;
- obj.tankLength = chosenMotor.Ltank;
- obj.time = chosenMotor.t;
- obj.thrust = chosenMotor.T;
- obj.fuelMass = chosenMotor.mFu;
- obj.oxidizerMass = chosenMotor.mOx;
- obj.propellantMass = chosenMotor.m;
- obj.structureMass = chosenMotor.mc;
- obj.xCg = chosenMotor.xcg;
- obj.pe = chosenMotor.Pe;
- obj.ae = chosenMotor.Ae;
- obj.fuselageMass = chosenMotor.mFus;
- obj.Ixx = chosenMotor.Ixx;
- obj.Iyy = chosenMotor.Iyy;
- obj.Izz = chosenMotor.Izz;
+ obj.coeffs = CoeffsTot;
+ obj.geometry = Geometry;
+ obj.state = State;
end
end
end
diff --git a/missions/classes/Environment.m b/missions/classes/Environment.m
index c8260d9..b02e6ea 100644
--- a/missions/classes/Environment.m
+++ b/missions/classes/Environment.m
@@ -32,7 +32,7 @@ classdef Environment < Config
function obj = Environment(mission, motor)
arguments
mission Mission = Mission()
- motor Motor = Mission()
+ motor Motor = Motor()
end
obj.mission = mission;
obj.motor = motor;
diff --git a/missions/classes/Inertia.asv b/missions/classes/Inertia.asv
deleted file mode 100644
index 429f4ef..0000000
--- a/missions/classes/Inertia.asv
+++ /dev/null
@@ -1,49 +0,0 @@
-classdef Inertia < Config
- %INERTIA Summary of this class goes here
- % Detailed explanation goes here
-
- properties
- rocketIxx % [kg*m^2] Inertia to x-axis
- rocketIyy % [kg*m^2] Inertia to y-axis
- rocketIzz % [kg*m^2] Inertia to z-axis
- end
-
- properties(Dependent)
- motorIxx % [kg*m^2] Inertia to x-axis
- motorIyy % [kg*m^2] Inertia to y-axis
- motorIzz % [kg*m^2] Inertia to z-axis
- Ixx
- Iyy
- Izz
- end
-
- properties(Access = protected)
- configName = 'inertiaConfig.m'
- mission Mission
- motor Motor
- mass Masses
- geometry Geome
- end
-
- methods
- function obj = Inertia(mission, motor)
- arguments (Input)
- mission Mission = Mission()
- motor Motor = Motor()
- mass Masses = Masses()
- end
- obj.mission = mission;
- obj.motor = motor;
- obj.mass = mass;
- if nargin == 0, return; end
- if nargin ~= 3, error('Wrong input arguments count: only specify 2 arguments'); end
- obj.loadConfig();
- end
-
- function outputArg = method1(obj,inputArg)
- %METHOD1 Summary of this method goes here
- % Detailed explanation goes here
- outputArg = obj.Property1 + inputArg;
- end
- end
-end
\ No newline at end of file
diff --git a/missions/classes/Inertia.m b/missions/classes/Inertia.m
index 8ea43c9..f7d800e 100644
--- a/missions/classes/Inertia.m
+++ b/missions/classes/Inertia.m
@@ -54,8 +54,19 @@ classdef Inertia < Config
Izz = obj.rocketIzz + ((obj.geometry.xCgNoMotor - obj.geometry.xCg).^2) .* obj.mass.noMotor + ...
(obj.motorIzz + ((obj.motor.xCg + obj.geometry.center.length - obj.geometry.xCg).^2 ) .* obj.motor.totalMass);
end
+
+ function Ixx = get.motorIxx(obj)
+ Ixx = obj.motor.Ixx;
+ end
+ function Iyy = get.motorIyy(obj)
+ Iyy = obj.motor.Iyy;
+ end
+ function Izz = get.motorIzz(obj)
+ Izz = obj.motor.Izz;
+ end
+
function Idot = get.Idot(obj)
- Idot = diff([obj.Ixx; obj.Iyy; obj.Izz])'./diff(obj.motor.time);
+ Idot = diff([obj.Ixx; obj.Iyy; obj.Izz]')'./diff(obj.motor.time);
Idot(:, end+1) = Idot(:, end);
end
end
diff --git a/missions/classes/Motor.asv b/missions/classes/Motor.asv
deleted file mode 100644
index 491fe8d..0000000
--- a/missions/classes/Motor.asv
+++ /dev/null
@@ -1,90 +0,0 @@
-classdef Motor < Config
- %MOTOR Summary of this class goes here
- % Detailed explanation goes here
-
- properties
- name {mustBeTextScalar} = '' % [-] Motor name
- motorLength double % [m] Motor length
- tankLength double % [m] Tank length
- ignitionTime double % [s] Ignition transient duration
- cutoffTime double % [s] Cutoff transient duration
- time double % [s] Engine time vector
- thrust double % [N] Engine thrust vector
- fuelMass double % [kg] Fuel (grain) initial mass
- oxidizerMass double % [kg] Oxidizer initial mass
- propellantMass double % [Kg] Propellant Mass (in time)
- structureMass double % [kg] Engine Structural Mass
- fuselageMass double % [kg] Fuselage of the engine only
- xCg double % [m] Engine xcg from tank tip
- pe double % [Pa] Eflux pressure
- ae double % [Pa] Eflux Area
- Ixx double
- Iyy double
- Izz double
- end
-
- properties(Dependent)
- totalMass double % [kg] Total motor mass
- totalLength double % [m] Motor + Tank lenght
- fuselageXCg double % [m] xcg of the engine fuselage only from tank tip
- end
-
- properties(Access = protected)
- configName = 'motorConfig.m'
- mission Mission
- end
-
- methods
- function obj = Motor(mission)
- arguments (Input)
- mission Mission = Mission()
- end
- obj.mission = mission;
- if nargin == 0, return; end
- obj.loadConfig();
- obj.loadData();
- end
-
- function set.name(obj, name)
- obj.name = name;
- obj.loadData();
- end
-
- function totalMass = get.totalMass(obj)
- totalMass = obj.propellantMass + ...
- obj.structureMass + obj.fuselageMass;
- end
-
- function fuselageXCg = get.fuselageXCg(obj)
- fuselageXCg = (obj.motorLength - ...
- obj.tankLength)/2 + obj.tankLength;
- end
-
- function totalLength = get.totalLength(obj)
- totalLength = obj.motorLength + obj.tankLength;
- end
- end
-
- methods (Access = private)
- function obj = loadData(obj)
- if isempty(obj.mission.name) || isempty(obj.name), return; end
- load(fullfile(obj.mission.dataPath, 'motors.mat'), 'motors');
- chosenMotor = motors(strcmp({motors.MotorName}, obj.name));
- if isempty(chosenMotor), error(strcat('Unable to find engine: ', obj.name)); end
-
- obj.motorLength = chosenMotor.L;
- obj.tankLength = chosenMotor.Ltank;
- obj.time = chosenMotor.t;
- obj.thrust = chosenMotor.T;
- obj.fuelMass = chosenMotor.mFu;
- obj.oxidizerMass = chosenMotor.mOx;
- obj.propellantMass = chosenMotor.m;
- obj.structureMass = chosenMotor.mc;
- obj.xCg = chosenMotor.xcg;
- obj.pe = chosenMotor.Pe;
- obj.ae = chosenMotor.Ae;
- obj.fuselageMass = chosenMotor.mFus;
- obj.Ixx
- end
- end
-end
\ No newline at end of file
diff --git a/tests/Test.m b/tests/Test.m
index fb2be57..0c70847 100644
--- a/tests/Test.m
+++ b/tests/Test.m
@@ -1,4 +1,4 @@
-classdef Test < handle
+classdef Test < Config
%TEST Summary of this class goes here
% Detailed explanation goes here
@@ -8,11 +8,21 @@ classdef Test < handle
c double
end
+ properties(SetAccess = protected)
+ d double
+ end
+
+ properties(Access = protected)
+ configName = '';
+ mission Mission
+ end
+
methods
- function obj = Test(a, b, c)
+ function obj = Test(a, b, c, d)
obj.a = a;
obj.b = b;
obj.c = c;
+ obj.d = d;
end
end
end
\ No newline at end of file
diff --git a/tests/TestOverride.m b/tests/TestOverride.m
new file mode 100644
index 0000000..1c954df
--- /dev/null
+++ b/tests/TestOverride.m
@@ -0,0 +1,22 @@
+classdef TestOverride
+ %TESTOVERRIDE Summary of this class goes here
+ % Detailed explanation goes here
+
+ properties
+ data
+ end
+
+ methods
+ function obj = TestOverride(config)
+ arguments(Input)
+ config Config {mustBeA(config,'Config')}
+ end
+ obj.data = config;
+ end
+
+ % function obj = set.data(obj, varargin)
+ % if nargin == 2, obj.data = varargin{1}; end
+ %
+ % end
+ end
+end
\ No newline at end of file
diff --git a/tests/argValidation.m b/tests/argValidation.m
index 96633f8..b20fd1f 100644
--- a/tests/argValidation.m
+++ b/tests/argValidation.m
@@ -1,11 +1,6 @@
function test(in)
arguments(Input)
- in Config {mustBeConfig}
+ in Config {mustBeA(in, 'Config')}
end
disp(in.a);
-end
-
-function mustBeConfig(obj)
- if ~isa(obj, 'Config'), error('suca');
- end
end
\ No newline at end of file
diff --git a/tests/ascent.asv b/tests/ascent.asv
deleted file mode 100644
index 406b1ac..0000000
--- a/tests/ascent.asv
+++ /dev/null
@@ -1,390 +0,0 @@
-function [dY, parout] = ascent(t, Y, geometry, mass, motor, inertia, environment, wind, settings)
-arguments
- t
- Y
- geometry Geometries
- mass Masses
- motor Motor
- inertia Inertia
- environment Environment
- wind Wind
- settings
-end
-% geometry, motor, environment
-
-%{
-ascent - ode function of the 6DOF Rigid Rocket Model
-
-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 [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 31/12/2014, Release, Ruben Di Battista
-
--#1 16/04/2016, Second version, Francesco Colombi
-
--#2 01/01/2021, Third version, Adriano Filippo Inno
-
-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 = geometry.center.crossSection; % [m^2] cross surface
-C = geometry.center.caliber; % [m] caliber
-g = environment.g0/(1 + (-z*1e-3/6371))^2; % [N/kg] module of gravitational field
-tb = motor.time(end); % [s] Burning Time
-local = [environment.z0, ...
- environment.temperature, ...
- environment.pressure, ...
- environment.rho]; % vector containing inputs for atmosphereData
-
-% if settings.stoch.N > 1
-% OMEGA = settings.stoch.OMEGA;
-% uncert = settings.stoch.uncert;
-% Day = settings.stoch.Day;
-% Hour = settings.stoch.Hour;
-% uw = settings.stoch.uw; vw = settings.stoch.vw; ww = settings.stoch.ww;
-% else
- omega = environment.omega;
- uncert = [0, 0];
-
- if not(settings.wind.input) && not(settings.wind.model)
- uw = wind.constWind(1); vw = wind.constWind(2); ww = wind.constWind(3);
- end
-% end
-
-%% INERTIAS
-
-I = [inertia.Ixx, inertia.Iyy, inertia.Izz];
-
-if t < tb
- if t < settings.timeEngineCut
- I = interpLinear(motor.time, I, t);
- Idot = interpLinear(motor.time, inertia.Idot, t);
- else
- I = I(end);
- 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
-
- % if settings.stoch.N > 1
- % [uw, vw, ww] = windMatlabGenerator(settings, z, t, Hour, Day);
- % else
- [uw, vw, ww] = windMatlabGenerator(settings, z, t); % prende wind, env
- % end
-
-elseif settings.wind.input
- [uw, vw, ww] = windInputGenerator(settings, z, uncert);
-elseif settings.wind.variable
- [uw, vw, ww] = windVariableGenerator(z, settings.stoch);
-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 + environment.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(motor.time, motor.totalMass, t);
- T = interpLinear(motor.time, motor.thrust, t);
- Pe = interpLinear(motor.time, motor.pe, t);
- T = T + motor.ae*(Pe - P);
- else
- % m = motor.totalMass(end) + mass.structure - motor.structureMass;
- % 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 = motor.totalMass(end) + mass.structure - motor.structureMass;
- 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 t >= settings.tControl && M <= settings.MachControl
- c = settings.control;
-else
- c = 1;
-end
-
-%% INTERPOLATE AERODYNAMIC COEFFICIENTS:
-
-[coeffsValues, angle0] = interpCoeffs(t, alpha, M, beta, absoluteAltitude,...
- c, settings); % Coeffs, state, engncut
-
-% 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);
-%--------------------------------------------
-%Clb = coeffsValues(16);
-%--------------------------------------------
-
-% 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 < environment.rampLength*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
-
- %-----------------------------------------------------
- %F = Fg + [-X+T*cos(chi), Y+T*sin(chi), -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
-
-if nargout == 2
- 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.Clb = Clb;
- %--------------------------------------------
- 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/tests/ascent.m b/tests/ascent.m
index 521606d..bd65b93 100644
--- a/tests/ascent.m
+++ b/tests/ascent.m
@@ -1,15 +1,16 @@
-function [dY, parout] = ascent(t, Y, geometry, mass, motor, inertia, environment, wind, settings)
-arguments
- t
- Y
- geometry Geometries
- mass Masses
- motor Motor
- inertia Inertia
- environment Environment
- wind Wind
- settings
-end
+function [dY, parout] = ascent(t, Y, geometry, mass, motor, inertia, environment, wind, aerodynamics, settings)
+% arguments
+% t
+% Y
+% geometry Geometries
+% mass Masses
+% motor Motor
+% inertia Inertia
+% environment Environment
+% wind Wind
+% aerodynamics Aerodynamics
+% settings
+% end
% geometry, motor, environment
%{
@@ -27,7 +28,7 @@ INPUTS:
OUTPUTS:
- dY, double [16, 1] state derivatives
-- parout, struct, interesting fligth quantities structure (aerodyn coefficients, forces and so on..)
+- parout, struct, interesting flight quantities structure (aero coefficients, forces and so on..)
CALLED FUNCTIONS: windMatlabGenerator, windInputGenerator, quatToDcm, interpCoeffs
@@ -69,11 +70,11 @@ q3 = Y(13);
S = geometry.center.crossSection; % [m^2] cross surface
C = geometry.center.caliber; % [m] caliber
g = environment.g0/(1 + (-z*1e-3/6371))^2; % [N/kg] module of gravitational field
-tb = motor.time(end); % [s] Burning Time
-local = [environment.z0, ...
- environment.temperature, ...
- environment.pressure, ...
- environment.rho]; % vector containing inputs for atmosphereData
+tb = settings.tb; % [s] Burning Time
+% local = [environment.z0, ...
+% environment.temperature, ...
+% environment.pressure, ...
+% environment.rho]; % vector containing inputs for atmosphereData
% if settings.stoch.N > 1
% OMEGA = settings.stoch.OMEGA;
@@ -92,7 +93,7 @@ local = [environment.z0, ...
%% INERTIAS
-I = [inertia.Ixx, inertia.Iyy, inertia.Izz];
+I = [inertia.Ixx; inertia.Iyy; inertia.Izz];
if t < tb
if t < settings.timeEngineCut
@@ -151,7 +152,10 @@ if -z < 0 % z is directed as the gravity vector
end
absoluteAltitude = -z + environment.z0;
-[~, a, P, rho] = atmosphereData(absoluteAltitude, g, local);
+%[~, a, P, rho] = atmosphereData(absoluteAltitude, g, local);
+a = 6.5;
+P = 101325;
+rho = 1.293;
M = V_norm/a;
M_value = M;
@@ -204,7 +208,7 @@ end
%% INTERPOLATE AERODYNAMIC COEFFICIENTS:
[coeffsValues, angle0] = interpCoeffs(t, alpha, M, beta, absoluteAltitude,...
- c, settings); % Coeffs, state, engncut
+ c, aerodynamics, settings); % Coeffs, state, engncut
% Retrieve Coefficients
diff --git a/tests/matlab.mat b/tests/matlab.mat
new file mode 100644
index 0000000..58b5214
--- /dev/null
+++ b/tests/matlab.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:033fcc1ba45d8acefed5f829e5b4d23efdb91ce771739d890918e72ed677fc19
+size 179
diff --git a/tests/speedBenchmark.m b/tests/speedBenchmark.m
index 9796390..36762e2 100644
--- a/tests/speedBenchmark.m
+++ b/tests/speedBenchmark.m
@@ -62,8 +62,14 @@ end
fprintf('elementWiseAssignment, struct: %f s\n', toc);
%% Integration test
-
-motor = Motor();
+mission = Mission(1);
+motor = Motor(mission);
+wind = Wind(mission);
+aero = Aerodynamics(mission);
+mass = Masses(mission, motor);
+env = Environment(mission, motor);
+geo = Geometries(mission, motor, mass);
+inertia = Inertia(mission, motor, mass, geo);
settings.wind.model = true;
settings.wind.input = false;
@@ -73,3 +79,41 @@ settings.timeEngineCut = inf;
settings.tControl = motor.time(end);
settings.control = [1 3];
settings.MachControl = 0.8;
+settings.tb = motor.time(end);
+
+Q0 = angleToQuat(env.phi, env.omega, 180*pi/180)';
+
+%%% State
+X0 = [0 0 0]';
+V0 = [0 0 0]';
+W0 = [0 0 0]';
+
+Y0a = [X0; V0; W0; Q0];
+
+settings.ode.optionsasc1 = odeset('Events', @eventApogee, 'InitialStep', 1);
+tf = 2000;
+
+t = 0;
+Y = zeros(1,13);
+Y(end) = 1;
+
+%tic
+[Ta, Ya] = ode113(@ascent, [0, tf], Y0a, settings.ode.optionsasc1, geo, mass, motor, inertia, ...
+ env, wind, aero, settings); % till the apogee
+%classTime = toc;
+
+% motor = struct(motor);
+% wind = struct(wind);
+% aero = struct(aero);
+% mass = struct(mass);
+% env = struct(env);
+% geo = struct(geo);
+% inertia = struct(inertia);
+% clc;
+%
+% tic
+% [~, ~] = ode113(@ascent, [0, tf], Y0a, settings.ode.optionsasc1, geo, mass, motor, inertia, ...
+% env, wind, aero, settings); % till the apogee
+% structTime = toc;
+%
+% fprintf('speed ratio (class/struct): %f\n', (classTime / structTime));
\ No newline at end of file
diff --git a/tests/testAccess.m b/tests/testAccess.m
new file mode 100644
index 0000000..6a2a640
--- /dev/null
+++ b/tests/testAccess.m
@@ -0,0 +1,10 @@
+% Obtain the meta.class object
+mc = meta.class.fromName('Test');
+
+% Modify the access of a property (e.g., MyPublicProperty)
+prop = mc.PropertyList(strcmp({mc.PropertyList.Name}, 'd'));
+prop.SetAccess = 'protected'; % Change to protected
+
+% Modify the access of a method (e.g., myPrivateMethod)
+method = mc.MethodList(strcmp({mc.MethodList.Name}, 'myPrivateMethod'));
+method.Access = 'public'; % Change to public
--
GitLab