From 187eaa07c04223219d1c8840d07e338943c0184c Mon Sep 17 00:00:00 2001
From: Mauco03 <marco.gaibotti@skywarder.eu>
Date: Tue, 27 Feb 2024 19:11:37 +0100
Subject: [PATCH] [tests][benchmark] Added test folder
---
tests/Test.m | 18 ++
tests/ascent.m | 378 +++++++++++++++++++++++++++++++++++++++++
tests/speedBenchmark.m | 35 +++-
3 files changed, 430 insertions(+), 1 deletion(-)
create mode 100644 tests/Test.m
create mode 100644 tests/ascent.m
diff --git a/tests/Test.m b/tests/Test.m
new file mode 100644
index 0000000..fb2be57
--- /dev/null
+++ b/tests/Test.m
@@ -0,0 +1,18 @@
+classdef Test < handle
+ %TEST Summary of this class goes here
+ % Detailed explanation goes here
+
+ properties
+ a double
+ b double
+ c double
+ end
+
+ methods
+ function obj = Test(a, b, c)
+ obj.a = a;
+ obj.b = b;
+ obj.c = c;
+ end
+ end
+end
\ No newline at end of file
diff --git a/tests/ascent.m b/tests/ascent.m
new file mode 100644
index 0000000..77219ad
--- /dev/null
+++ b/tests/ascent.m
@@ -0,0 +1,378 @@
+function [dY, parout] = ascent(t, Y, settings)
+% 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 = 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;
+ uncert = [0, 0];
+
+ if not(settings.wind.input) && not(settings.wind.model)
+ uw = settings.constWind(1); vw = settings.constWind(2); ww = settings.constWind(3);
+ end
+end
+
+%% 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
+
+ if settings.stoch.N > 1
+ [uw, vw, ww] = windMatlabGenerator(settings, z, t, Hour, Day);
+ else
+ [uw, vw, ww] = windMatlabGenerator(settings, z, t);
+ 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 + 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 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);
+
+% 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 < 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
+
+ %-----------------------------------------------------
+ %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/speedBenchmark.m b/tests/speedBenchmark.m
index 5c2c39e..1de77fa 100644
--- a/tests/speedBenchmark.m
+++ b/tests/speedBenchmark.m
@@ -1,5 +1,35 @@
clear; clc; close all;
+%% Low overhead
+
+testClass = Test(1, 2, 3);
+testStruct = struct;
+
+testStruct.a = 1;
+testStruct.b = 2;
+testStruct.c = 3;
+
+len = 3*10000;
+arr1 = zeros(1, len);
+arr2 = zeros(1, len);
+
+fields = fieldnames(testClass);
+
+tic
+for i = 1:len
+ arr1(i) = testStruct.(fields{mod(i,3) + 1});
+end
+fprintf('readTime, struct: %f s\n', toc);
+
+tic
+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
ms = Mission(true);
cl = Motor(ms);
@@ -29,4 +59,7 @@ tic
for i = 1:len
d(i) = st.xCg(i);
end
-fprintf('elementWiseAssignment, struct: %f s\n', toc);
\ No newline at end of file
+fprintf('elementWiseAssignment, struct: %f s\n', toc);
+
+%% Integration test
+
--
GitLab