diff --git a/missions/2024_Lyra_Roccaraso_September/config/environmentConfig.m b/missions/2024_Lyra_Roccaraso_September/config/environmentConfig.m
index 1f90f25bdb1ff9f1093f834625cccd88d6c700d4..8b6efbf34b41a527e53d72037e6812f7748bb6dc 100644
--- a/missions/2024_Lyra_Roccaraso_September/config/environmentConfig.m
+++ b/missions/2024_Lyra_Roccaraso_September/config/environmentConfig.m
@@ -2,6 +2,9 @@
%
% 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.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
diff --git a/missions/2024_Lyra_Roccaraso_September/config/inertiaConfig.m b/missions/2024_Lyra_Roccaraso_September/config/inertiaConfig.m
index dbb6cf0bc4816e127ed98ad78044050e77dc0446..e601a9684e45f77cdd2d44cd36216bea48280e35 100644
--- a/missions/2024_Lyra_Roccaraso_September/config/inertiaConfig.m
+++ b/missions/2024_Lyra_Roccaraso_September/config/inertiaConfig.m
@@ -12,24 +12,4 @@
%%% 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
-
-%%% Inertias for engine only
-inertia.engine.Ixx = HREmotors(iMotor).Ixx; % [kg*m^2] Inertia to x-axis
-inertia.engine.Iyy = HREmotors(iMotor).Iyy; % [kg*m^2] Inertia to y-axis
-inertia.engine.Izz = HREmotors(iMotor).Izz; % [kg*m^2] Inertia to z-axis
-
-%%% Total inertias
-inertia.Ixx = inertia.rocketIxx + engineIxx;
-
-inertia.Iyy = inertia.rocketIyy + ((inertia.xcgNoMot - inertia.xcg).^2) .* inertia.mNoMot + ...
- (engineIyy + ((inertia.motor.xcg + inertia.LcenterRocket - inertia.xcg).^2 ) .* (inertia.motor.expM + inertia.motor.mc));
-
-inertia.Izz = inertia.rocketIzz + ((inertia.xcgNoMot - inertia.xcg).^2) .* inertia.mNoMot + ...
- (engineIzz + ((inertia.motor.xcg + inertia.LcenterRocket - inertia.xcg).^2) .* (inertia.motor.expM + inertia.motor.mc));
-
-inertia.I = [inertia.Ixx; inertia.Iyy; inertia.Izz];
-
-%%% Inertias derivatives (avoiding calculations in ascent.m)
-% inertia.Idot = diff(inertia.I')'./diff(inertia.motor.expTime);
-% inertia.Idot(:, end+1) = inertia.Idot(:, end);
\ No newline at end of file
+inertia.rocket.Izz = 12.07701314; % [kg*m^2] Inertia to z-axis
\ No newline at end of file
diff --git a/missions/2024_Lyra_Roccaraso_September/config/windConfig.m b/missions/2024_Lyra_Roccaraso_September/config/windConfig.m
index de6c16f1b0bf632e54a54e9e653288d95a4c9d60..cc0c2f183c38209eee108757c39a0c9b95a9e4a3 100644
--- a/missions/2024_Lyra_Roccaraso_September/config/windConfig.m
+++ b/missions/2024_Lyra_Roccaraso_September/config/windConfig.m
@@ -16,4 +16,6 @@ wind.magDistribution = ["g", "u", "u"];
wind.mag = [7 2 10;
0.5 9 20];
wind.azDistribution = ["u", "u", "u"];
-wind.az = 90*pi/180 * ones(2,3);
\ No newline at end of file
+wind.az = 90*pi/180 * ones(2,3);
+
+wind.constwind = [10, 10, 10];
\ No newline at end of file
diff --git a/missions/classes/Config.m b/missions/classes/Config.m
index 0ae56e2f49a65a174c8047b9c316188770de4384..0a3c52504db3535b2ce1bc4c60a54ad29ba6d7ca 100644
--- a/missions/classes/Config.m
+++ b/missions/classes/Config.m
@@ -9,11 +9,7 @@ classdef(Abstract) Config < handle
properties(Access = protected, Abstract)
configName {mustBeTextScalar}
end
-
- properties(Access = public)
- a = 0;
- end
-
+
methods(Access = protected)
function loadConfig(obj)
fileName = obj.configName;
diff --git a/missions/classes/Environment.m b/missions/classes/Environment.m
index cd0f2f2e180a71ee4ad85aea98abd770e88a2b4e..c8260d9ee2ad9098cfd5e34898d464932359ae73 100644
--- a/missions/classes/Environment.m
+++ b/missions/classes/Environment.m
@@ -3,6 +3,8 @@ classdef Environment < Config
% Detailed explanation goes here
properties
+ omega double % [rad] Minimum Elevation Angle, user input in degrees (ex. 80)
+ phi double % [rad] Minimum Azimuth Angle from North Direction, user input in degrees (ex. 90)
lat0 double % [deg] Launchpad latitude
lon0 double % [deg] Launchpad longitude
z0 double % [m] Launchpad Altitude
diff --git a/missions/classes/Inertia.asv b/missions/classes/Inertia.asv
new file mode 100644
index 0000000000000000000000000000000000000000..429f4efa6fbc37c599dcd7a875164132ec443a2a
--- /dev/null
+++ b/missions/classes/Inertia.asv
@@ -0,0 +1,49 @@
+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
new file mode 100644
index 0000000000000000000000000000000000000000..8ea43c9e27bc9af4f58f74318ac6bdc4ea910aa4
--- /dev/null
+++ b/missions/classes/Inertia.m
@@ -0,0 +1,62 @@
+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
+ motorIyy
+ motorIzz
+ Ixx
+ Iyy
+ Izz
+ Idot
+ end
+
+ properties(Access = protected)
+ configName = 'inertiaConfig.m'
+ mission Mission
+ motor Motor
+ mass Masses
+ geometry Geometries
+ end
+
+ methods
+ function obj = Inertia(mission, motor, mass, geometry)
+ arguments (Input)
+ mission Mission = Mission()
+ motor Motor = Motor()
+ mass Masses = Masses()
+ geometry Geometries = Geometries()
+ end
+ obj.mission = mission;
+ obj.motor = motor;
+ obj.mass = mass;
+ obj.geometry = geometry;
+ if nargin == 0, return; end
+ if nargin ~= 4, error('Wrong input arguments count: only specify 4 arguments'); end
+ obj.loadConfig();
+ end
+
+ function Ixx = get.Ixx(obj)
+ Ixx = obj.rocketIxx + obj.motor.Ixx;
+ end
+ function Iyy = get.Iyy(obj)
+ Iyy = obj.rocketIyy + ((obj.geometry.xCgNoMotor - obj.geometry.xCg).^2) .* obj.mass.noMotor + ...
+ (obj.motorIyy + ((obj.motor.xCg + obj.geometry.center.length - obj.geometry.xCg).^2 ) .* obj.motor.totalMass);
+ end
+ function Izz = get.Izz(obj)
+ 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 Idot = get.Idot(obj)
+ Idot = diff([obj.Ixx; obj.Iyy; obj.Izz])'./diff(obj.motor.time);
+ Idot(:, end+1) = Idot(:, end);
+ end
+ end
+end
\ No newline at end of file
diff --git a/missions/classes/Motor.asv b/missions/classes/Motor.asv
new file mode 100644
index 0000000000000000000000000000000000000000..491fe8d3653b19bfe0ba4cc57862abe49837881f
--- /dev/null
+++ b/missions/classes/Motor.asv
@@ -0,0 +1,90 @@
+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/missions/classes/Motor.m b/missions/classes/Motor.m
index 7c3b769ec2f5ef83e7132bed323175d144a8b69e..07ab9bbe6f637c8076ee8201e4f9f0eb1d2b40b3 100644
--- a/missions/classes/Motor.m
+++ b/missions/classes/Motor.m
@@ -18,6 +18,9 @@ classdef Motor < Config
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)
@@ -69,18 +72,21 @@ classdef Motor < Config
chosenMotor = motors(strcmp({motors.MotorName}, obj.name));
if isempty(chosenMotor), error(strcat('Unable to find engine: ', obj.name)); end
- obj.motorLength = chosenMotor.L;
+ 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.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;
end
end
end
\ No newline at end of file
diff --git a/missions/classes/Wind.m b/missions/classes/Wind.m
new file mode 100644
index 0000000000000000000000000000000000000000..057539e43af9e519896d9ee3501d9dde9ea93307
--- /dev/null
+++ b/missions/classes/Wind.m
@@ -0,0 +1,32 @@
+classdef Wind < Config
+ %WIND Summary of this class goes here
+ % Detailed explanation goes here
+
+ properties
+ magnitude % [m/s] Magnitude, [min, max]
+ elevationMin % [rad] Elevation, user input in degrees, [min, max <= 90]
+ altitudes % [rad] Azimuth, user input in degrees, [min, max]
+ magDistribution
+ mag
+ azDistribution
+ az
+ constwind
+ end
+
+ properties(Access = protected)
+ configName = 'windConfig.m'
+ mission Mission
+ end
+
+
+ methods
+ function obj = Wind(mission)
+ arguments (Input)
+ mission Mission = Mission()
+ end
+ obj.mission = mission;
+ if nargin == 0, return; end
+ obj.loadConfig();
+ end
+ end
+end
\ No newline at end of file
diff --git a/tests/ascent.asv b/tests/ascent.asv
new file mode 100644
index 0000000000000000000000000000000000000000..406b1ac8ce9983b48e5622a472ddbd6a2ad008fd
--- /dev/null
+++ b/tests/ascent.asv
@@ -0,0 +1,390 @@
+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 77219ad16cd8dcb26547d26d7e948c3385cca22e..68c48c4645183f108918da7dd05339a03d8e226b 100644
--- a/tests/ascent.m
+++ b/tests/ascent.m
@@ -1,4 +1,15 @@
-function [dY, parout] = ascent(t, Y, settings)
+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
%{
@@ -39,7 +50,6 @@ SPDX-License-Identifier: GPL-3.0-or-later
%}
-
% recalling the states
% x = Y(1);
% y = Y(2);
@@ -56,50 +66,53 @@ 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
-
-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 = settings.OMEGA;
+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 = settings.constWind(1); vw = settings.constWind(2); ww = settings.constWind(3);
+ uw = wind.constWind(1); vw = wind.constWind(2); ww = wind.constWind(3);
end
-end
+% end
%% INERTIAS
+I = [inertia.Ixx, inertia.Iyy, inertia.Izz];
+
if t < tb
if t < settings.timeEngineCut
- I = interpLinear(settings.motor.expTime, settings.I, t);
- Idot = interpLinear(settings.motor.expTime, settings.Idot, t);
+ I = interpLinear(motor.time, I, t);
+ Idot = interpLinear(motor.time, inertia.Idot, t);
else
- I = settings.IengineCut;
+ I = I(:, end);
Idot = zeros(3, 1);
end
else
- if settings.timeEngineCut < tb
- I = settings.IengineCut;
- else
- I = settings.I(:, end);
- end
+ % if settings.timeEngineCut < tb
+ % I = I(:, end);
+ % else
+ I = 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);
@@ -107,16 +120,15 @@ 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);
- end
+ % 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
@@ -137,7 +149,7 @@ if -z < 0 % z is directed as the gravity vector
z = 0;
end
-absoluteAltitude = -z + settings.z0;
+absoluteAltitude = -z + environment.z0;
[~, a, P, rho] = atmosphereData(absoluteAltitude, g, local);
M = V_norm/a;
@@ -147,21 +159,21 @@ M_value = M;
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);
+ 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 = settings.expMengineCut + settings.ms;
- T = 0;
+ % 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;
+ %m = settings.ms + settings.expMengineCut;
else
- m = settings.ms + settings.motor.expM(end);
+ m = motor.totalMass(end) + mass.structure - motor.structureMass;
end
T = 0;
@@ -192,7 +204,7 @@ end
%% INTERPOLATE AERODYNAMIC COEFFICIENTS:
[coeffsValues, angle0] = interpCoeffs(t, alpha, M, beta, absoluteAltitude,...
- c, settings);
+ c, settings); % Coeffs, state, engncut
% Retrieve Coefficients
@@ -240,9 +252,9 @@ end
%%
-if -z < settings.lrampa*sin(OMEGA) % No torque on the launchpad
+if -z < environment.rampLength*sin(omega) % No torque on the launchpad
- Fg = m*g*sin(OMEGA); % [N] force due to the gravity
+ 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;
diff --git a/tests/speedBenchmark.m b/tests/speedBenchmark.m
index 1de77fa1f17d6155a5584f5a5a5cf5b7d3bb8dfa..979639065b8109f2dfd2f99db8f03db6958dd85a 100644
--- a/tests/speedBenchmark.m
+++ b/tests/speedBenchmark.m
@@ -26,7 +26,6 @@ for i = 1:len
arr2(i) = testClass.(fields{mod(i,3) + 1});
end
fprintf('readTime, class: %f s\n', toc);
-
fprintf('\n=======================\n\n');
%% High overhead
@@ -48,6 +47,7 @@ fprintf('stdAssignment, class: %f s\n', toc);
tic
c = st.xCg;
fprintf('stdAssignment, struct: %f s\n', toc);
+fprintf('\n=======================\n\n');
tic
for i = 1:len
@@ -63,3 +63,13 @@ fprintf('elementWiseAssignment, struct: %f s\n', toc);
%% Integration test
+motor = Motor();
+
+settings.wind.model = true;
+settings.wind.input = false;
+settings.wind.variable = false;
+settings.timeEngineCut = inf;
+
+settings.tControl = motor.time(end);
+settings.control = [1 3];
+settings.MachControl = 0.8;