Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • skyward/matlab-dependencies/common
1 result
Select Git revision
Show changes
Commits on Source (16)
Showing
with 603 additions and 179 deletions
......@@ -51,8 +51,10 @@ classdef Rocket < Config
pitot Pitot % [-] Pitot component
parachutes cell % [-] (nParachutes, nStages) Parachutes onboard
coefficients struct % [-] Aerodynamic coefficients
coefficientsHighAOA struct % [-] Aerodynamic coefficients at high angle of attack
dynamicDerivatives (1, 1) logical % [-] True if dynamic derivatives will be loaded
coefficients Coefficient % [-] Aerodynamic coefficients
coefficientsHighAOA Coefficient % [-] Aerodynamic coefficients at high angle of attack
end
properties(Access = private)
......@@ -138,22 +140,62 @@ classdef Rocket < Config
obj.parachutes = paraLoader.parachutes;
if options.loadCoefficients
obj = obj.loadData();
if obj.dynamicDerivatives
coeffName = obj.motor.name;
else
coeffName = 'generic';
end
obj.coefficients = Coefficient( ...
fullfile(mission.dataPath, 'aeroCoefficients.mat'), ...
coeffName);
obj.coefficientsHighAOA = Coefficient( ...
fullfile(mission.dataPath, 'aeroCoefficientsHighAOA.mat'), ...
coeffName);
answer = '';
if isempty(obj.coefficients) || isempty(obj.coefficientsHighAOA)
if isempty(obj.coefficients.static) || isempty(obj.coefficientsHighAOA.static)
answer = questdlg(['Coefficient matrices not found. ' ...
'Do you want to create new matrices?']);
elseif options.checkGeometry && ~all(obj.checkGeometry())
elseif options.checkGeometry
[check, summary] = obj.checkGeometry();
if ~all(check)
disp(summary);
answer = questdlg(['Coefficient matrices differ from rocket geometry. ' ...
'Do you want to create new matrices?']);
end
end
switch answer
case 'Yes'
% Check for required tools
isMsa = exist(mission.msaPath, "file");
isDissile = exist(mission.dissilePath, "file");
if ~all([isMsa, isDissile])
msg = 'The following toolkits are required to update coefficients:';
if ~isMsa
msg = sprintf(...
'%s\n - <a href="https://git.skywarder.eu/afd/msa/msa-toolkit">msa-toolkit</a>:\n\t In %s\n', ...
msg, mission.msaPath);
end
if ~isDissile
msg = sprintf('%s\n - <a href="https://git.skywarder.eu/afd/msa/dissilematcom">dissilematcom</a>:\n\t In %s\n', ...
msg, mission.dissilePath);
end
error(msg);
end
% Start coefficient creation
parserPath = fullfile(mission.msaPath, 'autoMatricesProtub');
addpath(genpath(parserPath));
[obj.coefficients, obj.coefficientsHighAOA] = ...
mainAutoMatProtub(obj);
obj = obj.loadData();
case 'Cancel'
error('Rocket creation aborted')
otherwise
......@@ -164,7 +206,7 @@ classdef Rocket < Config
methods
[coeffsValues, angle0] = interpCoeffs(obj, t, alpha, mach, beta, alt, c)
checks = checkGeometry(obj)
[checks, summary] = checkGeometry(obj)
img = plot(obj)
end
......
function checks = checkGeometry(obj)
function [checks, summary] = checkGeometry(obj)
% checkGeometry - This methods checks if the rocket geometry
% is consistent with the geometry of the
% aerodynamic matrices
......@@ -6,6 +6,11 @@ function checks = checkGeometry(obj)
% OUTPUTS:
% - checks (n fields of geometry, 1): boolean value of the geometry checks
names = ["Initial xcg"; "Final xcg"; "Ogive type"; "Diameter"; "Nose Length";
"Nose p Mod"; "Nose c Mod"; "Center Length"; "Fins root chord";
"Fins free chord"; "Fins height"; "Fins XLE"; "Fins number";
"Boat Length"; "Boat final diameter"; "Boat type"];
xCgRocket = round([ ...
obj.xcg(1); ...
obj.xcg(end) ...
......@@ -16,10 +21,12 @@ function checks = checkGeometry(obj)
obj.coefficients.geometry.xcg(end);
], 3);
if (obj.parafoil.noseCMod & obj.parafoil.nosePMod) == 0 % KARMAN ogive case, no modified p and c coefficients
noseTypeCheck = strcmp(obj.coefficients.geometry.ogType, obj.parafoil.noseType);
geometryRocket = round([
obj.diameter;
obj.parafoil.noseLength;
0; 0; % Placeholders for nose p and c mod
obj.lengthCenter;
obj.rear.finsRootChord;
obj.rear.finsFreeChord;
......@@ -33,6 +40,7 @@ function checks = checkGeometry(obj)
geometryTest = round([
obj.coefficients.geometry.diameter;
obj.coefficients.geometry.lNose;
0; 0; % Placeholders for nose p and c mod
obj.coefficients.geometry.lCenter;
obj.coefficients.geometry.chord1;
obj.coefficients.geometry.chord2;
......@@ -42,42 +50,48 @@ function checks = checkGeometry(obj)
obj.coefficients.geometry.boatL;
obj.coefficients.geometry.boatD;
], 3);
else % MHAAK ogive case, modified p and c coefficients
geometryRocket = round([
obj.diameter;
obj.parafoil.noseLength;
obj.lengthCenter;
if (obj.parafoil.noseCMod && obj.parafoil.nosePMod)
geometryRocket([3, 4]) = round([ ...
obj.parafoil.noseCMod;
obj.parafoil.nosePMod;
obj.rear.finsRootChord;
obj.rear.finsFreeChord;
obj.rear.finsHeight;
obj.rear.finsDeltaXFreeChord;
obj.rear.nPanel;
obj.rear.boatLength;
obj.rear.boatFinalDiameter;
], 3);
obj.parafoil.nosePMod;], 4);
geometryTest = round([
obj.coefficients.geometry.diameter;
obj.coefficients.geometry.lNose;
obj.coefficients.geometry.lCenter;
obj.coefficients.geometry.cMod;
obj.coefficients.geometry.pMod;
obj.coefficients.geometry.chord1;
obj.coefficients.geometry.chord2;
obj.coefficients.geometry.height;
obj.coefficients.geometry.deltaXLE;
obj.coefficients.geometry.nPanel;
obj.coefficients.geometry.boatL;
obj.coefficients.geometry.boatD;
], 3);
geometryTest([3, 4]) = round([ ...
obj.parafoil.noseCMod;
obj.parafoil.nosePMod;], 4);
end
boatTypeCheck = strcmp(obj.coefficients.geometry.boatType, obj.rear.boatType);
checks = [
xCgRocket == xCgTest;
strcmp(obj.coefficients.geometry.ogType, obj.parafoil.noseType);
~obj.dynamicDerivatives | (xCgRocket == xCgTest);
noseTypeCheck;
geometryRocket == geometryTest;
strcmp(obj.coefficients.geometry.boatType, obj.rear.boatType)
boatTypeCheck;
];
if nargout == 2
% Summarize data and leave placeholders
ogiveType = ["CONE", "OGIVE", "POWER", "HAACK", "KARMAN", "MHAACK"];
boatType = ["CONE", "OGIVE"];
% Convert text to DATCOM standard. !TODO: Would be more robust if
% enums were used
boatRocket = find(strcmp(obj.rear.boatType, boatType), 1, 'first') - 1;
boatTest = find(strcmp(obj.coefficients.geometry.boatType, boatType), 1, 'first') - 1;
ogRocket = find(strcmp(obj.parafoil.noseType, ogiveType), 1, 'first') - 1;
ogTest = find(strcmp(obj.coefficients.geometry.ogType, ogiveType), 1, 'first') - 1;
rocketData = [xCgRocket; ogRocket; geometryRocket; boatRocket];
testData = [xCgTest; ogTest; geometryTest; boatTest];
checkNames = names(~checks);
if isempty(checkNames)
summary = table;
else
summary = table(rocketData(~checks), testData(~checks), ...
'RowNames', checkNames, 'VariableNames', ["Rocket", "Coefficients"]);
end
end
end
\ No newline at end of file
......@@ -23,6 +23,9 @@ classdef Mission < Config
configPath char
dataPath char
msaPath char
prpPath char
gncPath char
dissilePath char
end
properties(Access = private)
......@@ -98,6 +101,12 @@ classdef Mission < Config
obj.dataPath = fullfile(obj.currentPath, obj.name, 'data');
obj.msaPath = trimPath( ...
fullfile(obj.currentPath, '..', '..', '..', 'msa-toolkit'));
obj.prpPath = trimPath( ...
fullfile(obj.currentPath, '..', '..', '..', 'prp-toolkit'));
obj.gncPath = trimPath( ...
fullfile(obj.currentPath, '..', '..', '..', 'matlab-simulator'));
obj.dissilePath = trimPath( ...
fullfile(obj.currentPath, '..', '..', '..', 'dissilematcom'));
end
function updateMatlabPath(obj)
......
classdef CoeffType
enumeration
AlphaTotPhi
AlphaBeta
end
end
File moved
classdef Coefficient
% Coefficient: Manages the coefficients associated with a rocket
%
% Constructor:
% - Coefficient: Creates an instance of the Coefficient class.
% Loaded config: -
% Loaded data: aeroCoefficients.mat
% Arguments:
% - filePath: char, path to coefficient file to load
% - name: (optional) coefficient name. Used when dynamic
% derivatives are also needed
%
% Coefficients are retrieved from DATCOM in the following format:
% Coefficient matrix: double (15, nAlpha, nMach, nBeta, nAlt, nAbk, nXcg)
% Where the 15 coefficients are
% - CA Axial force
% - CYB Side force derivative wrt beta
% - CY0 Side force
% - CNA Normal force derivative wrt alpha
% - CN0 Normal force
% - Cl Rolling moment
% - Clp Rolling moment derivative wrt roll rate
% - Cma Pitching moment derivative wrt alpha
% - Cm0 Pitching moment
% - Cmad Pitching moment derivative wrt alpha dot
% - Cmq Pitching moment derivative wrt pitch rate
% - Cnb Yawing moment derivative wrt beta
% - Cn0 Yawing moment
% - Cnr Yawing moment derivative wrt yaw rate
% - Cnp Yawing moment derivative wrt roll rate
%
% Coefficients are then stored in the following format:
% static: double (6, nAlpha, nMach, nBeta, nAlt, nAbk)
% Where the 6 coefficients are:
% - [CA, CY, CN, Cl, Cm, Cn]
% dynamic: double (5, nAlpha, nMach, nBeta, nAlt, nAbk, nXcg)
% - [Clp, Cmad, Cmq, Cnr, Cnp]
%
% NOTE: When dynamic derivatives are not loaded, the matrix will be
% 6-D, as there will be no need to interpolate wrt to the xcg: moment
% transport formulas will be used instead
%
% NOTE: Coefficients in an interpolants are stored on the last
% dimension, instead of the first, to improve performance
properties(Dependent)
total double % Coefficients stored in DATCOM format
geometry struct % Reference geometry
state struct % Flight envelope
end
properties
finsCN double % Fins-only CN
end
properties(Dependent, SetAccess = private)
static double % Static coefficients: [CA, CY, CN, Cl, Cm, Cn]
dynamic double % Dynamic derivatives: [Clp, Cmad, Cmq, Cnr, Cnp]
end
properties(SetAccess = private)
type CoeffType % Wheter based on alphaTot-phi or alpha-beta
symmetry (1, 1) logical % Apply axial simmetry: valid for alphaTot-phi
isReady (1, 1) logical % Whether all coefficients are loaded and of the right size
isDynamic (1, 1) logical % Whether to load dynamic derivatives. Adds XCG dependece
end
properties(Access = private)
TOTAL % Cached variable for total
GEOMETRY % Cached variable for geometry
STATE % Cached variable for state
end
properties(Access = private)
staticInterpolant (1, 1) griddedInterpolant
dynamicInterpolant (1, 1) griddedInterpolant
end
methods
function obj = Coefficient(filePath, name)
arguments
filePath (1, :) char = ''
name (1, :) char = 'generic'
end
if isempty(filePath), return; end
obj = obj.loadData(filePath, name);
obj.isReady = obj.checkProperties();
obj = obj.updateInterpolants();
end
function [coefficients] = get(obj, alpha, mach, beta, altitude, airbakes, xcg)
% GET: Retrieve aerodynamic coefficients for specified flight conditions
%
% This method interpolates the aerodynamic coefficients based on
% the provided flight conditions and adjusts them for the given
% center of gravity (xcg).
%
% Arguments:
% - alpha double, angle of attack [rad]
% - mach double, Mach number [-]
% - beta double, sideslip angle [rad]
% - altitude double, altitude [m]
% - airbakes double, airbrake deployment [-]
% - xcg double, center of gravity position [m]
%
% Returns:
% - coefficients: double (11, 1), aerodynamic coefficients
% [CA, CY, CN, Cl, Cm, Cn, Clp, Cmad, Cmq, Cnr, Cnp]
%
% NOTE: Static coefficients are interpolated and adjusted for
% the specified xcg. Dynamic derivatives are included only if
% loadDynamic is true.
if ~obj.isReady, error('Cannot interpolate coefficients: check that all dimensions match'); end
% Interpolating coefficients
coefficients = zeros(11, 1);
if obj.type == CoeffType.AlphaBeta
coefficients(1:6) = obj.staticInterpolant(alpha, mach, beta, altitude, airbakes);
else
[alphaTot, phi] = getAlphaPhi(alpha, beta);
phi = wrapTo2Pi(phi);
finAngle = 2*pi / obj.GEOMETRY.nPanel; % Angle between two fin panels
n = floor(phi / finAngle);
deltaPhi = n*finAngle;
psi = phi - deltaPhi; % Angle wrapped to finAngle
coefficients(1:6) = obj.staticInterpolant(alphaTot, mach, psi, altitude, airbakes); % Get coeffs in limited range
if n > 1 % If necessary, perform rotation on different frame
R = [cos(deltaPhi), -sin(deltaPhi);
sin(deltaPhi), cos(deltaPhi)];
coefficients([2, 3]) = R * coefficients([2, 3]);
coefficients([5, 6]) = R' * coefficients([5, 6]);
end
end
% Transporting static coefficients to new xcg
% C_B = C_A + d / c * [0; -CN; CY]_A <- NOTE: Non torna il meno su CN
d = xcg - obj.GEOMETRY.xcg(1);
l = obj.GEOMETRY.diameter;
CY = coefficients(2); CN = coefficients(3);
forceCoeffs = [0; CN; CY];
coefficients(4:6) = coefficients(4:6) + d/l * forceCoeffs;
if ~obj.isDynamic, return; end
if obj.type == CoeffType.AlphaTotPhi
error('Alpha-Phi coefficients do not support dynamic derivatives yet')
end
coefficients(7:11) = obj.dynamicInterpolant(alpha, mach, beta, altitude, airbakes, xcg);
end
end
methods % Getters
function value = get.total(obj)
value = obj.TOTAL;
end
function value = get.static(obj)
if isempty(obj.TOTAL), value = []; return; end
value = obj.TOTAL([1, 3, 5, 6, 9, 13], :, :, :, :, :, 1);
end
function value = get.dynamic(obj)
if ~obj.isDynamic || isempty(obj.TOTAL), value = []; return; end
value = obj.TOTAL([7, 10, 11, 14, 15], :, :, :, :, :, :);
end
function value = get.geometry(obj)
value = obj.GEOMETRY;
end
function value = get.state(obj)
value = obj.STATE;
end
end
methods % Setters
function obj = set.total(obj, value)
obj.TOTAL = value;
obj.isDynamic = size(value, 7) ~= 1;
obj.isReady = obj.checkProperties();
obj = obj.updateInterpolants();
end
function obj = set.geometry(obj, value)
if ~isempty(obj.STATE) && isfield(obj.STATE, 'phis') && ~isempty(obj.STATE.phis)
obj.symmetry = (360/obj.STATE.phis(end) == value.nPanel);
end
obj.GEOMETRY = value;
obj.isReady = obj.checkProperties();
obj = obj.updateInterpolants();
end
function obj = set.state(obj, value)
if isfield(value, 'phis') && ~isempty(value.phis)
if ~isempty(value.betas)
error('Cannot set both alpha, beta and phi');
end
if ~isempty(obj.GEOMETRY)
obj.symmetry = (360/value.phis(end) == obj.GEOMETRY.nPanel);
end
obj.type = CoeffType.AlphaTotPhi;
else
obj.type = CoeffType.AlphaBeta;
end
obj.STATE = value;
obj.isReady = obj.checkProperties();
obj = obj.updateInterpolants();
end
end
methods(Access = private)
function obj = loadData(obj, filePath, name)
% Load the coefficients from the specified file
% filePath: char, path to coefficient file to load
% name: (optional) coefficient name. Used when dynamic
% derivatives are also needed
%
% NOTE: Coefficient alpha and beta derivatives are not loaded
% as coefficients are already interpolated linearly
if ~exist(filePath, 'file'), return; end
warning('off'); % Disable warning in case coefficients are not present
dataCoeffs = load(filePath, name);
warning('on');
if ~isfield(dataCoeffs, name), return; end
dataCoeffs = dataCoeffs.(name);
obj.finsCN = dataCoeffs.finsCN;
obj.GEOMETRY = dataCoeffs.geometry;
obj.STATE = dataCoeffs.state;
% Load coefficients
obj.total = dataCoeffs.total;
end
function ready = checkProperties(obj)
% Check if STATIC, DYNAMIC, GEOMETRY, and STATE are non-empty
ready = ~isempty(obj.static) && ...
~isempty(obj.GEOMETRY) && ~isempty(obj.STATE) && ...
~(isempty(obj.dynamic) && obj.isDynamic);
if ~ready, return; end
alpha = obj.STATE.alphas;
mach = obj.STATE.machs;
altitude = obj.STATE.altitudes;
airbakes = obj.STATE.hprot;
xcg = obj.GEOMETRY.xcg;
if obj.type == CoeffType.AlphaBeta
theta = obj.STATE.betas; % Converting to rad
else
theta = obj.STATE.phis; % Converting to rad
end
gridVecs = {alpha, mach, theta, altitude, airbakes, xcg};
dims = cellfun(@(x) length(x), gridVecs);
dims(dims == 0) = 1; % Empty case is the same as scalar case
ready = all(size(obj.total, 2:7) == dims);
end
function obj = updateInterpolants(obj)
% Update static and dynamic interpolants
if ~obj.isReady, return; end
%% Retrieve flight conditions
alpha = obj.STATE.alphas*pi/180; % Converting to rad
mach = obj.STATE.machs;
altitude = obj.STATE.altitudes;
xcg = obj.GEOMETRY.xcg;
if isempty(obj.STATE.hprot), airbakes = 0;
else, airbakes = obj.STATE.hprot/obj.STATE.hprot(end); % Normalizing on height
end
if obj.type == CoeffType.AlphaBeta
theta = obj.STATE.betas*pi/180; % Converting to rad
else
theta = obj.STATE.phis*pi/180; % Converting to rad
end
gridVecs = {alpha, mach, theta, altitude, airbakes, xcg};
% Find singleton dims (last dimension is coefficients and will not be of length 1)
singletonDims = cellfun(@(x) isscalar(x) || isempty(x), gridVecs);
gridVecs(singletonDims) = deal({[0, 1]});
%% Create interpolants
staticCoeffs = permute(obj.static, [2, 3, 4, 5, 6, 1]);
staticCoeffs = repmat(staticCoeffs, singletonDims(1:end-1) + 1); % Exclude xcg for static coefficients
obj.staticInterpolant = griddedInterpolant(gridVecs(1:end-1), staticCoeffs, 'linear', 'nearest');
if ~obj.isDynamic, return; end
%% Handle dynamic derivatives
dynamicCoeffs = permute(obj.dynamic, [2, 3, 4, 5, 6, 7, 1]);
% Flip xcg in case of descending order
if ~isscalar(xcg) && xcg(2) < xcg(1)
gridVecs{end} = flip(xcg);
dynamicCoeffs = flip(dynamicCoeffs, 7);
end
obj.dynamicInterpolant = griddedInterpolant(gridVecs, dynamicCoeffs, 'linear', 'nearest');
end
end
end
\ No newline at end of file
......@@ -9,6 +9,7 @@ classdef Wind < Config
% - mission: Mission, mission object
% - varIn: (optional) config source. Alternative to config.m
% file
properties(Dependent)
altitudes (1, :)
......@@ -25,6 +26,9 @@ classdef Wind < Config
end
properties(Access = private)
isReady (1, 1) logical % true if all properties are set
dataInerpolant (1, 1) griddedInterpolant % [-] Interpolates [magnitude, azimuth, elevation]
ALTITUDES
MAGNITUDE_DISTRIBUTION = ''
......@@ -42,76 +46,70 @@ classdef Wind < Config
magnitude
end
properties(Access = private)
checks (1, 3) logical % true if #[magnitude, azimuth, elevation] == #altitudes
end
properties(Constant, Access = protected)
configName = 'windConfig.m'
variableName = 'windCustom'
end
methods % Getters / Setters
function obj = set.altitudes(obj, value)
if length(value) ~= length(obj.ALTITUDES)
obj.ALTITUDES = value;
obj.magnitude = obj.updateDistribution(obj.MAGNITUDE_DISTRIBUTION, obj.MAGNITUDE_PARAMETERS);
obj.azimuth = obj.updateDistribution(obj.AZIMUTH_DISTRIBUTION, obj.AZIMUTH_PARAMETERS);
obj.elevation = obj.updateDistribution(obj.ELEVATION_DISTRIBUTION, obj.ELEVATION_PARAMETERS);
else
obj.ALTITUDES = value;
methods % Getters
function value = get.altitudes(obj), value = obj.ALTITUDES; end
function value = get.magnitudeDistribution(obj), value = obj.MAGNITUDE_DISTRIBUTION; end
function value = get.azimuthDistribution(obj), value = obj.AZIMUTH_DISTRIBUTION; end
function value = get.elevationDistribution(obj), value = obj.ELEVATION_DISTRIBUTION; end
function value = get.magnitudeParameters(obj), value = obj.MAGNITUDE_PARAMETERS; end
function value = get.azimuthParameters(obj), value = obj.AZIMUTH_PARAMETERS; end
function value = get.elevationParameters(obj), value = obj.ELEVATION_PARAMETERS; end
end
obj.checks(1) = length(obj.magnitude) == length(obj.ALTITUDES);
obj.checks(2) = length(obj.azimuth) == length(obj.ALTITUDES);
obj.checks(3) = length(obj.elevation) == length(obj.ALTITUDES);
methods % Setters
function obj = set.altitudes(obj, value)
obj.ALTITUDES = value;
obj = obj.updateAll();
end
function obj = set.magnitudeDistribution(obj, value)
obj.MAGNITUDE_DISTRIBUTION = value;
obj.isReady = obj.checkProperties();
obj.magnitude = obj.updateDistribution(value, obj.MAGNITUDE_PARAMETERS);
obj.checks(1) = length(obj.magnitude) == length(obj.ALTITUDES);
obj = obj.updateInterpolant();
end
function obj = set.magnitudeParameters(obj, value)
obj.MAGNITUDE_PARAMETERS = value;
obj.isReady = obj.checkProperties();
obj.magnitude = obj.updateDistribution(obj.MAGNITUDE_DISTRIBUTION, value);
obj.checks(1) = length(obj.magnitude) == length(obj.ALTITUDES);
obj = obj.updateInterpolant();
end
function obj = set.azimuthDistribution(obj, value)
obj.AZIMUTH_DISTRIBUTION = value;
obj.isReady = obj.checkProperties();
obj.azimuth = obj.updateDistribution(value, obj.AZIMUTH_PARAMETERS);
obj.checks(2) = length(obj.azimuth) == length(obj.ALTITUDES);
obj = obj.updateInterpolant();
end
function obj = set.azimuthParameters(obj, value)
obj.AZIMUTH_PARAMETERS = value;
obj.isReady = obj.checkProperties();
obj.azimuth = obj.updateDistribution(obj.AZIMUTH_DISTRIBUTION, value);
obj.checks(2) = length(obj.azimuth) == length(obj.ALTITUDES);
obj = obj.updateInterpolant();
end
function obj = set.elevationDistribution(obj, value)
obj.ELEVATION_DISTRIBUTION = value;
obj.isReady = obj.checkProperties();
obj.elevation = obj.updateDistribution(value, obj.ELEVATION_PARAMETERS);
obj.checks(3) = length(obj.elevation) == length(obj.ALTITUDES);
obj = obj.updateInterpolant();
end
function obj = set.elevationParameters(obj, value)
obj.ELEVATION_PARAMETERS = value;
obj.isReady = obj.checkProperties();
obj.elevation = obj.updateDistribution(obj.ELEVATION_DISTRIBUTION, value);
obj.checks(3) = length(obj.elevation) == length(obj.ALTITUDES);
obj = obj.updateInterpolant();
end
function value = get.altitudes(obj), value = obj.ALTITUDES; end
function value = get.magnitudeDistribution(obj), value = obj.MAGNITUDE_DISTRIBUTION; end
function value = get.azimuthDistribution(obj), value = obj.AZIMUTH_DISTRIBUTION; end
function value = get.elevationDistribution(obj), value = obj.ELEVATION_DISTRIBUTION; end
function value = get.magnitudeParameters(obj), value = obj.MAGNITUDE_PARAMETERS; end
function value = get.azimuthParameters(obj), value = obj.AZIMUTH_PARAMETERS; end
function value = get.elevationParameters(obj), value = obj.ELEVATION_PARAMETERS; end
end
methods
......@@ -123,18 +121,8 @@ classdef Wind < Config
obj@Config(mission, varIn);
end
function obj = updateAll(obj)
obj.magnitude = obj.updateDistribution(obj.MAGNITUDE_DISTRIBUTION, obj.MAGNITUDE_PARAMETERS);
obj.azimuth = obj.updateDistribution(obj.AZIMUTH_DISTRIBUTION, obj.AZIMUTH_PARAMETERS);
obj.elevation = obj.updateDistribution(obj.ELEVATION_DISTRIBUTION, obj.ELEVATION_PARAMETERS);
obj.checks(1) = length(obj.magnitude) == length(obj.ALTITUDES);
obj.checks(2) = length(obj.azimuth) == length(obj.ALTITUDES);
obj.checks(3) = length(obj.elevation) == length(obj.ALTITUDES);
end
function [uw, vw, ww] = getVels(obj, z)
if ~all(obj.checks)
if ~obj.isReady
error(['Parameters and distributions must be the same ' ...
'size as altitudes or scalar']);
end
......@@ -144,18 +132,52 @@ classdef Wind < Config
vw = round( - obj.magnitude * sin(obj.azimuth) * cos(obj.elevation), 6);
ww = round( - obj.magnitude * (-sin(obj.elevation)), 6);
else
mag = interpLinear(obj.altitudes, obj.magnitude, z, true);
az = interpLinear(obj.altitudes, obj.azimuth, z, true);
el = interpLinear(obj.altitudes, obj.elevation, z, true);
data = obj.dataInerpolant(z);
mag = data(1); az = data(2); el = data(3);
uw = round( - mag * cos(az) * cos(el), 6);
vw = round( - mag * sin(az) * cos(el), 6);
ww = round( - mag * (-sin(el)), 6);
end
end
function obj = updateAll(obj)
obj.magnitude = obj.updateDistribution(obj.MAGNITUDE_DISTRIBUTION, obj.MAGNITUDE_PARAMETERS);
obj.azimuth = obj.updateDistribution(obj.AZIMUTH_DISTRIBUTION, obj.AZIMUTH_PARAMETERS);
obj.elevation = obj.updateDistribution(obj.ELEVATION_DISTRIBUTION, obj.ELEVATION_PARAMETERS);
obj.isReady = obj.checkProperties();
obj = obj.updateInterpolant();
end
function obj = updateInterpolant(obj)
% UPDATEINTERPOLANT: Interpolates [magnitude; azimuth; elevation] based on altitude
%
% Returns:
% - interpolant: griddedInterpolant for [magnitude; azimuth; elevation]
if ~obj.isReady, return; end
if isscalar(obj.ALTITUDES), return; end
% Prepare data for interpolation
data = [obj.magnitude; obj.azimuth; obj.elevation];
gridVec = obj.ALTITUDES;
% Create gridded interpolant
obj.dataInerpolant = griddedInterpolant(gridVec, data', 'linear', 'nearest');
end
end
methods(Access = private)
function ready = checkProperties(obj)
% Check if STATIC, DYNAMIC, GEOMETRY, and STATE are non-empty
ready = ...
length(obj.magnitude) == length(obj.ALTITUDES) && ...
length(obj.azimuth) == length(obj.ALTITUDES) && ...
length(obj.elevation) == length(obj.ALTITUDES);
end
function vec = updateDistribution(obj, dist, parameters)
s = length(obj.ALTITUDES);
vec = [];
......
......@@ -49,12 +49,12 @@ if not(any(vars.alpha == 1))
error('vars.alpha does not contains 1');
end
if not(any(vars.alpha == -1))
error('vars.alpha does not contains -1');
if isempty(vars.phi) && not(any(vars.alpha == -1))
error('vars.alpha does not contain -1 in alpha-beta configuration');
end
if not(isequal(vars.alpha, -fliplr(vars.alpha)))
error('vars.alpha is not symmetric');
if isempty(vars.phi) && not(isequal(vars.alpha, -fliplr(vars.alpha)))
error('vars.alpha is not symmetric in alpha-beta configuration');
end
% if any(vars.abk > 1)
......@@ -69,6 +69,7 @@ input.fltcon.about = 'Flight conditions quantities';
input.fltcon.MACH = vars.mach;
input.fltcon.ALPHA = vars.alpha;
input.fltcon.BETA = vars.beta;
input.fltcon.PHI = vars.phi;
input.fltcon.ALT = vars.alt;
%% REFQ
......
......@@ -96,11 +96,13 @@ if theta < tb
m = interpLinear(rocket.motor.time, rocket.mass, theta);
T = interpLinear(rocket.motor.time, rocket.motor.thrust, theta);
Pe = interpLinear(rocket.motor.time, rocket.motor.pe, theta);
xcg = interpLinear(rocket.motor.time, rocket.xcg, theta);
T = T + rocket.motor.ae*(Pe - P);
I = interpLinear(rocket.motor.time, rocket.inertia, theta);
Idot = interpLinear(rocket.motor.time, rocket.inertiaDot, theta);
else % for theta >= tb the flight condition is the empty one (no interpolation needed)
xcg = rocket.xcg(end);
m = rocket.cutoffMass;
T = 0;
......@@ -143,53 +145,20 @@ betaOut = beta;
%% INTERPOLATE AERODYNAMIC COEFFICIENTS:
if abs(alpha)>25*pi/180 || abs(beta)>25*pi/180
coeffsValues = interpN( rocket.coefficientsHighAOA.total,...
{rocket.coefficientsHighAOA.state.alphas, rocket.coefficientsHighAOA.state.machs, ...
rocket.coefficientsHighAOA.state.betas, rocket.coefficientsHighAOA.state.altitudes}, ...
[alpha, mach, beta, absoluteAltitude]);
angle0 = [alpha beta];
if abs(alpha) > rocket.coefficients.state.alphas(end)*pi/180 || ...
abs(beta) > rocket.coefficients.state.betas(end)*pi/180
coeffsValues = rocket.coefficientsHighAOA.get(alpha, mach, beta, absoluteAltitude, extension, xcg);
else
%[coeffsValues, angle0] = rocket.interpCoeffs(theta, alpha, mach, beta, absoluteAltitude, opening);
[coeffsValues, angle0] = interpCoeffs(rocket.coefficients, theta, rocket.motor.cutoffTime, ...
alpha, mach, beta, absoluteAltitude, extension);
coeffsValues = rocket.coefficients.get(alpha, mach, beta, absoluteAltitude, extension, xcg);
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);
% 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
CA = coeffsValues(1); CY = coeffsValues(2); CN = coeffsValues(3);
Cl = coeffsValues(4); Cm = coeffsValues(5); Cn = coeffsValues(6);
Clp = coeffsValues(7); Cmad = coeffsValues(8); Cmq = coeffsValues(9);
Cnr = coeffsValues(10); Cnp = coeffsValues(11);
% if Cn == 0 && CY == 0
% XCPlat = -5;
% end
XCPtot = NaN;
%% RAMP / FREE FLIGHT
if altitude < environment.effectiveRampAltitude % No torque on the launchpad
......@@ -209,14 +178,10 @@ if altitude < environment.effectiveRampAltitude % No
betaOut = NaN;
fY = 0;
fZ = 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
......@@ -255,8 +220,6 @@ else
CMaTot = cos(phi)*Cm - sin(phi)*Cn;
if CFaTot ~= 0
XCPtot = CMaTot/CFaTot;
else
XCPtot = XCPlon;
end
dAngle = 0;
......@@ -331,19 +294,12 @@ if nargout == 2 || ~isempty(wrapper)
parout.accelerometer.body_acc = ([-fX+T, fY, -fZ]')/m;
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;
parout.uncertanty = [];
......
......@@ -16,6 +16,11 @@ rocket.inertiaNoMotor = [0.0786; 12.972; 12.972]; % [kg*m^2] OVERRIDE iner
rocket.xCgNoMotor = 1.6025; % [m] OVERRIDE xCg without motor
rocket.lengthCenterNoMot = 1.7640; % [m] OVERRIDE Center length - no nose, no motor
% If dynamic derivatives are loaded, coefficient will depend on rocket xcg
% When false, coefficients are saved with current motor's name
% When true, coefficients are saved as 'generic'
rocket.dynamicDerivatives = false; % [-] True if dynamic derivatives will be loaded
%% PLD - Includes Parafoil + Nose
parafoil = Parafoil();
......
......@@ -16,6 +16,11 @@ rocket.inertiaNoMotor = [0.06; 7.511; 7.512]; % [kg*m^2] OVERRIDE inertia
rocket.xCgNoMotor = 1.1986; % [m] OVERRIDE xCg without motor
rocket.lengthCenterNoMot = 1.7840; % [m] OVERRIDE Center length - no nose, no motor
% If dynamic derivatives are loaded, coefficient will depend on rocket xcg
% When false, coefficients are saved with current motor's name
% When true, coefficients are saved as 'generic'
rocket.dynamicDerivatives = false; % [-] True if dynamic derivatives will be loaded
%% PLD - Includes Parafoil + Nose
parafoil = Parafoil();
......
......@@ -16,6 +16,11 @@ rocket.inertiaNoMotor = [0.06303; 9.62497; 9.62524]; % [kg*m^2] OVERRIDE i
rocket.xCgNoMotor = 1.2239; % [m] OVERRIDE xCg without motor
rocket.lengthCenterNoMot = 1.4470; % [m] OVERRIDE Center length - no nose, no motor
% If dynamic derivatives are loaded, coefficient will depend on rocket xcg
% When false, coefficients are saved with current motor's name
% When true, coefficients are saved as 'generic'
rocket.dynamicDerivatives = false; % [-] True if dynamic derivatives will be loaded
%% PLD - Includes Parafoil + Nose
parafoil = Parafoil();
......
......@@ -18,6 +18,11 @@ rocket.xCgNoMotor = 1.2957; % [m]
% settings.Lcenter - motors(iMotor).L/1000
rocket.lengthCenterNoMot = 1.61; % [m] OVERRIDE Center length - no nose, no motor
% If dynamic derivatives are loaded, coefficient will depend on rocket xcg
% When false, coefficients are saved with current motor's name
% When true, coefficients are saved as 'generic'
rocket.dynamicDerivatives = false; % [-] True if dynamic derivatives will be loaded
%% PLD - Includes Parafoil + Nose
parafoil = Parafoil();
......
......@@ -15,6 +15,12 @@ rocket.inertiaNoMotor = [0.06535397; 17.21019828; 17.21056483]; % [kg*m^2]
rocket.xCgNoMotor = 1.254; % [m] OVERRIDE xCg without motor
rocket.lengthCenterNoMot = 1.517; % [m] OVERRIDE Center length - no nose, no motor
% If dynamic derivatives are loaded, coefficient will depend on rocket xcg
% When false, coefficients are saved with current motor's name
% When true, coefficients are saved as 'generic'
rocket.dynamicDerivatives = false; % [-] True if dynamic derivatives will be loaded
%% PLD - Includes Parafoil + Nose
parafoil = Parafoil();
......
......@@ -15,6 +15,11 @@ rocket.inertiaNoMotor = [0.06535397; 12.07664659; 12.07701314]; % [kg*m^2]
rocket.xCgNoMotor = 1.149; % [m] OVERRIDE xCg without motor
rocket.lengthCenterNoMot = 1.517; % [m] OVERRIDE Center length - no nose, no motor
% If dynamic derivatives are loaded, coefficient will depend on rocket xcg
% When false, coefficients are saved with current motor's name
% When true, coefficients are saved as 'generic'
rocket.dynamicDerivatives = false; % [-] True if dynamic derivatives will be loaded
%% PLD - Includes Parafoil + Nose
parafoil = Parafoil();
......
......@@ -9,6 +9,11 @@ rocket.inertiaNoMotor = []; % [kg*m^
rocket.xCgNoMotor = []; % [m] OVERRIDE xCg without motor
rocket.lengthCenterNoMot = []; % [m] OVERRIDE Center length - no nose, no motor
% If dynamic derivatives are loaded, coefficient will depend on rocket xcg
% When false, coefficients are saved with current motor's name
% When true, coefficients are saved as 'generic'
rocket.dynamicDerivatives = false; % [-] True if dynamic derivatives will be loaded
%% PLD - Includes Parafoil + Nose
parafoil = Parafoil();
......
......@@ -15,6 +15,11 @@ rocket.lengthCenterNoMot = []; % [m]
% rocket.xCgNoMotor = 1.28; % [m] OVERRIDE xCg without motor
% rocket.lengthCenterNoMot = 1.778; % [m] OVERRIDE Center length - no nose, no motor
% If dynamic derivatives are loaded, coefficient will depend on rocket xcg
% When false, coefficients are saved with current motor's name
% When true, coefficients are saved as 'generic'
rocket.dynamicDerivatives = false; % [-] True if dynamic derivatives will be loaded
%% PLD - Includes Parafoil + Nose
parafoil = Parafoil();
......
......@@ -9,6 +9,11 @@ rocket.diameter = 0.15; % [m]
% rocket.xcgNoMotor = []; % [m] OVERRIDE xCg without motor
% rocket.lengthCenterNoMot = []; % [m] OVERRIDE Center length - no nose, no motor
% If dynamic derivatives are loaded, coefficient will depend on rocket xcg
% When false, coefficients are saved with current motor's name
% When true, coefficients are saved as 'generic'
rocket.dynamicDerivatives = false; % [-] True if dynamic derivatives will be loaded
%% PRF - Includes Parafoil + Nose
parafoil = Parafoil();
......@@ -69,7 +74,7 @@ airbrakes.servoTau = 0.0374588;
%% MOTOR
motor = Motor();
motor.name = 'HRE_ARM_OPT_3_Update';
motor.name = 'HRE_ARM_EuRoC_2024';
motor.cutoffTime = []; % [s] OVERRIDE Cutoff time
motor.ignitionTransient = 0.3; % [s] Ignition transient
......
......@@ -15,6 +15,11 @@ rocket.lengthCenterNoMot = []; % [m]
% rocket.xCgNoMotor = 1.28; % [m] OVERRIDE xCg without motor
% rocket.lengthCenterNoMot = 1.778; % [m] OVERRIDE Center length - no nose, no motor
% If dynamic derivatives are loaded, coefficient will depend on rocket xcg
% When false, coefficients are saved with current motor's name
% When true, coefficients are saved as 'generic'
rocket.dynamicDerivatives = false; % [-] True if dynamic derivatives will be loaded
%% PLD - Includes Parafoil + Nose
parafoil = Parafoil();
......