diff --git a/RoccarasoFlight/postprocessing/ASD/.gitkeep b/RoccarasoFlight/postprocessing/ASD/.gitkeep deleted file mode 100644 index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000 diff --git a/RoccarasoFlight/postprocessing/ASD/classes/Bay.m b/RoccarasoFlight/postprocessing/ASD/classes/Bay.m new file mode 100644 index 0000000000000000000000000000000000000000..b7e177537439bdf53a6a95443e04992d00bc8022 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/Bay.m @@ -0,0 +1,19 @@ +classdef Bay < Component & matlab.mixin.Heterogeneous +% Bay: Represents an abstraction layer for all bays +% Arrays of bay subclasses will be represented as [mxn Bay] arrays + + properties(Abstract) + position % [m] OVERRIDE Position, relative to previous component - + length % [m] Total bay length + mass % [kg] Total bay mass + inertia % [kg*m^2] Total bay inertia (Body reference) + xCg % [m] Cg relative to bay upper side + end + + methods (Static, Sealed, Access=protected) + function default = getDefaultScalarElement + default = Payload; + end + end +end + diff --git a/RoccarasoFlight/postprocessing/ASD/classes/Component.m b/RoccarasoFlight/postprocessing/ASD/classes/Component.m new file mode 100644 index 0000000000000000000000000000000000000000..d99c1c59474ca0cf765e0bcca9428f43c96aa5c4 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/Component.m @@ -0,0 +1,115 @@ +classdef Component < Config & matlab.mixin.Copyable +% Component: Represents an abstraction layer for all components +% Components standardize object construction, config and data loading +% +% Constructor: +% - Component: Creates an instance of the Component class. +% Arguments: +% - mission: Mission, mission object, used to find config +% files +% - varsIn: If present, data is read from this variable +% instead of config file +% +% Options: +% +% - elementWise: +% - if true, data is copied element-by-element +% into target class +% - if false, the loaded data substitutes target class + + properties(Abstract, Access = protected) + mission Mission + end + + methods + function obj = Component(mission, varsIn, options) + % Component: Contains and manages all bays + % + % Constructor: + % - Rocket: Creates an instance of the Rocket class. + % Loaded config: rocketConfig.m > rocket (for overrides) + % Loaded data: - + % Arguments: + % - mission: Mission, mission object + % - varIn: (optional) config source. Alternative to config.m + % file + + arguments (Input) + mission Mission = Mission() + varsIn = [] + options.elementWise logical = true + end + + if isempty(mission.name) + return; + end + + if isempty(varsIn) + varsIn = obj.getConfig(mission); + end + + fileName = obj.configName; + varName = obj.variableName; + if isempty(varName) + varName = strtok(fileName, 'C'); + end + + if ~isfield(varsIn, varName) + error(['Class not found inside the config folder: %s\n' ... + 'With name: %s\n' ... + 'Check that the correct mission is set in your component ' ... + 'and the config file is correct'], fileName, varName); + end + + configObj = varsIn.(varName); + + if options.elementWise + fields = configObj.getProperties('writable'); + for field = fields + if isempty([configObj.(field)]), continue; end + obj.(field) = configObj.(field); + end + else + obj = configObj; + end + + [obj.mission] = deal(mission); + end + end + + methods(Access = protected) + function cp = copyElement(obj) + fields = obj.getProperties('readable', NonCopyable = 0); + cp = copyElement@matlab.mixin.Copyable(obj); %shallow copy of all elements + + for field = fields + if isa(obj.(field), 'matlab.mixin.Copyable') + cp.(field) = copy(obj.(field)); + end + end + end + + function outputVariables = getConfig(obj, mission) + fileName = obj.configName; + filePath = mission.configPath; + file = fullfile(filePath, fileName); + + if ~isfile(file) + error(['File not found inside the config folder: %s\n' ... + 'Check that the correct mission is set in your component ' ... + 'and the config file exists'], fileName); + end + + variables = [who(); "variables"]; + run(file); + loadedVariables = who(); + loadedVariables(ismember(loadedVariables, variables)) = []; + + outputVariables = struct(); + for i = 1:length(loadedVariables) + outputVariables.(loadedVariables{i}) = eval(loadedVariables{i}); + end + end + end +end + diff --git a/RoccarasoFlight/postprocessing/ASD/classes/Config.m b/RoccarasoFlight/postprocessing/ASD/classes/Config.m new file mode 100644 index 0000000000000000000000000000000000000000..cec711a4737ba8e95eb0bc2c54ae6732da94de76 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/Config.m @@ -0,0 +1,136 @@ +classdef(Abstract) Config < handle +% Config: Represents an abstraction layer for all config-dependent Classes +% Config standardizes property management and declares fundemental +% methods (see doc for more info) + + + properties(Abstract, Access = protected) + configName {mustBeTextScalar} + variableName char + end + + methods + function structOut = toStruct(obj) + structOut = struct(); + fields = obj.getProperties('readable'); + for j = 1:size(fields, 2), structOut.(fields{j}) = obj.(fields{j}); end + end + + function [propertiesOut, valuesOut] = getProperties(obj, preset, options, attributes) + % This function returns properties that match specific attributes + % You can specify presets to include combinations of properties: + % - 'readable': GetAccess = 'public', Abstract = false, + % Hidden = false + % - 'writable': SetAccess = 'public', Dependent = false, + % Constant = false, Abstract = false, Hidden = false + % + % Presets have the priority over options, in case of conflicts + % + % WARNING: If second output is used, only values with public + % getAccess will be read + + arguments + obj + preset char {mustBeMember(preset, {'readable', 'writable', ''})} = '' + options.Superclass char = '' + options.Heterogeneous logical = 0 + attributes.GetAccess char {mustBeMember(attributes.GetAccess, {'public', 'private', ''})} = '' + attributes.SetAccess char {mustBeMember(attributes.SetAccess, {'public', 'private', ''})} = '' + attributes.Dependent logical {mustBeMember(attributes.Dependent, [0, 1])} = [] + attributes.Constant logical {mustBeMember(attributes.Constant, [0, 1])} = [] + attributes.Abstract logical {mustBeMember(attributes.Abstract, [0, 1])} = [] + attributes.Transient logical {mustBeMember(attributes.Transient, [0, 1])} = [] + attributes.Hidden logical {mustBeMember(attributes.Hidden, [0, 1])} = [] + attributes.AbortSet logical {mustBeMember(attributes.AbortSet, [0, 1])} = [] + attributes.NonCopyable logical {mustBeMember(attributes.NonCopyable, [0, 1])} = [] + end + + if ~isempty(preset) + switch preset + case 'readable' + attributes.GetAccess = 'public'; + attributes.Abstract = 0; + attributes.Hidden = 0; + case 'writable' + attributes.SetAccess = 'public'; + attributes.Dependent = 0; + attributes.Constant = 0; + attributes.Abstract = 0; + attributes.Hidden = 0; + end + end + if nargout == 2, attributes.GetAccess = 'public'; end + + fields = fieldnames(attributes); + attributeCheck = @(field, prop) isempty(attributes.(field)) ... + || isequal(prop.(field), attributes.(field)); + + optionCheck = @(prop) isempty(options.Superclass) ... + || ~isempty(prop.Validation) ... + && ~isempty(prop.Validation.Class) ... + && ~isempty(prop.Validation.Class.SuperclassList) ... + && any(strcmp({prop.Validation.Class.SuperclassList.Name}, options.Superclass)); + + mc = metaclass(obj); + ii = 0; + nProperties = length(mc.PropertyList); + properties = cell(1,nProperties); + for c = 1:nProperties + mp = mc.PropertyList(c); + checks = [cellfun(@(field) attributeCheck(field, mp), fields); + optionCheck(mp)]; + if all(checks) + ii = ii + 1; + properties(ii) = {mp.Name}; + end + end + propertiesOut = string(properties(1:ii)); + + if nargout == 2 + if isempty(options.Superclass) || ~options.Heterogeneous + valuesOut = cell(1, ii); + for i = 1:ii, valuesOut{i} = obj.(propertiesOut(i)); end + else + valuesOut(ii) = obj.(propertiesOut(1)); + for i = 1:ii, valuesOut(i) = obj.(propertiesOut(i)); end + end + end + end + end + + methods(Hidden) + function lh = addlistener(varargin) + lh = addlistener@handle(varargin{:}); + end + function notify(varargin) + notify@handle(varargin{:}); + end + function delete(varargin) + delete@handle(varargin{:}); + end + function Hmatch = findobj(varargin) + Hmatch = findobj@handle(varargin{:}); + end + function p = findprop(varargin) + p = findprop@handle(varargin{:}); + end + function TF = eq(varargin) + TF = eq@handle(varargin{:}); + end + function TF = ne(varargin) + TF = ne@handle(varargin{:}); + end + function TF = lt(varargin) + TF = lt@handle(varargin{:}); + end + function TF = le(varargin) + TF = le@handle(varargin{:}); + end + function TF = gt(varargin) + TF = gt@handle(varargin{:}); + end + function TF = ge(varargin) + TF = ge@handle(varargin{:}); + end + end +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/classes/DataWrapper.m b/RoccarasoFlight/postprocessing/ASD/classes/DataWrapper.m new file mode 100644 index 0000000000000000000000000000000000000000..6d484e266694e71cfd7d196163a7167963f06402 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/DataWrapper.m @@ -0,0 +1,122 @@ +classdef DataWrapper < handle + % OutputWrapper: A fast way to exchange data by pointer + % + % Constructor: + % - OutputWrapper: Creates an instance of the Wrapper class. + % Arguments: + % - reference, data used to preallocate cache with empty fields + % - chunkSize, default: 1000, space to preallocate each + % time cache is full + % - shift, default: 0, leaves n empty fields at the + % beginning of cache if needed + % + % WARNING: Intended for data exchange where output arguments are inaccessible. + % AVOID using such method if output arguments can be read to + % improve code maintainability + + properties(SetAccess = private) + data struct % Cached data + dummy struct % Empty struct for field reference + counter int32 % Current size of data stored + end + + properties(Access = private) + shift int32 % Number of fields to skip in cache, useful to preallocate initial values / starting conditions + chunkSize int32 % Size of preallocated memory + dataLen int32 % Current preallocated size of cache + end + + methods + function obj = DataWrapper(reference, chunkSize, shift) + arguments + reference struct + chunkSize {mustBeInteger, mustBePositive} = 1000 + shift = 0 + end + obj.chunkSize = chunkSize; + obj.shift = shift; + + fields = fieldnames(reference)'; fields{2, 1} = []; + obj.dummy = struct(fields{:}); + + obj.reset(); + end + + function setCache(obj, data) + % setCache - Saves data as last element in cache (unpromoted) + % When calling setCache, any unpromoted data will be + % overridden + % To avoid overwriting, promote data using setData + + obj.data(obj.counter + 1) = data; + end + + function out = getData(obj, startIdx) + % detData - Reads data from wrapper + % This action will NOT reset cache + arguments + obj + startIdx = 1 + end + out = obj.data(startIdx:obj.counter); + end + + function setData(obj) + % setData - Promoted current data stored in cache + % promoted data will be saved until wrapper is reset + + obj.counter = obj.counter + 1; + if obj.counter == obj.dataLen + len = obj.dataLen + obj.chunkSize; + obj.dataLen = len; + obj.data(len) = obj.dummy; + end + end + + function out = popData(obj) + % popData - Reads data from wrapper + % This action WILL reset cache + + out = obj.packageData(obj.data(1:obj.counter)); + obj.reset(); + end + + function obj = reset(obj) + % reset - Resets cache to starting conditions + + obj.data = obj.dummy; % Set struct dataType + obj.data(obj.chunkSize) = obj.dummy; % Preallocate struct + + obj.dataLen = obj.chunkSize; + obj.counter = obj.shift; + end + end + + methods(Static) + function out = packageData(data) + % Converts struct array to struct with array fields + % + % Main assumption: fields to concatenate are colums vectors or + % scalars + + fields = fieldnames(data); + + for i = 1:length(fields) + currentField = fields{i}; + currentValue = [data.(currentField)]; + + if isstruct(currentValue) + out.(currentField) = DataWrapper.packageData(currentValue); + else + sz = size(data(1).(currentField)); + if all(sz > 1) + szIn = num2cell(sz); szIn{end+1} = []; %#ok<AGROW> + currentValue = reshape(currentValue, szIn{:}); + end + out.(currentField) = currentValue; + end + end + end + end +end + diff --git a/RoccarasoFlight/postprocessing/ASD/classes/Mission.m b/RoccarasoFlight/postprocessing/ASD/classes/Mission.m new file mode 100644 index 0000000000000000000000000000000000000000..f505cd8248e23e677770f3a0954924074d295f9d --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/Mission.m @@ -0,0 +1,90 @@ +classdef Mission < Config + % Mission: Contains names and paths needed to access mission-dependent + % config files + % + % Constructor: + % - Mission: Creates an instance of the Mission class. + % Loaded config: missionConfig.m + % Loaded data: - + % Arguments: + % - loadConfig: bool, Whether to load config file or return + % an empty mission + + properties + name char % Mission name, used to access <mission> folder + end + + properties(SetAccess = protected) + currentPath char + configPath char + dataPath char + msaPath char + end + + properties(Access = protected) + configName = 'missionConfig.m' + variableName = '' + end + + methods + function obj = Mission(mission) + arguments + mission {mustBeA(mission, {'logical', 'string', 'char'})} = false + end + + if nargin == 0 + return; + end + + filePath = fullfile(fileparts(mfilename("fullpath")), '..', 'missions'); + filePath = trimPath(filePath); + + obj.currentPath = filePath; + if isa(mission, 'char') || isa(mission, 'string') + obj.name = mission; + else + obj.loadConfig(); + end + end + + function set.name(obj, name) + if ~isempty(obj.currentPath) && ... + ~isfolder(fullfile(obj.currentPath, name)) %#ok<MCSUP> + error('Invalid mission: %s', name) + end + + obj.name = name; + obj.updatePath(); + end + end + + methods(Access = protected) + function updatePath(obj) + if isempty(obj.name) || isempty(obj.currentPath), return; end + + obj.configPath = fullfile(obj.currentPath, obj.name, 'config'); + obj.dataPath = fullfile(obj.currentPath, obj.name, 'data'); + obj.msaPath = trimPath( ... + fullfile(obj.currentPath, '..', '..', '..', 'msa-toolkit') ... + ); + end + + function loadConfig(obj) + fileName = obj.configName; + filePath = obj.currentPath; + + if ~isfile(fullfile(filePath, fileName)) + error(strcat("File not found inside the config folder: ", filePath)); + end + + varName = strtok(fileName,'C'); + run(fileName); + configObj = eval(varName); + obj.name = configObj.name; + end + end + + methods(Static, Access = protected) + function loadData(), end + end +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/classes/Rocket.m b/RoccarasoFlight/postprocessing/ASD/classes/Rocket.m new file mode 100644 index 0000000000000000000000000000000000000000..3f3e895253d0cd6b737a808f94809a54708b082e --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/Rocket.m @@ -0,0 +1,398 @@ +classdef Rocket < Component + % Rocket: Contains and manages all bays + % + % Constructor: + % - Rocket: Creates an instance of the Rocket class. + % Loaded config: rocketConfig.m > rocket (for overrides) + % Loaded data: - + % Arguments: + % - mission: Mission, mission object + % - varIn: (optional) config source. Alternative to config.m + % file + + properties + payload Payload % [-] Payload bay + recovery Recovery % [-] Recovery bay + electronics Electronics % [-] Electronics bay + airbrakes Airbrakes % [-] Airbrakes bay + motor Motor % [-] Motor bay + rear Rear % [-] Rear bay + + pitot Pitot % [-] Pitot component + parachutes Para % [-] (nParachutes, nStages) Parachutes onboard + + length double % [m] Total length + lengthCenter double % [m] Center length - no nose, boat + lengthCenterNoMot double % [m] Center length - no nose, no motor + diameter double % [m] Diameter + mass double % [kg] Total mass + massNoMotor double % [kg] Mass without motor + inertia double % [kg*m^2] Total inertia - Axibody reference + inertiaNoMotor double % [kg*m^2] Inertia without motor + xCg double % [m] Total xCg + xCgNoMotor double % [m] xCg without motor + + stagesMass double % [kg] Mass of stage 2, 3, .. without chutes + + coefficients % [-] Aerodynamic coefficients + coefficientsHighAOA % [-] Aerodynamic coefficients for high angle of attacks + end + + properties(SetAccess = protected) + crossSection % [m^2] Rocket cross sectional area + inertiaDot double % [kg*m^2/s]Total inertia time derivative + cutoffInertia double % [kg*m^2] Total inertia at motor cutoff + cutoffMass double % [m] Total mass at motor cutoff + absolutePositions dictionary ... % [m] Bay start positions - 0 is set at nose base + = dictionary() + bays dictionary ... % [Bay] All bays + = dictionary() + end + + properties(Access = protected, Hidden) + configName = 'rocketConfig.m' + variableName = '' + mission Mission + end + + methods % Updaters + function updateAbsolutePositions(obj) + bDictionay = obj.bays; + b = values(bDictionay); + + if ~isempty(obj.lengthCenterNoMot) + obj.absolutePositions = dictionary(["motor", "rear"], ... + [obj.lengthCenterNoMot, obj.lengthCenterNoMot + obj.rear.position]); + return + end + % Populate shift vectors (how much each bay shifts the next one back) + % Shift overrides are applyed to the previous bay + shift = zeros(size(b)); + overrides = ~cellfun(@isempty, {b.position}); + overridesShift = logical([overrides(2:end), 0]); + + shift(overridesShift) = b(overrides).position; + shift(~overridesShift) = [b(~overridesShift).length]; + + absPositions = cumsum([-shift(1); shift(1:end-1)]); + absPositions = absPositions + ... + (obj.payload.length - obj.payload.noseLength); + obj.absolutePositions = dictionary(keys(bDictionay), absPositions); + + end + function updateGeometry(obj) + if ~isempty(obj.lengthCenterNoMot) + obj.length = obj.lengthCenterNoMot + obj.motor.length + obj.payload.noseLength; + else + obj.length = (obj.absolutePositions("rear") + obj.bays("rear").length) + obj.bays("payload").noseLength; + obj.lengthCenterNoMot = obj.absolutePositions("motor"); + end + if isempty(obj.crossSection) + obj.crossSection = 0.25*pi*obj.diameter^2; + end + obj.lengthCenter = obj.length - obj.payload.noseLength - obj.rear.boatLength; + end + + function updateMassNoMotor(obj) + if ~isempty(obj.massNoMotor) + return; + end + bNoMotor = values(remove(obj.bays, 'motor')); + obj.massNoMotor = sum([bNoMotor.mass]); + end + + function updateMass(obj) + obj.mass = obj.massNoMotor + obj.motor.mass; + end + + function updateXCgNoMotor(obj) + if ~isempty(obj.xCgNoMotor) + return; + end + bNoMotor = values(remove(obj.bays, 'motor'))'; + aPosNoMotor = values(remove(obj.absolutePositions, 'motor'))'; + obj.xCgNoMotor = ((aPosNoMotor + [bNoMotor.xCg])* ... + [bNoMotor.mass]')/obj.massNoMotor; + end + + function updateXCg(obj) + obj.xCg = (obj.xCgNoMotor*obj.massNoMotor + ... + obj.motor.mass.*(obj.motor.xCg+obj.lengthCenterNoMot))./ ... + obj.mass; + end + + function updateInertiaNoMotor(obj) + if ~isempty(obj.inertiaNoMotor) + return; + end + bNoMotor = values(remove(obj.bays, 'motor'))'; + aPosNoMotor = values(remove(obj.absolutePositions, 'motor'))'; + % [3x1] Ix, Iy, Iz + % Assumption: xCgs are close to rocket axis + inertias = [bNoMotor.inertia]; + temp = inertias(2:3, :) + ... + (obj.xCgNoMotor' - (aPosNoMotor + [bNoMotor.xCg])).^2 .* [bNoMotor.mass]; + obj.inertiaNoMotor = [sum(inertias(1, :)); sum(temp, 2)]; + end + + function updateInertia(obj) + % [3xN] Ix, Iy, Iz + % Assumption: xCgs are close to rocket axis + motorInertia = obj.motor.inertia(2:3, :) + ... + (obj.xCg - (obj.absolutePositions('motor') + obj.motor.xCg)).^2 .* obj.motor.mass; + baysInertia = obj.inertiaNoMotor(2:3, :) + ... + (obj.xCg - obj.xCgNoMotor).^2 .* obj.massNoMotor; + + obj.inertia = [obj.motor.inertia(1, :) + obj.inertiaNoMotor(1, :); + motorInertia + baysInertia]; + + obj.inertiaDot = diff(obj.inertia, 1, 2)./diff(obj.motor.time); + obj.inertiaDot(:, end + 1) = obj.inertiaDot(:, end); + end + + function updateCutoff(obj) + obj.cutoffInertia = interpLinear(obj.motor.time, obj.inertia, obj.motor.cutoffTime); + obj.cutoffMass = interpLinear(obj.motor.time, obj.mass, obj.motor.cutoffTime); + end + + function updateStagesMass(obj) + stage1 = obj.cutoffMass - (obj.payload.mass + obj.parachutes(1,2).mass + obj.parachutes(2,2).mass); + % Everything at cut off without payload, payload drogue and + % payload airfoil + stage2 = obj.payload.mass + obj.parachutes(1,2).mass + obj.parachutes(2,2).mass; + % only payload: payload, payload drogue and payload airfoil + obj.stagesMass = [stage1 stage2]; + end + + function updateAll(obj) + % Note: properties without motor must be updated first + obj.updateAbsolutePositions; + obj.updateGeometry; + obj.updateMassNoMotor; + obj.updateMass; + obj.updateXCgNoMotor; + obj.updateXCg; + obj.updateInertiaNoMotor; + obj.updateInertia; + obj.updateCutoff; + obj.updateStagesMass; + end + end + + methods(Access = protected) % Loaders + function obj = loadData(obj) + if isempty(obj.motor) + return; + end + varNames = {'total', 'geometry', 'state', 'finsCN'}; + motorName = obj.motor.name; + aeroPath = fullfile(obj.mission.dataPath, 'aeroCoefficients.mat'); + aeroHighAOAPath = fullfile(obj.mission.dataPath, 'aeroCoefficientsHighAOA.mat'); + + if ~(exist(aeroPath, 'file') && exist(aeroHighAOAPath,'file')) + return; + end + dataCoeffs = load(aeroPath); + dataCoeffsHighAOA = load(aeroHighAOAPath); + + if isfield(dataCoeffs, motorName) && isfield(dataCoeffsHighAOA, motorName) + dataCoeffs = dataCoeffs.(motorName); + dataCoeffsHighAOA = dataCoeffsHighAOA.(motorName); + else + return; + end + + if isfield(dataCoeffs, "finsCN") && ... + isfield(dataCoeffsHighAOA, "finsCN") + obj.coefficients.finsCN = dataCoeffs.finsCN; + obj.coefficientsHighAOA.finsCN = dataCoeffsHighAOA.finsCN; + end + + varNames = varNames{1:end - 1}; + if all(and(isfield(dataCoeffs, varNames), ... + isfield(dataCoeffsHighAOA, varNames))) + obj.coefficients.total = dataCoeffs.total; + obj.coefficientsHighAOA.total = dataCoeffsHighAOA.total; + obj.coefficients.geometry = dataCoeffs.geometry; + obj.coefficientsHighAOA.geometry = dataCoeffsHighAOA.geometry; + obj.coefficients.state = dataCoeffs.state; + obj.coefficientsHighAOA.state = dataCoeffsHighAOA.state; + end + end + end + + methods + function obj = Rocket(mission, varIn, options) + arguments + mission Mission = Mission() + varIn struct = [] + options.loadCoefficients logical = true + options.checkGeometry logical = true + end + obj@Component(mission, varIn); + + %% Loading data + if isempty(mission.name) + return; + end + + vars = obj.getConfig(mission); % Load config once and apply to other bays + obj.payload = Payload(mission, vars); + obj.recovery = Recovery(mission, vars); + obj.electronics = Electronics(mission, vars); + obj.motor = Motor(mission, vars); + obj.airbrakes = Airbrakes(mission, obj.motor, vars); + obj.rear = Rear(mission, vars); + obj.pitot = Pitot(mission, vars); + paras = Para(mission); + obj.parachutes = paras(:, :, 1); + + [bayNames, bayValues] = obj.getProperties(Superclass='Bay', Heterogeneous=1); + obj.bays = dictionary(bayNames, bayValues); + obj.updateAll(); + + if options.loadCoefficients + obj.loadData(); + answer = ''; + if isempty(obj.coefficients) || isempty(obj.coefficientsHighAOA) + answer = questdlg(['Coefficient matrices not found. ' ... + 'Do you want to create new matrices?']); + elseif options.checkGeometry && ~all(obj.checkGeometry()) + answer = questdlg(['Coefficient matrices differ from rocket geometry. ' ... + 'Do you want to create new matrices?']); + end + + switch answer + case 'Yes' + parserPath = fullfile(mission.msaPath, 'autoMatricesProtub'); + addpath(genpath(parserPath)); + mainAutoMatProtub(obj); + obj.loadData(); + case 'Cancel' + error('Rocket creation aborted') + otherwise + end + end + end + + function [coeffsValues, angle0] = interpCoeffs(obj, t, alpha, mach, beta, alt, c) + % interpCoeffs - interpolation of aerodynamic coefficients. + % + % INPUTS: + % - t, double [1,1], integration time, [s]; + % - alpha, double[1,1], angle of attack, []; + % - mach, double[1,1], mach number, []; + % - beta, double[1,1], sideslip angle, []; + % - alt, double[1,1], altitude, [m]; + % - c, double[1,1], aerobrakes control variable, []; + % + % OUTPUTS: + % - coeffsValues, array [16,1], aerodynamic coefficients; + % - angle0, array [2,1], array of the reference aerodynamic angles. + + if mach > obj.airbrakes.maxMach + c = 0; + end + h = c*obj.coefficients.state.hprot(end); + + [nearAlpha, indexAlpha] = nearestValSorted(obj.coefficients.state.alphas*pi/180, alpha); + [~, indexMach] = nearestValSorted(obj.coefficients.state.machs, mach); + [nearBeta, indexBeta] = nearestValSorted(obj.coefficients.state.betas*pi/180, beta); + [~, indexAlt] = nearestValSorted(obj.coefficients.state.altitudes, alt); + [~, indexControl] = nearestValSorted(obj.coefficients.state.hprot, h); + + angle0 = [nearAlpha; nearBeta]; + + if t >= obj.motor.cutoffTime + % Interpolate on airbrake control + % Uses end xCg, even if tCutOff < tBurnTime + num = length(obj.coefficients.state.hprot); %#ok<CPROPLC> + index1 = indexControl; + index2 = indexControl + 1; + if index2 > num + if num == 1 + coeffsValues = obj.coefficients(:, indexAlpha, indexMach, indexBeta, indexAlt, 1, end); + return; + else + index2 = index2 - 1; + index1 = index1 - 1; + end + end + + VF = obj.coefficients.total(:, indexAlpha, indexMach, indexBeta, indexAlt, index1, end); + VE = obj.coefficients.total(:, indexAlpha, indexMach, indexBeta, indexAlt, index2, end); + + deltaH = obj.coefficients.state.hprot(index2) - obj.coefficients.state.hprot(index1); + coeffsValues = ( (h - obj.coefficients.state.hprot(index1))/deltaH)*(VE - VF) + VF; + else + % Interpolate on xCg, close airbrakes + num = 1:1:length(obj.coefficients.state.xcgTime)-1; %#ok<CPROPLC> + indexXcg = t >= obj.coefficients.state.xcgTime(1:end-1) & t < obj.coefficients.state.xcgTime(2:end); + index1 = num(indexXcg); + index2 = index1 + 1; + + VF = obj.coefficients.total(:, indexAlpha, indexMach, indexBeta, indexAlt, indexControl, index1); + VE = obj.coefficients.total(:, indexAlpha, indexMach, indexBeta, indexAlt, indexControl, index2); + + deltaT = obj.coefficients.state.xcgTime(index2) - obj.coefficients.state.xcgTime(index1); + coeffsValues = ( (t - obj.coefficients.state.xcgTime(index1))/deltaT)*(VE - VF) + VF; + end + end + + function checks = checkGeometry(obj) + % checkGeometry - This methods checks if the rocket geometry + % is consistent with the geometry of the + % aerodynamic matrices + % + % OUTPUTS: + % - checks (n fields of geometry, 1): boolean value of the geometry checks + + xCgRocket = round([ ... + obj.xCg(1); ... + obj.xCg(end) ... + ], 3); + + xCgTest = round([ + obj.coefficients.geometry.xcg(1); + obj.coefficients.geometry.xcg(end); + ], 3); + + geometryRocket = round([ + obj.diameter; + obj.payload.noseLength; + obj.lengthCenter; + obj.payload.noseCMod; + obj.payload.nosePMod; + obj.rear.finsRootChord; + obj.rear.finsFreeChord; + obj.rear.finsHeight; + obj.rear.finsDeltaXFreeChord; + obj.rear.nPanel; + obj.rear.boatLength; + obj.rear.boatFinalDiameter; + ], 3); + + 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); + + checks = [ + xCgRocket == xCgTest; + strcmp(obj.coefficients.geometry.ogType, obj.payload.noseType); + geometryRocket == geometryTest; + strcmp(obj.coefficients.geometry.boatType, obj.rear.boatType) + ]; + end + end +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/classes/Settings.m b/RoccarasoFlight/postprocessing/ASD/classes/Settings.m new file mode 100644 index 0000000000000000000000000000000000000000..9afdc2c433dea2d28665cdf3ec6cc080e5eb0919 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/Settings.m @@ -0,0 +1,85 @@ +classdef Settings < Config & dynamicprops +% Settings: Provides standardized way of loading config files +% Looks for files in missions > settings for global settings +% Looks for files in caller's folder otherwise +% +% Constructor: +% - Settings: Creates an instance of the Settings class. +% Loaded config: - +% Loaded data: - +% Arguments: +% - configNames: config file names. Accepts shortened version +% e.g: ode -> odeConfig.m + + properties(Access = protected) + configPath = '' + configName = '' + variableName = '' + end + + methods + function obj = Settings(configNames) + arguments (Input, Repeating) + configNames char + end + + configNames = string(configNames); + + noConfig = ~endsWith(configNames, ["Config.m", "Config"]); + noExtension = ~endsWith(configNames, ".m"); + + configNames(noConfig) = strcat(configNames(noConfig), 'Config'); + configNames(noExtension) = strcat(configNames(noExtension), '.m'); + + varsIn = obj.getConfig(configNames); + obj.loadConfig(varsIn); + end + end + + methods(Access = protected) + function loadConfig(obj, varsIn) + props = fields(varsIn); + for j = 1:length(props) + obj.addprop(props{j}); + obj.(props{j}) = varsIn.(props{j}); + end + end + end + + methods(Static) + function outputVariables = getConfig(configNames) + outputVariables = struct(); + for i = 1:length(configNames) + fileName = configNames{i}; + filePath = fullfile(fileparts(mfilename("fullpath")), '..', 'settings'); + + if ~isfile(fullfile(filePath, fileName)) % Get caller function path + callers = dbstack("-completenames"); + if length(callers) > 2 + filePath = fileparts(callers(3).file); + end + end + + if ~isfile(fullfile(filePath, fileName)) + error(['File not found inside the settings folder, nor in local folder: %s\n' ... + 'Check that the correct mission is set in your component ' ... + 'and the config file exists'], fileName); + end + + variables = [who(); "variables"]; + run(fileName); + loadedVariables = who(); + loadedVariables(ismember(loadedVariables, variables)) = []; + + for j = 1:length(loadedVariables) + outputVariables.(loadedVariables{j}) = eval(loadedVariables{j}); + end + end + end + end + + + methods(Static, Access = protected) + function loadData(), end + end +end diff --git a/RoccarasoFlight/postprocessing/ASD/classes/bays/Airbrakes.m b/RoccarasoFlight/postprocessing/ASD/classes/bays/Airbrakes.m new file mode 100644 index 0000000000000000000000000000000000000000..dbb8478838af27fba4ea60be0c2c1bb4c559b74e --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/bays/Airbrakes.m @@ -0,0 +1,87 @@ +classdef Airbrakes < Bay +% Airbrakes: Represents airbrakes configuration for a rocket. +% +% Constructor: +% - Airbrakes: Creates an instance of the Airbrakes class. +% Loaded config: rocketConfig.m > airbrakes +% Loaded data: - +% Arguments: +% - mission: Mission, mission object +% - varIn: (optional) config source. Alternative to config.m +% file + + properties + enabled logical % If false, airbrakes will be kept closed during flight + identification logical % If false, + + n double % [-] number of brakes + thickness double % [m] brakes thickness + width double % [m] brakes width + xDistance double % [m] axial position from nosecone base + + servoMaxAngle double + servoMinAngle double + servoTau double + servoOmega double % [rad/s] Servo-motor angular velocity + + extension double % aerobrakes, 1-2-3 for 0%, 50% or 100% opened + deltaTime double % aerobrakes, configurations usage time + maxMach double % [-] Maximum Mach at which airbrakes can be used + height double % [m] Block airbrakes opening coordinate ( First entry must be 0! ) + angleFunction function_handle % [-] Relation between angle and extension height + + position % [m] Offset with respect to other bays (negative -> shift forward) + length % [m] Total bay length + mass % [kg] Total bay mass + inertia % [kg*m^2] Total bay inertia (Body reference) + xCg % [m] Cg relative to bay upper side + end + + properties(Access = protected) + configName = 'rocketConfig.m' + variableName = 'airbrakes' + mission Mission = Mission() + motor Motor = Motor() + end + + properties(Access = public) + time double + end + + methods + function obj = Airbrakes(mission, motor, varIn) + arguments(Input) + mission Mission = Mission() + motor Motor = Motor() + varIn = [] + end + obj@Bay(mission, varIn); + obj.motor = motor; + end + + function set.motor(obj, motor) + obj.motor = motor; + obj.updateTime; + end + end + + methods + function updateTime(obj) + obj.time = cumsum(obj.deltaTime) + obj.motor.cutoffTime; + end + + function extension = angle2Extension(obj, angle) + extension = obj.angleFunction(angle) / obj.height(end); + end + + function extension = getExtension(obj, time, mach) + if time < obj.time(1) || ~obj.enabled || mach > obj.maxMach + extension = 0; + return; + end + + idx = sum(time >= obj.time); + extension = obj.extension(idx); + end + end +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/classes/bays/Electronics.m b/RoccarasoFlight/postprocessing/ASD/classes/bays/Electronics.m new file mode 100644 index 0000000000000000000000000000000000000000..a9448aa834efea6ca915ddce4ebf3b4b97b7e4aa --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/bays/Electronics.m @@ -0,0 +1,36 @@ +classdef Electronics < Bay +% Electronics: Represents electronics configuration for a rocket. +% +% Constructor: +% - Electronics: Creates an instance of the Electronics class. +% Loaded config: rocketConfig.m > electronics +% Loaded data: - +% Arguments: +% - mission: Mission, mission object +% - varIn: (optional) config source. Alternative to config.m +% file + properties + offset = 0 + position % [m] Absolute position, relative to nose base + length % [m] Total bay length + mass % [kg] Total bay mass + inertia % [kg*m^2] Total bay inertia (Body reference) + xCg % [m] Cg relative to bay upper side + end + + properties(Access = protected) + configName = 'rocketConfig.m' + variableName = 'electronics' + mission Mission = Mission() + end + + methods + function obj = Electronics(mission, varIn) + arguments(Input) + mission Mission = Mission() + varIn = [] + end + obj@Bay(mission, varIn); + end + end +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/classes/bays/Motor.m b/RoccarasoFlight/postprocessing/ASD/classes/bays/Motor.m new file mode 100644 index 0000000000000000000000000000000000000000..3a8af1ffa1d9645fb741fcd2a51d53f0412686a4 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/bays/Motor.m @@ -0,0 +1,120 @@ +classdef Motor < Bay + % Motor: Represents the motor configuration for a rocket. + % + % Constructor: + % - Motor: Creates an instance of the Motor class. + % Loaded config: rocketConfig.m > motor + % Loaded data: motors.mat + % Arguments: + % - mission: Mission, mission object + % - varIn: (optional) config source. Alternative to config.m file + + properties + name {mustBeTextScalar} = '' % [-] Motor name + position % [m] Absolute position, relative to nose base + length % [m] Total length (motor + tank) + inertia % [kg*m^2] Total Motor inertia (Body reference) + cutoffTime double % [s] Shutdown time + tankLength double % [m] Tank length + ignitionTransient double % [s] Ignition transient duration + cutoffTransient 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 % [m] Engine xcg from tank tip + pe double % [Pa] Eflux pressure + ae double % [Pa] Eflux Area + end + + properties + mass % [kg] Total Motor mass + fuselageXCg double % [m] xcg of the engine fuselage only from tank tip + end + + properties(Access = protected) + configName = 'rocketConfig.m' + variableName = 'motor' + mission Mission = Mission() + end + + properties(SetAccess = protected) + isHRE logical % [-] Flag relateed to the type of motor: true if HRE + end + + methods + function obj = Motor(mission, varIn) + arguments(Input) + mission Mission = Mission() + varIn = [] + end + obj@Bay(mission, varIn); + obj.loadData(); + obj.updateAll(); + end + + function set.name(obj, name) + obj.name = name; + obj.loadData(); + obj.updateAll(); + end + + function updateAll(obj) + obj.updateMass(); + end + + function updateMass(obj) + obj.fuselageXCg = (obj.length - ... + obj.tankLength)/2 + obj.tankLength; + obj.mass = obj.propellantMass + obj.structureMass; + end + end + + methods (Access = protected) + 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(obj.name) || isvarname(obj.name)) + error('Motor name is not a valid MATLAB variable name') + end + + if isempty(chosenMotor) + error(strcat('Unable to find engine: ', obj.name)); + end + + obj.isHRE = contains(obj.name, 'HRE'); + + if obj.isHRE + obj.length = chosenMotor.L; + obj.tankLength = chosenMotor.Ltank; + obj.time = chosenMotor.t; + obj.thrust = chosenMotor.T; + obj.fuelMass = chosenMotor.mFu; + obj.oxidizerMass = chosenMotor.mOx; + obj.structureMass = chosenMotor.mc; + obj.propellantMass = chosenMotor.m; + obj.inertia = [chosenMotor.Ixx;chosenMotor.Iyy;chosenMotor.Izz]; + obj.xCg = chosenMotor.xcg; + obj.pe = chosenMotor.Pe; + obj.ae = chosenMotor.Ae; + obj.fuselageMass = chosenMotor.mFus; + else + obj.length = chosenMotor.L; + obj.time = chosenMotor.t; + obj.thrust = chosenMotor.T; + obj.propellantMass = chosenMotor.m; + obj.structureMass = chosenMotor.mc; + end + if isempty(obj.cutoffTime) || obj.cutoffTime > obj.time(end) + obj.cutoffTime = obj.time(end); + end + end + end +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/classes/bays/Payload.m b/RoccarasoFlight/postprocessing/ASD/classes/bays/Payload.m new file mode 100644 index 0000000000000000000000000000000000000000..0f9f7edb3f5528648e3142d3cb85cfe508438e49 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/bays/Payload.m @@ -0,0 +1,55 @@ +classdef Payload < Bay + % Payload: Represents nose + payload configuration for a rocket. + % + % Constructor: + % - Payload: Creates an instance of the Payload class. + % Loaded config: rocketConfig.m > payload + % Loaded data: - + % Arguments: + % - mission: Mission, mission object + % - varIn: (optional) config source. Alternative to config.m + % file + + properties + noseLength double % [m] Nosecone length + noseType {mustBeMember(noseType, { ... + '', 'CONE', 'OGIVE', ... + 'POWER', 'HAACK', 'KARMAN', ... + 'MHAACK'})} = '' % [-] Nosecone shape + nosePower double % [-] Nosecone power type parameter + nosePMod double % [-] P coefficient for modified nosecone shapes + noseCMod double % [-] C coefficient for modified nosecone shapes + + position % [m] Absolute position, relative to nose base + length % [m] Total bay length + mass % [kg] Total bay mass + inertia % [kg*m^2] Total bay inertia (Body reference) + xCg % [m] Cg relative to bay upper side + end + + properties (Dependent) + inverseInertia double + end + + properties(Access = protected) + configName = 'rocketConfig.m' + variableName = 'payload' + mission Mission = Mission() + end + + methods + function obj = Payload(mission, varIn) + arguments(Input) + mission Mission = Mission() + varIn = [] + end + obj@Bay(mission, varIn); + end + end + + methods + function inverseInertia = get.inverseInertia(obj) + inverseInertia = inv(obj.inertiaMatrix); + end + end +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/classes/bays/Rear.m b/RoccarasoFlight/postprocessing/ASD/classes/bays/Rear.m new file mode 100644 index 0000000000000000000000000000000000000000..65d561d73eb9973703fe104184df2805f0981cd7 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/bays/Rear.m @@ -0,0 +1,53 @@ +classdef Rear < Bay +% Rear: Represents fincan + boat configuration for a rocket. +% +% Constructor: +% - Rear: Creates an instance of the Rear class. +% Loaded config: rocketConfig.m > rear +% Loaded data: - +% Arguments: +% - mission: Mission, mission object +% - varIn: (optional) config source. Alternative to config.m +% file + + properties + position % [m] Absolute position, relative to nose base + length % [m] Total bay length + mass % [kg] Total bay mass + inertia % [kg*m^2] Total bay inertia (Body reference) + xCg % [m] Cg relative to bay upper side + + boatType {mustBeMember(boatType, { ... + '', 'CONE', 'OGIVE'})} = '' % [] Boat type + boatLength double % [m] Boat length + boatFinalDiameter double % [m] Boat end diameter + + finsDeltaXRootChord double % [m] start of Chord 1 measured from beginning of bay + finsRootChord double % [m] attached chord length + finsFreeChord double % [m] free chord length + finsHeight double % [m] fin heigth + finsDeltaXFreeChord double % [m] start of Chord 2 measured from start of Chord 1 + nPanel double % [m] number of fins + finsLeadingEdgeRadius double % [deg] Leading edge radius at each span station + finsAxialDistance double % [m] distance between end of root chord and end of center body + finsSemiThickness double % [m] fin semi-thickness + finsMaxThicknessPosition double % [m] Fraction of chord from leading edge to max thickness + end + + properties(Access = protected) + configName = 'rocketConfig.m' + variableName = 'rear' + mission Mission = Mission() + end + + methods + function obj = Rear(mission, varIn) + arguments(Input) + mission Mission = Mission() + varIn = [] + end + obj@Bay(mission, varIn); + end + end +end + diff --git a/RoccarasoFlight/postprocessing/ASD/classes/bays/Recovery.m b/RoccarasoFlight/postprocessing/ASD/classes/bays/Recovery.m new file mode 100644 index 0000000000000000000000000000000000000000..b5d10de236fb9b2668a0c71493442c1cf19b903b --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/bays/Recovery.m @@ -0,0 +1,36 @@ +classdef Recovery < Bay +% Recovery: Represents recovery configuration for a rocket. +% +% Constructor: +% - Recovery: Creates an instance of the Recovery class. +% Loaded config: rocketConfig.m > recovery +% Loaded data: - +% Arguments: +% - mission: Mission, mission object +% - varIn: (optional) config source. Alternative to config.m +% file + + properties + position % [m] Absolute position, relative to nose base + length % [m] Total bay length + mass % [kg] Total bay mass + inertia % [kg*m^2] Total bay inertia (Body reference) + xCg % [m] Cg relative to bay upper side + end + + properties(Access = protected) + configName = 'rocketConfig.m' + variableName = 'recovery' + mission Mission = Mission() + end + + methods + function obj = Recovery(mission, varIn) + arguments(Input) + mission Mission = Mission() + varIn = [] + end + obj@Bay(mission, varIn); + end + end +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/classes/components/Environment.m b/RoccarasoFlight/postprocessing/ASD/classes/components/Environment.m new file mode 100644 index 0000000000000000000000000000000000000000..e9b47cc6aa304805dea5d0c326ec95379bfcbc75 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/components/Environment.m @@ -0,0 +1,94 @@ +classdef Environment < Component +% Environment: Represents launch site dependent variables. +% +% Constructor: +% - Environment: Creates an instance of the Environment class. +% Loaded config: environmentConfig.m +% Loaded data: - +% Arguments: +% - mission: Mission, mission object +% - motor: Motor, used to compute pin distance +% - varIn: (optional) config source. Alternative to config.m +% file + + properties + lat0 double % [deg] Launchpad latitude + lon0 double % [deg] Launchpad longitude + z0 double % [m] Launchpad Altitude + omega double % [deg] Launchpad Elevation + phi double % [deg] Launchpad Azimuth + pin1Length double % [m] Distance from the upper pin to the upper tank cap + pin2Length double % [m] Distance from the lower pin to the lower tank cap + rampLength double % [m] Total launchpad length + + temperature double = 288.15 % [K] Ground temperature + pressure double % [Pa] Ground pressure + rho double % [Kg/m^3] Ground air density + gamma double = 1.4 % [-] Gas constant + end + + properties(SetAccess = private) + g0 double % [-] Gravity costant at launch latitude and altitude + local double % [-] Vector conatining inputs for atmosphereData + pinDistance double % [m] Distance of the upper pin from the rail base (upper pin-boat + boat-rail base) + effectiveRampLength double % [m] Total launchpad length + effectiveRampAltitude double % [m] Projection of effectiveRampLength on z axis + end + + properties(Access = protected) + configName = 'environmentConfig.m' + variableName = '' + mission Mission + motor Motor + end + + methods + function obj = Environment(mission, motor, varIn) + arguments(Input) + mission Mission = Mission() + motor Motor = Motor() + varIn = [] + end + obj@Component(mission, varIn); + obj.motor = motor; + + %% Update angles + obj.omega = deg2rad(obj.omega); + obj.phi = deg2rad(obj.phi); + + obj.updateAll(); + end + end + + % Updaters + methods + function updateAll(obj) + obj.updateG0; + obj.updatePinDistance; + obj.updateRamp; + obj.updateLocal; + end + + function obj = updateG0(obj) + obj.g0 = gravitywgs84(obj.z0, obj.lat0); + end + + function obj = updatePinDistance(obj) + obj.pinDistance = obj.pin1Length + obj.pin2Length ... + + obj.motor.tankLength; + end + + function obj = updateRamp(obj) + obj.effectiveRampLength = obj.rampLength - obj.pinDistance; + obj.effectiveRampAltitude = obj.effectiveRampLength*sin(obj.omega); + end + + function obj = updateLocal(obj) + if isempty(obj.temperature), obj.temperature = 288.15; end + if isempty(obj.gamma), obj.gamma = 1.4; end + + obj.local = [obj.z0, obj.temperature, ... % vector containing inputs for atmosphereData + obj.pressure, obj.rho]; + end + end +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/classes/components/Para.m b/RoccarasoFlight/postprocessing/ASD/classes/components/Para.m new file mode 100644 index 0000000000000000000000000000000000000000..ef74338bcccc5f29d0943eaa1efb635286f06495 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/components/Para.m @@ -0,0 +1,52 @@ +classdef Para < matlab.mixin.Heterogeneous & Component + % Parachute: Represents a parachute component. + % + % Constructor: + % - Parachute: Creates an instance of the Parachute class. + % Loaded config: rocketConfig.m > parachute + % Loaded data: - + % Arguments: + % - mission: Mission, mission object + % - varIn: (optional) config source. Alternative to config.m + % file + % + % WARNING: To preserve object type when creating an array, a third + % dimension is added to force the creation of an heterogeneous array. + % the third dimension can be trimmed away. + % Do NOT consider the third dimension as valid + + properties(Access = protected) + configName = 'paraConfig.m' + variableName = 'para' + mission Mission + end + + methods + function obj = Para(mission, varIn) + arguments(Input) + mission Mission = Mission() + varIn = [] + end + obj@Component(mission, varIn, "elementWise", false); + + %% Forcing creation of heterogeneous array + if any(size(obj) > 1) + switch class(obj) + % In case of homogeneous type, forces conversion to + % heterogeneous Para array + case 'Parachute' + obj(1, 1, 2) = Parafoil(); + case 'Parafoil' + obj(1, 1, 2) = Parachute(); + otherwise + end + end + end + end + + methods (Static,Sealed,Access=protected) + function default = getDefaultScalarElement + default = Parachute; + end + end +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/classes/components/Parachute.m b/RoccarasoFlight/postprocessing/ASD/classes/components/Parachute.m new file mode 100644 index 0000000000000000000000000000000000000000..9144fd0edf5ea2ffd57e9817ccbd55937d475ca3 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/components/Parachute.m @@ -0,0 +1,43 @@ +classdef Parachute < Para +% Parachute: Represents a parachute component. +% +% Constructor: +% - Parachute: Creates an instance of the Parachute class. +% Loaded config: rocketConfig.m > parachute +% Loaded data: - +% Arguments: +% - mission: Mission, mission object +% - varIn: (optional) config source. Alternative to config.m +% file + + properties + name = '' + surface double % [m^2] Surface + mass double % [kg] Parachute Mass + openingTime double % [s] Parachute opening time (delay + time to get fully open) + finalAltitude double % [m] Final altitude of the parachute + chordLength double % [m] Shock Chord Length + chordK double % [N/m^2] Shock Chord Elastic Constant + chordC double % [Ns/m] Shock Chord Dynamic Coefficient + expulsionSpeed double + cx double % [/] Parachute Longitudinal Drag Coefficient + cd double % [/] Parachute Drag Coefficient + cl double % [/] Parachute Lift Coefficient + m double % [m^2/s] Coefficient of the surface vs. time opening model + nf double % [/] Adimensional Opening Time + forceCoefficient double % [-] Empirical coefficient to obtain correct peak force at deployment + end + + methods + function obj = Parachute(mission, varIn) + arguments(Input) + mission Mission = Mission() + varIn = [] + end + + obj@Para(mission, varIn); + obj = obj(:, :, 1); + end + end +end + diff --git a/RoccarasoFlight/postprocessing/ASD/classes/components/Parafoil.m b/RoccarasoFlight/postprocessing/ASD/classes/components/Parafoil.m new file mode 100644 index 0000000000000000000000000000000000000000..d02cc7bd2089750af7c23c2995b3b427843378de --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/components/Parafoil.m @@ -0,0 +1,77 @@ +classdef Parafoil < Para +% Parachute: Represents a parachute component. +% +% Constructor: +% - Parachute: Creates an instance of the Parachute class. +% Loaded config: rocketConfig.m > parachute +% Loaded data: - +% Arguments: +% - mission: Mission, mission object +% - varIn: (optional) config source. Alternative to config.m +% file + + properties + name = '' + surface double % [m^2] Surface + deltaSMax double % [-] Aerodynamic control coefficient - symmetric, max value + + uMax + uMin + identification + deltaATau + maxSpeed + + mass double % [kg] Parafoil Mass + openingTime double % [s] Parafoil opening time (delay + time to get fully open) + + inertia % [kg*m^2] 3x3 Inertia Matrix (the diagonal is made of the inertia vector values) + finalAltitude double % [m] Final altitude of the parafoil + semiWingSpan double % [m] Parafoil semiwingspan + MAC double % [m] Parafoil mean aero chord + + cd0 double % [-] aerodynamic control coefficient - asymetric, CD0 + cdAlpha double % [-] aerodynamic control coefficient - asymetric, CDAlpha + cdSurface double % [-] aerodynamic control coefficient - asymetric, CDSurface + + cl0 double % [-] aerodynamic control coefficient - asymetric, CL0 + clAlpha double % [-] aerodynamic control coefficient - asymetric, CLAlpha + clSurface double % [-] aerodynamic control coefficient - asymetric, CLSurface + + cLP double % [-] aerodynamic control coefficient - asymetric, Clp + cLPhi double % [-] aerodynamic control coefficient - asymetric, ClPhi + cLSurface + + cM0 double % [-] aerodynamic control coefficient - asymetric, Cm0 + cMAlpha double % [-] aerodynamic control coefficient - asymetric, CmAlpha + cMQ double % [-] aerodynamic control coefficient - asymetric, Cmq + + cNR double % [-] aerodynamic control coefficient - asymetric, Cnr + cNSurface double % [-] aerodynamic control coefficient - asymetric, CnDeltaA + end + + properties(SetAccess = private) + inverseInertia + end + + methods + function obj = Parafoil(mission, varIn) + arguments(Input) + mission Mission = Mission() + varIn = [] + end + + obj@Para(mission, varIn); + obj = obj(:, :, 1); + end + + function set.inertia(obj, inertia) + obj.inertia = inertia; + obj.updateAll(); + end + + function updateAll(obj) + if isempty(obj.inertia), return; end + obj.inverseInertia = obj.inertia\eye(3); + end + end +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/classes/components/Pitot.m b/RoccarasoFlight/postprocessing/ASD/classes/components/Pitot.m new file mode 100644 index 0000000000000000000000000000000000000000..1d8cd4885a8e90c1f0fcb17f681edfdbf77f5118 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/components/Pitot.m @@ -0,0 +1,37 @@ +classdef Pitot < Component +% Pitot: Represents a pitot tube. +% +% Constructor: +% - Pitot: Creates an instance of the Pitot class. +% Loaded config: rocketConfig.m > pitot +% Loaded data: - +% Arguments: +% - mission: Mission, mission object +% - varIn: (optional) config source. Alternative to config.m +% file + + properties + length double % [m] Pitot tube length + diameter double % [m] Pitot tube diameter + initialConeLength double % [m] Pitot initial conic section length + finalConeLength double % [m] Pitot final conic section length + initialConeDiameter double % [m] Pitot initial conic section diameter + finalConeDiameter double % [m] Pitot final conic section diameter + end + + properties(Access = protected) + configName = 'rocketConfig.m' + variableName = 'pitot' + mission Mission = Mission() + end + + methods + function obj = Pitot(mission, varIn) + arguments(Input) + mission Mission = Mission() + varIn = [] + end + obj@Component(mission, varIn); + end + end +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/classes/components/WindCustom.m b/RoccarasoFlight/postprocessing/ASD/classes/components/WindCustom.m new file mode 100644 index 0000000000000000000000000000000000000000..99d26ba500d01891a7ca153c383c10281282ed64 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/components/WindCustom.m @@ -0,0 +1,92 @@ +classdef WindCustom < Component +% WindCustom: Represents Skyward's custom wind model. +% +% Constructor: +% - WindCustom: Creates an instance of the WindCustom class. +% Loaded config: windConfig.m > windCustom +% Loaded data: - +% Arguments: +% - mission: Mission, mission object +% - varIn: (optional) config source. Alternative to config.m +% file + properties + altitudes (1, :) + + magnitudeDistribution (1, :) ... + {mustBeMember(magnitudeDistribution, {'u', 'g'})} = 'u' % [-] Magnitude distrubution: uniform, gaussian + azimuthDistribution (1, :) ... + {mustBeMember(azimuthDistribution, {'u', 'g'})} = 'u' % [-] Azimuth distrubution: uniform, gaussian + magnitudeParameters (2, :) % [m/s] Distribution parameters: [min; max], [mu; sigma] + azimuthParameters (2, :) % [deg] Distribution parameters: [min; max], [mu; sigma] + end + + properties(SetAccess = private) + azimuth + magnitude + end + + properties(Access = protected) + configName = 'windConfig.m' + variableName = 'windCustom' + mission Mission + end + + methods + function obj = WindCustom(mission, varIn) + arguments(Input) + mission Mission = Mission() + varIn = [] + end + obj@Component(mission, varIn); + obj.updateAll(); + end + + function [uw, vw, ww] = getVels(obj, z) + s = length(obj.altitudes); % reference vectors length + if s==1 + uw = round( - obj.magnitude * cos(obj.azimuth), 6); + vw = round( - obj.magnitude * sin(obj.azimuth), 6); + ww = 0; + else + mag = interpLinear(obj.altitudes, obj.magnitude, z, true); + az = interpLinear(obj.altitudes, obj.azimuth, z, true); + + uw = round( - mag * cos(az), 6); + vw = round( - mag * sin(az), 6); + ww = 0; + end + end + + function updateAll(obj) + s = length(obj.altitudes); + magVector = nan(1,s); + azVector = nan(1,s); + isUniformMag = strcmp(obj.magnitudeDistribution, "u"); + isUniformAz = strcmp(obj.azimuthDistribution, "u"); + + uniformDist = @(val1, val2) rand(1)*(val2 - val1) + val1; + gaussianDist = @(mean, sigma) normrnd(mean, sigma, 1); + + for i = 1:s + if isUniformMag(i) + magVector(i) = uniformDist(obj.magnitudeParameters(1,i), ... + obj.magnitudeParameters(2,i)); + else + magVector(i) = gaussianDist(obj.magnitudeParameters(1,i), ... + obj.magnitudeParameters(2,i)); + end + + if isUniformAz(i) + azVector(i) = uniformDist(obj.azimuthParameters(1,i), ... + obj.azimuthParameters(2,i)); + else + azVector(i) = gaussianDist(obj.azimuthParameters(1,i), ... + obj.azimuthParameters(2,i)); + end + end + + obj.magnitude = magVector; + obj.azimuth = azVector; + end + end +end diff --git a/RoccarasoFlight/postprocessing/ASD/classes/components/WindMatlab.m b/RoccarasoFlight/postprocessing/ASD/classes/components/WindMatlab.m new file mode 100644 index 0000000000000000000000000000000000000000..b47c45d6b4cf8f56fb1c72a65b41beb5796d75cc --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/classes/components/WindMatlab.m @@ -0,0 +1,62 @@ +classdef WindMatlab < Component +% WindMatlab: Represents matlab wind variables. +% +% Constructor: +% - WindMatlab: Creates an instance of the WindMatlab class. +% Loaded config: windConfig.m > windMatlab +% Loaded data: - +% Arguments: +% - mission: Mission, mission object +% - environment: Environment, used to get launch coordinates +% - varIn: (optional) config source. Alternative to config.m +% file + + properties + DayMin % [d] Minimum Day of the launch + DayMax % [d] Maximum Day of the launch + HourMin % [h] Minimum Hour of the day + HourMax % [h] Maximum Hour of the day + ww % [m/s] Vertical wind speed + end + + properties(Access = protected) + configName = 'windConfig.m' + variableName = 'windMatlab' + mission Mission + environment Environment + end + + methods + function obj = WindMatlab(mission, environment, varIn) + arguments(Input) + mission Mission = Mission() + environment Environment = Environment() + varIn = [] + end + if nargin > 0 && nargin < 2, error('MATLAB:narginchk:notEnoughInputs', ... + 'Not enough input arguments.'); end + obj@Component(mission, varIn); + obj.environment = environment; + end + + function [uw, vw, ww] = getVels(obj, z, t, Hour, Day) + h = z + obj.environment.z0; + if h < 0 + h = 0; + end + + if nargin == 3 + if obj.HourMin == obj.HourMax + Day = obj.DayMin; + Hour = obj.HourMin; + end + end + + Seconds = Hour*3600; + % horizontal wind + [uw,vw] = atmoshwm(obj.environment.lat0, obj.environment.lon0,h,'day',Day,... + 'seconds',Seconds+t,'model','quiet','version','14'); % NED reference + ww = obj.ww; + end + end +end diff --git a/RoccarasoFlight/postprocessing/ASD/data/geminiFlightRoccaraso.mat b/RoccarasoFlight/postprocessing/ASD/data/geminiFlightRoccaraso.mat new file mode 100644 index 0000000000000000000000000000000000000000..6393467cdbaf599b623f4e7928b82799ebb2e402 Binary files /dev/null and b/RoccarasoFlight/postprocessing/ASD/data/geminiFlightRoccaraso.mat differ diff --git a/RoccarasoFlight/postprocessing/ASD/data/geminiRocketRoccaraso.mat b/RoccarasoFlight/postprocessing/ASD/data/geminiRocketRoccaraso.mat new file mode 100644 index 0000000000000000000000000000000000000000..9df327cd9f3181733f097e7f07c7c2b4243a2a8c Binary files /dev/null and b/RoccarasoFlight/postprocessing/ASD/data/geminiRocketRoccaraso.mat differ diff --git a/RoccarasoFlight/postprocessing/ASD/plots.m b/RoccarasoFlight/postprocessing/ASD/plots.m new file mode 100644 index 0000000000000000000000000000000000000000..79288cb13410a1548c8807644c4d605792b8f449 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/plots.m @@ -0,0 +1,42 @@ +% plots - Plots the progress in time of the measured and simulated +% internal loads during parachute's deployment. +% +% REVISIONS: +% - 0 16-07-2024, Release, Filippo Cataldo + +%% PLOT INTERNAL LOADS +for i = 1:nBays + + iMeasuredLoads = measuredLoads(:,i,:); + iSimulatedLoads = simulatedLoads(:,i,:); + + figure() + hold on; + grid on; + plot(time, iMeasuredLoads(1,:), 'Color', '#0072BD', 'LineWidth', 2); + plot(simulationTime, iSimulatedLoads(1,:), 'Color',... + '#0072BD', 'LineWidth', 2, 'LineStyle', '--'); + plot(time, iMeasuredLoads(2,:), 'Color', "#D95319", 'LineWidth', 2); + plot(simulationTime, vecnorm(iSimulatedLoads(2,:), 2, 1), 'Color',... + "#D95319", 'LineWidth', 2, 'LineStyle', '--'); + plot(time, iMeasuredLoads(3,:), 'Color', "#EDB120", 'LineWidth', 2); + plot(simulationTime, vecnorm(iSimulatedLoads(3,:), 2, 1), 'Color',... + "#EDB120", 'LineWidth', 2, 'LineStyle', '--'); + + if i > 1 + title(join(['Internal loads' newline settings.parachuteName... + aeroDistribution 'pressure distribution |' namesBays(i) ' - '... + namesBays(i-1)]), 'FontSize', 12); + else + title(join(['Internal loads' newline settings.parachuteName... + aeroDistribution 'pressure distribution |' namesBays(i)... + ' - ']), 'FontSize', 12); + end + + xlim('Tight'); + xlabel('Time [s]', 'FontSize', 12); + ylabel('Force [N], Moment [Nm]', 'FontSize', 12); + legend('Measured N', 'Simulated N', 'Measured S',... + 'Simulated S', 'Measured M', 'Simulated M', 'FontSize', 10); + +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/postProcessDrogue.m b/RoccarasoFlight/postprocessing/ASD/postProcessDrogue.m new file mode 100644 index 0000000000000000000000000000000000000000..2f251648307442a919cac3e739a97197d63144bf --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/postProcessDrogue.m @@ -0,0 +1,222 @@ +% postProcessMain - Processes flight data of Gemini's launch at +% Roccaraso 2023 and gets estimates of DROGUE parachute's +% loads at deployment. Then compares them with simulated +% loads. +% NAS time is used - with lower sampling frequency. +% No specification of reference frame means 2D inertial +% frame. Body frame can be 2D or 3D depending on the +% context. Anyway, frames are only used to rotate +% vectors: velocities or acceleration are always with +% respect to an inertial frame. +% +% CALLED FUNCTIONS: eqResidual.m, aeroResultant.m, internalLoads.m +% +% REVISIONS: +% - 0 16-07-2024, Release, Filippo Cataldo + +clear +close all +clc + +path = genpath(pwd); +addpath(path); +% Necessary classes are needed in the path to correctly read Rocket object +load("data\geminiRocketRoccaraso.mat"); +load("data\geminiFlightRoccaraso.mat"); + +%% Set parameters for the post process +% Pressure distribution type, only 'uniform' supported for now +aeroDistribution = 'uniform'; + +indexOpen = 975; % Drogue opening esimtated index in quaternions time +gNED = [0; 0; 9.81]; % Gravity acceleration in NED [m/s^2] +xImu = 0.1854; % [m] x position of IMU from top of ELC bay + +% Options for non-linear system +guess = [200*ones(2,1); ones(2,1)]; % Initial guess +nonlinOptions = optimoptions('fsolve'); +nonlinOptions.Algorithm = 'levenberg-marquardt'; + +%% Set parameters for simulation +% "settings" struct needed to run the simulation +settings.nameFirstBay = 'recovery'; +settings.parachuteName = 'DROGUE chute'; +settings.odeOptions = odeset('RelTol', 1e-8, 'AbsTol', 1e-8); +settings.g = [0; -9.81]; % Gravity acceleration vector [m/s^2] +settings.time = [NAS.time(indexOpen) NAS.time(indexOpen) + 2.5]; % [s] + +%% Get important rocket's parameters +bays = values(rocket.bays); +namesBays = keys(rocket.bays); +indexFirstBay = find(namesBays == settings.nameFirstBay); +bays = bays(indexFirstBay:end); +namesBays = namesBays(indexFirstBay:end); +nBays = length(bays); + +% Position of each bay from the base of the nosecone +absolutePositions = values(rocket.absolutePositions); +absolutePositions = absolutePositions(indexFirstBay:end); % [m] + +diameter = rocket.diameter; % [m] +xCg = rocket.xCg(end); % [m] +armImu = xCg - (rocket.absolutePositions('electronics') +... + xImu); % Accelerometer position from xCg [m] + +%% Get flight data +% Time [s] +nasTime = NAS.time'; % Time span of NAS +imuTime = IMU.accTime'; % Time span of IMU +time = nasTime; % NAS time is used - with lower sampling frequency +nTime = length(time); +dt = time(2) - time(1); +indexEnd = indexOpen + ceil((settings.time(end) - settings.time(1))/dt); + +% Quaternions [-] +quat = zeros(4,1,nTime); +quat(:,:,:) = [interp1(nasTime, NAS.states(:,10)', time) + interp1(nasTime, NAS.states(:,7)', time) + interp1(nasTime, NAS.states(:,8)', time) + interp1(nasTime, NAS.states(:,9)', time)]; +quat = squeeze(quat); +Rbn = quat2rotm(quat'); % 3D rotation from body to NED +Rnb = pagetranspose(Rbn); % 3D rotation from NED to body + +% Position [m] +posNED = zeros(3,nTime); +posNED(:,:,:) = [interp1(nasTime, NAS.states(:,1)', time) + interp1(nasTime, NAS.states(:,2)', time) + interp1(nasTime, NAS.states(:,3)', time)]; + +% Velocity [m/s] +velNED = zeros(3,1,nTime); +velNED(:,:,:) = [interp1(nasTime, NAS.states(:,4)', time) + interp1(nasTime, NAS.states(:,5)', time) + interp1(nasTime, NAS.states(:,6)', time)]; +velBody = pagemtimes(Rnb, velNED); +velBody = squeeze(velBody); +velNED = squeeze(velNED); + +% Acceleration [m/s^2] +accImu = IMU.accelerometerMeasures'; +accImu = [interp1(imuTime, accImu(1,:), time) + interp1(imuTime, accImu(2,:), time) + interp1(imuTime, accImu(3,:), time)]; + +% Angular velocity [rad/s] +omegaBody = IMU.gyroMeasures'; +omegaBody = [interp1(imuTime, omegaBody(1,:), time) + interp1(imuTime, omegaBody(2,:), time) + interp1(imuTime, omegaBody(3,:), time)]; + +% Angular acceleration [rad/s^2] +omegaTemp = zeros(3,1,nTime); +omegaTemp(:,:,:) = [omegaBody(1,:); omegaBody(2,:); omegaBody(3,:)]; +omegaNED = pagemtimes(Rbn, omegaTemp); +omegaNED = squeeze(omegaNED); % Transform to NED to differentiate + +omegaDotNED = [zeros(3,1) diff(omegaNED, [], 2)/dt]; +omegaDotTemp = zeros(3,1,nTime); +omegaDotTemp(:,:,:) = [omegaDotNED(1,:); omegaDotNED(2,:); omegaDotNED(3,:)]; +omegaDotBody = pagemtimes(Rnb, omegaDotTemp); +omegaDotBody = squeeze(omegaDotBody); + +% Use time of parachute's opening +time = time(indexOpen:indexEnd); +nTime = length(time); +posNED = posNED(:,indexOpen:indexEnd); +velBody = velBody(:,indexOpen:indexEnd); +velNED = velNED(:,indexOpen:indexEnd); +accImu = accImu(:,indexOpen:indexEnd); +omegaNED = omegaNED(:,indexOpen:indexEnd); +omegaBody = omegaBody(:,indexOpen:indexEnd); +omegaDotNED = omegaDotNED(:,indexOpen:indexEnd); +omegaDotBody = omegaDotBody(:,indexOpen:indexEnd); +Rbn = Rbn(:,:,indexOpen:indexEnd); +Rnb = Rnb(:,:,indexOpen:indexEnd); + +%% Compute measured internal loads +dStateBody = zeros(12, nTime); +accBody = zeros(3, nTime); % CG acceleration [m/s^2] +measuredForces = zeros(3, nBays + 1, nTime); +measuredMoments = zeros(3, nBays + 1, nTime); +measuredLoads = zeros(3, nBays + 1, nTime); + +for i = 1:nTime + + % CHECK FOR SIGN OF RADIAL DISTANCE; SHOULD BE POSITIVE IN CENTRIFUGAL + accBody(:,i) = accImu(:,i) - cross(omegaDotBody(:,i),... + [armImu; 0; diameter/2]) - cross(omegaBody(:,i),... + cross(omegaBody(:,i), [armImu; 0; diameter/2])); + + dStateBody(:,i) = [velBody(:,i) + accBody(:,i) + omegaBody(:,i) + omegaDotBody(:,i)]; + + iRnb = Rnb(:,:,i); + + % Solve non-linear system to impose global equilibrium through + % aerodynamic forces + sol = fsolve(@(var) eqResidual(var, dStateBody(:,i), iRnb*gNED,... + aeroDistribution, rocket, settings),... + guess, nonlinOptions); + aeroAmplitude = [0; sol(3:4)]; + + for j = nBays:-1:1 + + % Extract parameters of current bay + absolutePositionBay = absolutePositions(j); + + if j < nBays + lengthBay = absolutePositions(j+1) - absolutePositionBay; + else + lengthBay = rocket.length - absolutePositionBay +... + absolutePositions(1); + end + + xCgBay = bays(j).xCg(end); + armBay = xCg - (absolutePositionBay + xCgBay); + mBay = bays(j).mass(end); + inertiaBay = diag(bays(j).inertia); % Bay's inertia tensor [kg*m^2] + + accBayBody = accBody(:,i) + cross(omegaDotBody(:,i),... + [armBay; 0; 0]) + cross(omegaBody(:,i),... + cross(omegaBody(:,i), [armBay; 0; 0])); + + [aeroForce, aeroMoment] = aeroResultant(lengthBay, xCgBay,... + aeroAmplitude, aeroDistribution); + + % Compute loads + measuredForces(:,j,i) = measuredForces(:,j+1,i) +... + mBay*(accBayBody - iRnb*gNED) - aeroForce; + measuredMoments(:,j,i) = measuredMoments(:,j+1,i) +... + cross(-[lengthBay - xCgBay; 0; 0], measuredForces(:,j+1,i)) -... + cross([xCgBay; 0; 0], measuredForces(:,j,i)) +... + inertiaBay*omegaDotBody(:,i) +... + cross(omegaBody(:,i), inertiaBay*omegaBody(:,i)) - aeroMoment; + + measuredLoads(1,j,i) = measuredForces(1,j,i); + measuredLoads(2,j,i) = norm(measuredForces(2:3,j,i)); + measuredLoads(3,j,i) = norm(measuredMoments(2:3,j,i)); + + end + +end + +% Free interface is not of interest +measuredForces = measuredForces(:,1:nBays,:); +measuredMoments = measuredMoments(:,1:nBays,:); +measuredLoads = measuredLoads(:,1:nBays,:); + +%% Simulate internal loads +alt0 = -posNED(3,1); +velX0 = norm(velNED(1:2,1)); +velY0 = -velNED(3,1); +theta0 = -asin(Rnb(1,3,1)); +AoA0 = theta0 - atan2(velY0, velX0); +state0 = [alt0; velX0; velY0; AoA0]; +[~, simulatedLoads, simulationTime] =... + internalLoads(state0, rocket, settings); + +%% Plot results +plots; \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/postProcessMain.m b/RoccarasoFlight/postprocessing/ASD/postProcessMain.m new file mode 100644 index 0000000000000000000000000000000000000000..ab5c93302936bbca8a800b0d30e7b1a9043d96ec --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/postProcessMain.m @@ -0,0 +1,222 @@ +% postProcessMain - Processes flight data of Gemini's launch at +% Roccaraso 2023 and gets estimates of MAIN parachute's +% loads at deployment. Then compares them with simulated +% loads. +% NAS time is used - with lower sampling frequency. +% No specification of reference frame means 2D inertial +% frame. Body frame can be 2D or 3D depending on the +% context. Anyway, frames are only used to rotate +% vectors: velocities or acceleration are always with +% respect to an inertial frame. +% +% CALLED FUNCTIONS: eqResidual.m, aeroResultant.m, internalLoads.m +% +% REVISIONS: +% - 0 16-07-2024, Release, Filippo Cataldo + +clear +close all +clc + +path = genpath(pwd); +addpath(path); +% Necessary classes are needed in the path to correctly read Rocket object +load("data\geminiRocketRoccaraso.mat"); +load("data\geminiFlightRoccaraso.mat"); + +%% Set parameters for the post process +% Pressure distribution type, only 'uniform' supported for now +aeroDistribution = 'uniform'; + +indexOpen = 3585; % Main opening esimtated index in quaternions time +gNED = [0; 0; 9.81]; % Gravity acceleration in NED [m/s^2] +xImu = 0.1854; % [m] x position of IMU from top of ELC bay + +% Options for non-linear system +guess = [200*ones(2,1); ones(2,1)]; % Initial guess +nonlinOptions = optimoptions('fsolve'); +nonlinOptions.Algorithm = 'levenberg-marquardt'; + +%% Set parameters for simulation +% "settings" struct needed to run the simulation +settings.nameFirstBay = 'recovery'; +settings.parachuteName = 'MAIN chute'; +settings.odeOptions = odeset('RelTol', 1e-8, 'AbsTol', 1e-8); +settings.g = [0; -9.81]; % Gravity acceleration vector [m/s^2] +settings.time = [NAS.time(indexOpen) NAS.time(indexOpen) + 2.5]; % [s] + +%% Get important rocket's parameters +bays = values(rocket.bays); +namesBays = keys(rocket.bays); +indexFirstBay = find(namesBays == settings.nameFirstBay); +bays = bays(indexFirstBay:end); +namesBays = namesBays(indexFirstBay:end); +nBays = length(bays); + +% Position of each bay from the base of the nosecone +absolutePositions = values(rocket.absolutePositions); +absolutePositions = absolutePositions(indexFirstBay:end); % [m] + +diameter = rocket.diameter; % [m] +xCg = rocket.xCg(end); % [m] +armImu = xCg - (rocket.absolutePositions('electronics') +... + xImu); % Accelerometer position from xCg [m] + +%% Get flight data +% Time [s] +nasTime = NAS.time'; % Time span of NAS +imuTime = IMU.accTime'; % Time span of IMU +time = nasTime; % NAS time is used - with lower sampling frequency +nTime = length(time); +dt = time(2) - time(1); +indexEnd = indexOpen + ceil((settings.time(end) - settings.time(1))/dt); + +% Quaternions [-] +quat = zeros(4,1,nTime); +quat(:,:,:) = [interp1(nasTime, NAS.states(:,10)', time) + interp1(nasTime, NAS.states(:,7)', time) + interp1(nasTime, NAS.states(:,8)', time) + interp1(nasTime, NAS.states(:,9)', time)]; +quat = squeeze(quat); +Rbn = quat2rotm(quat'); % 3D rotation from body to NED +Rnb = pagetranspose(Rbn); % 3D rotation from NED to body + +% Position [m] +posNED = zeros(3,nTime); +posNED(:,:,:) = [interp1(nasTime, NAS.states(:,1)', time) + interp1(nasTime, NAS.states(:,2)', time) + interp1(nasTime, NAS.states(:,3)', time)]; + +% Velocity [m/s] +velNED = zeros(3,1,nTime); +velNED(:,:,:) = [interp1(nasTime, NAS.states(:,4)', time) + interp1(nasTime, NAS.states(:,5)', time) + interp1(nasTime, NAS.states(:,6)', time)]; +velBody = pagemtimes(Rnb, velNED); +velBody = squeeze(velBody); +velNED = squeeze(velNED); + +% Acceleration [m/s^2] +accImu = IMU.accelerometerMeasures'; +accImu = [interp1(imuTime, accImu(1,:), time) + interp1(imuTime, accImu(2,:), time) + interp1(imuTime, accImu(3,:), time)]; + +% Angular velocity [rad/s] +omegaBody = IMU.gyroMeasures'; +omegaBody = [interp1(imuTime, omegaBody(1,:), time) + interp1(imuTime, omegaBody(2,:), time) + interp1(imuTime, omegaBody(3,:), time)]; + +% Angular acceleration [rad/s^2] +omegaTemp = zeros(3,1,nTime); +omegaTemp(:,:,:) = [omegaBody(1,:); omegaBody(2,:); omegaBody(3,:)]; +omegaNED = pagemtimes(Rbn, omegaTemp); +omegaNED = squeeze(omegaNED); % Transform to NED to differentiate + +omegaDotNED = [zeros(3,1) diff(omegaNED, [], 2)/dt]; +omegaDotTemp = zeros(3,1,nTime); +omegaDotTemp(:,:,:) = [omegaDotNED(1,:); omegaDotNED(2,:); omegaDotNED(3,:)]; +omegaDotBody = pagemtimes(Rnb, omegaDotTemp); +omegaDotBody = squeeze(omegaDotBody); + +% Use time of parachute's opening +time = time(indexOpen:indexEnd); +nTime = length(time); +posNED = posNED(:,indexOpen:indexEnd); +velBody = velBody(:,indexOpen:indexEnd); +velNED = velNED(:,indexOpen:indexEnd); +accImu = accImu(:,indexOpen:indexEnd); +omegaNED = omegaNED(:,indexOpen:indexEnd); +omegaBody = omegaBody(:,indexOpen:indexEnd); +omegaDotNED = omegaDotNED(:,indexOpen:indexEnd); +omegaDotBody = omegaDotBody(:,indexOpen:indexEnd); +Rbn = Rbn(:,:,indexOpen:indexEnd); +Rnb = Rnb(:,:,indexOpen:indexEnd); + +%% Compute measured internal loads +dStateBody = zeros(12, nTime); +accBody = zeros(3, nTime); % CG acceleration [m/s^2] +measuredForces = zeros(3, nBays + 1, nTime); +measuredMoments = zeros(3, nBays + 1, nTime); +measuredLoads = zeros(3, nBays + 1, nTime); + +for i = 1:nTime + + % CHECK FOR SIGN OF RADIAL DISTANCE; SHOULD BE POSITIVE IN CENTRIFUGAL + accBody(:,i) = accImu(:,i) - cross(omegaDotBody(:,i),... + [armImu; 0; diameter/2]) - cross(omegaBody(:,i),... + cross(omegaBody(:,i), [armImu; 0; diameter/2])); + + dStateBody(:,i) = [velBody(:,i) + accBody(:,i) + omegaBody(:,i) + omegaDotBody(:,i)]; + + iRnb = Rnb(:,:,i); + + % Solve non-linear system to impose global equilibrium through + % aerodynamic forces + sol = fsolve(@(var) eqResidual(var, dStateBody(:,i), iRnb*gNED,... + aeroDistribution, rocket, settings),... + guess, nonlinOptions); + aeroAmplitude = [0; sol(3:4)]; + + for j = nBays:-1:1 + + % Extract parameters of current bay + absolutePositionBay = absolutePositions(j); + + if j < nBays + lengthBay = absolutePositions(j+1) - absolutePositionBay; + else + lengthBay = rocket.length - absolutePositionBay +... + absolutePositions(1); + end + + xCgBay = bays(j).xCg(end); + armBay = xCg - (absolutePositionBay + xCgBay); + mBay = bays(j).mass(end); + inertiaBay = diag(bays(j).inertia); % Bay's inertia tensor [kg*m^2] + + accBayBody = accBody(:,i) + cross(omegaDotBody(:,i),... + [armBay; 0; 0]) + cross(omegaBody(:,i),... + cross(omegaBody(:,i), [armBay; 0; 0])); + + [aeroForce, aeroMoment] = aeroResultant(lengthBay, xCgBay,... + aeroAmplitude, aeroDistribution); + + % Compute loads + measuredForces(:,j,i) = measuredForces(:,j+1,i) +... + mBay*(accBayBody - iRnb*gNED) - aeroForce; + measuredMoments(:,j,i) = measuredMoments(:,j+1,i) +... + cross(-[lengthBay - xCgBay; 0; 0], measuredForces(:,j+1,i)) -... + cross([xCgBay; 0; 0], measuredForces(:,j,i)) +... + inertiaBay*omegaDotBody(:,i) +... + cross(omegaBody(:,i), inertiaBay*omegaBody(:,i)) - aeroMoment; + + measuredLoads(1,j,i) = measuredForces(1,j,i); + measuredLoads(2,j,i) = norm(measuredForces(2:3,j,i)); + measuredLoads(3,j,i) = norm(measuredMoments(2:3,j,i)); + + end + +end + +% Free interface is not of interest +measuredForces = measuredForces(:,1:nBays,:); +measuredMoments = measuredMoments(:,1:nBays,:); +measuredLoads = measuredLoads(:,1:nBays,:); + +%% Simulate internal loads +alt0 = -posNED(3,1); +velX0 = norm(velNED(1:2,1)); +velY0 = -velNED(3,1); +theta0 = -asin(Rnb(1,3,1)); +AoA0 = theta0 - atan2(velY0, velX0); +state0 = [alt0; velX0; velY0; AoA0]; +[~, simulatedLoads, simulationTime] =... + internalLoads(state0, rocket, settings); + +%% Plot results +plots; \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/src/aeroResultant.m b/RoccarasoFlight/postprocessing/ASD/src/aeroResultant.m new file mode 100644 index 0000000000000000000000000000000000000000..48e38e47a8597f5fcd6e673953838d244a7845ac --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/src/aeroResultant.m @@ -0,0 +1,34 @@ +function [aeroForce, aeroMoment] =... + aeroResultant(length, armTip, aeroAmplitude, aeroDistribution) + +% aeroResultant - Computes the 3D resultant in force and moment wrt CG of a +% portion of rocket, given a selected pressure distribution +% +% INPUTS: +% - length, double [1], length of the portion [m]; +% - armTip, double [1], distance from upper +% (towards nosecone) interface +% of the portion to CG [m]; +% - aeroAmplitude, double [3, 1], amplitudes of the pressure +% distribution in x, y and z +% body axes [N/m]; +% - aeroDistribution, char, pressure distribution type +% (only 'uniform' for now); +% +% OUTPUTS: +% - aeroForce, double [3, 1], force resultant in body axes [N]; +% - aeroMoment, double [3, 1], moment resultant in body axes [Nm]. +% +% REVISIONS: +% - 0 16-07-2024, Release, Filippo Cataldo + +switch aeroDistribution + case 'uniform' + + aeroForce = integral(@(x) -aeroAmplitude, -(length - armTip),... + armTip, 'ArrayValued', true); + aeroMoment = integral(@(x) cross([x; 0; 0], -aeroAmplitude),... + -(length - armTip), armTip, 'ArrayValued', true); +end + +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/src/eqResidual.m b/RoccarasoFlight/postprocessing/ASD/src/eqResidual.m new file mode 100644 index 0000000000000000000000000000000000000000..1c819c0f55fee8b9eaee2dc439e55cb4330063dd --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/src/eqResidual.m @@ -0,0 +1,60 @@ +function res = eqResidual(var, dStateBody, gBody, aeroDistribution, rocket, settings) + +% eqResidual - Computes the residual in the 3D equations of equilibrium of +% the whole rocket, along y and z body axes +% +% INPUTS: +% - var, double [4, 1], variable of the system: +% y and z parachute force's components, double [2, 1], [N]; +% y and z amplitudes of pressure, double [2, 1], [N/m]; +% - dStateBody, double [12, 1], 3D state's derivative: +% CG velocity, double [3, 1], [m/s]; +% CG acceleration, double [3, 1], [m/s^2]; +% angular velocity, double [3, 1], [rad/s]; +% angular acceleration, double [3, 1], [rad/s^2]; +% - gBody, double [3, 1], gravity acceleration [m/s^2]; +% - aeroDistribution, char, pressure distribution type +% (only 'uniform' for now); +% - rocket, object struct, rocket's geometry and mass; +% - settings, object struct, simulation's parameters; +% +% OUTPUTS: +% - res, double [4, 1], residual in displacement and +% rotation equilibrium [N]. +% +% REVISIONS: +% - 0 16-07-2024, Release, Filippo Cataldo + +%%% Get data +nameFirstBay = settings.nameFirstBay; % Name of tip's bay +indexFirstBay = find(keys(rocket.bays) == nameFirstBay); + +% Position of each bay from the base of the nosecone +absolutePositions = values(rocket.absolutePositions); +absolutePositions = absolutePositions(indexFirstBay:end); % [m] + +length = rocket.length; % [m] +xCg = rocket.xCg(end); % [m] +armTip = xCg - absolutePositions(1); % [m] +m = rocket.mass(end); % [kg] +inertia = diag(rocket.inertia(:,end)); % [kg*m^2] + +accBody = dStateBody(4:6); +omegaBody = dStateBody(7:9); +omegaDotBody = dStateBody(10:12); + +% Parachute's force and pressure amplitudes, x component is not used +% since it has no arm in moment equilibrium +force = [0; var(1:2)]; % [N] +aeroAmplitude = [0; var(3:4)]; % [N/m] + +%%% Evaluate equations of equilibrium +[aeroForce, aeroMoment] = aeroResultant(length, armTip, aeroAmplitude,... + aeroDistribution); + +resForce = force + aeroForce - m*(accBody - gBody); +resMoment = cross([armTip; 0; 0], force) + aeroMoment - ... + inertia*omegaDotBody - cross(omegaBody, inertia*omegaBody); +res = [resForce(2:3); resMoment(2:3)]; % Select only y and z components + +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/src/internalLoads.m b/RoccarasoFlight/postprocessing/ASD/src/internalLoads.m new file mode 100644 index 0000000000000000000000000000000000000000..c4ee8be6993ba366d78dc0ff76720114818cdd90 --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/src/internalLoads.m @@ -0,0 +1,134 @@ +function [maxLoads, loads, time, state, dState] =... + internalLoads(state0, rocket, settings) + +% internalLoads - Integrates the ODE of 2D parachute's deployment. +% Then obtains the internal loads at every interface of the +% rocket. Quantities are in inertial frame. +% +% INPUTS: +% - state0, double [4, 1], initial state of the rocket: +% altitude, double [1], [m]; +% CG velocity, double [2, 1], [m/s]; +% attack angle, double [1], [rad]; +% - rocket, object struct, geometry and mass of the rocket; +% - settings, object struct, simulation's parameters; +% +% OUTPUTS: +% - maxLoads, double [3, nBays], max loads at interfaces: +% axial, double [1, nBays], [N]; +% shear, double [1, nBays], [N]; +% moment, double [1, nBays], [Nm]; +% - loads, double [3, nBays, nTime], loads in time at interfaces +% with same structure as +% maxLoads; +% - time, double [1, nTime], ODE time vector [s]; +% - state, double [5, nTime], states at every step: +% altitude, double [1, nTime], [m]; +% CG velocity, double [2, nTime], [m/s]; +% pitch angle, double [1, nTime], [rad]; +% angular velocity, double [1, nTIme], [rad/s]; +% - dState, double [5, nTime], state's time derivative. +% +% CALLED FUNCTIONS: parachuteDeployment.m +% +% REVISIONS: +% - 0 Release, Luca Erbetta, Lorenzo Ciuti +% - 1 Refactoring, Filippo Cataldo +% Changes: Refactored the structure of the code. +% - 2 16-07-2024, Refactoring, Filippo Cataldo, +% Antonio Crisalli, Sasha Dell'Orto +% Changes: Refactored the structure of the code. + +%%% Get data +odeOptions = settings.odeOptions; +g = settings.g; % Gravity acceleration vector [m/s^2] +time = settings.time; % Time array of the simulation [s] + +% Array with all the rocket's bays +nameFirstBay = settings.nameFirstBay; % Name of tip's bay +indexFirstBay = find(keys(rocket.bays) == nameFirstBay); +bays = values(rocket.bays); +bays = bays(indexFirstBay:end); +nBays = length(bays); + +% Position of each bay from the base of the nosecone +absolutePositions = values(rocket.absolutePositions); +absolutePositions =... % [m] + absolutePositions(indexFirstBay:end); + +% Position of rocket's CG from the base of the nosecone +xCg = rocket.xCg(end); % [m] + +%%% Initial state +alt0 = state0(1); % Initial altitude [m] +vel0 = state0(2:3); % Initial CG velocity [m/s] +AoA0 = state0(4); % Initial attack angle [rad] + +theta0 = atan2(vel0(2), vel0(1)) + AoA0; % Initial pitch angle [rad] +omega0 = g(2)*vel0(1)/norm(vel0)^2; % Initial angular velocity [rad/s] + +state0 = [alt0; vel0; theta0; omega0]; % Rearrange initial state + +%%% Solve differential equation +[time, state] = ode113(@(time, state)... + parachuteDeployment(time, state, rocket, settings),... + time, state0, odeOptions); +time = time'; +state = state'; +dState = zeros(size(state)); + +%%% Compute internal loads +nTime = length(time); +loads = zeros(3, nBays + 1, nTime); + +for i = 1:nTime + + iTheta = state(4,i); + iOmega = state(5,i); + + % Rotation matrix from inertial to body frame + iRib = [cos(iTheta) sin(iTheta) + -sin(iTheta) cos(iTheta)]; + iRbi = iRib'; % Rotation matrix from body to inertial frame + + dState(:,i) = parachuteDeployment(time(i), state(:,i), rocket,... + settings); + iAcc = dState(2:3,i); + iOmegaDot = dState(5,i); % refactoring magic! + + for j = nBays:-1:1 + + % Extract parameters of current bay + absolutePositionBay = absolutePositions(j); + + if j < nBays + lBay = absolutePositions(j+1) - absolutePositionBay; + else + lBay = rocket.length - absolutePositionBay +... + absolutePositions(1); + end + + xCgBay = bays(j).xCg(end); + armBay = xCg - (absolutePositionBay + xCgBay); + mBay = bays(j).mass(end); + IBay = bays(j).inertia(3,end); + accBay = iAcc + iRbi*armBay*[-iOmega^2; iOmegaDot]; + + % Compute loads + loads(1:2,j,i) = iRib*mBay*(accBay - g) +... % Force [N] + loads(1:2,j+1,i); + loads(3,j,i) = IBay*iOmegaDot + loads(3,j+1,i) -... % Moment [Nm] + (lBay - xCgBay)*loads(2,j+1,i) - xCgBay*loads(2,j,i); + + end + +end + +% There are (nBays + 1) interfaces, but we don't need the last one (free +% interface) +loads = loads(:,1:nBays,:); + +% Maximum loads for every bay and every combination +maxLoads = max(abs(loads), [], 3); + +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/src/parachuteDeployment.m b/RoccarasoFlight/postprocessing/ASD/src/parachuteDeployment.m new file mode 100644 index 0000000000000000000000000000000000000000..040469bcc8c994aa778ecd709a905781db345a8f --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/src/parachuteDeployment.m @@ -0,0 +1,67 @@ +function dState = parachuteDeployment(time, state, rocket, settings) + +% parachuteDeployment - ODE of 2D parachute's deployment. +% Quantities are in inertial frame. +% +% INPUTS: +% - time, double [1], time instant from deployment [s]; +% - state, double [5, 1], state of the rocket: +% altitude, double [1], [m]; +% CG velocity, double [2], [m/s]; +% pitch angle, double [1], [rad]; +% angular velocity, double [1], [rad/s]; +% - rocket, object struct, geometry and mass of the rocket; +% - settings, object struct, simulation's parameters; +% +% OUTPUTS: +% - dState, double [5, nTime], state's time derivative. +% +% CALLED FUNCTIONS: parachuteForce.m +% +% REVISIONS: +% - 0 Release, Luca Erbetta, Lorenzo Ciuti +% - 1 Refactoring, Filippo Cataldo +% Changes: Refactored the structure of the code. +% - 2 16-07-2024, Refactoring, Filippo Cataldo, +% Antonio Crisalli, Sasha Dell'Orto +% Changes: Refactored the structure of the code. + +%%% Get data +g = settings.g; % Gravity acceleration vector [m/s^2] +nameFirstBay = settings.nameFirstBay; % Name of tip's bay +indexFirstBay = find(keys(rocket.bays) == nameFirstBay); + +% Position of each bay from the base of the nosecone +absolutePositions = values(rocket.absolutePositions); +absolutePositions =... % [m] + absolutePositions(indexFirstBay:end); + +% Position of rocket's CG from the base of the nosecone +xCg = rocket.xCg(end); % [m] + +% Position of rocket's tip wrt CG +armTip = xCg - absolutePositions(1); % [m] + +m = rocket.mass(end); % [kg] +I = rocket.inertia(3,end); % Inertia wrt 3rd body axis [kg*m^2] + +%%% Extract state's variables +velY = state(3); % CG vertical velocity [m/s] +theta = state(4); % Pitch angle [rad] +omega = state(5); % Angular velocity [rad/s] + +Rib = [cos(theta) sin(theta) % Rotation matrix from inertial to body frame + -sin(theta) cos(theta)]; + +%%% Compute state's derivative +force = parachuteForce(time, state, rocket, settings); + +% Moment wrt CG, using second component of parachute's force in body axes +moment = armTip*Rib(2,:)*force; + +acc = g + force/m; +omegaDot = moment/I; + +dState = [velY; acc; omega; omegaDot]; + +end \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/ASD/src/parachuteForce.m b/RoccarasoFlight/postprocessing/ASD/src/parachuteForce.m new file mode 100644 index 0000000000000000000000000000000000000000..5a3e1a949036b33100c5a8c8dda558c8aa1b372e --- /dev/null +++ b/RoccarasoFlight/postprocessing/ASD/src/parachuteForce.m @@ -0,0 +1,81 @@ +function force = parachuteForce(time, state, rocket, settings) + +% parachuteForce - Parachute's force at given state and time. Quantities +% are in inertial frame. +% +% INPUTS: +% - time, double [1], time instant [s]; +% - state, double [5, 1], state of the rocket: +% altitude, double [1], [m]; +% CG velocity, double [2], [m/s]; +% pitch angle, double [1], [rad]; +% angular velocity, double [1], [rad/s]; +% - rocket, object struct, geometry and mass of the rocket; +% - settings, object struct, simulation's parameters; +% +% OUTPUTS: +% - force, double [2, 1], parachute's force [N]. +% +% REVISIONS: +% - 0 Release, Luca Erbetta, Lorenzo Ciuti +% - 1 Refactoring, Filippo Cataldo +% Changes: Refactored the structure of the code. +% - 2 16-07-2024, Refactoring, Filippo Cataldo, +% Antonio Crisalli, Sasha Dell'Orto +% Changes: Refactored the structure of the code. + +%%% Get data +time0 = settings.time(1); % Initial time of the simulation [s] +nameFirstBay = settings.nameFirstBay; % Name of tip's bay +indexFirstBay = find(keys(rocket.bays) == nameFirstBay); + +% Position of each bay from the base of the nosecone +absolutePositions = values(rocket.absolutePositions); +absolutePositions =... % [m] + absolutePositions(indexFirstBay:end); + +% Position of rocket's CG from the base of the nosecone +xCg = rocket.xCg(end); % [m] + +% Position of rocket's tip wrt CG +armTip = xCg - absolutePositions(1); % [m] + +% Parachute's parameters +parachuteName = settings.parachuteName; + +switch parachuteName + case 'DROGUE chute' + cd = rocket.parachutes(1,1).cd; + forceCoefficient = rocket.parachutes(1,1).forceCoefficient; + openingTime = rocket.parachutes(1,1).openingTime; + surfaceFinal = rocket.parachutes(1,1).surface; + case 'MAIN chute' + cd = rocket.parachutes(2,1).cd; + forceCoefficient = rocket.parachutes(2,1).forceCoefficient; + openingTime = rocket.parachutes(2,1).openingTime; + surfaceFinal = rocket.parachutes(2,1).surface; +end + +%%% Extract state's variables +alt = state(1); % Altitude [m] +vel = state(2:3); % CG velocity [m/s] +theta = state(4); % Pitch angle [rad] +omega = state(5); % Angular velocity [rad/s] + +%%% Compute parachute's force with aerodynamic drag formula +% Parachute's surface increases linearly with time until surfaceFinal +if (time - time0) < openingTime + surface = (time - time0)/openingTime*surfaceFinal; +else + surface = surfaceFinal; +end + +% Air density of standard atmosphere +[~, ~, ~, rho] = atmosisa(alt); % [kg/m^3] + +% Rocket's tip velocity (rigid body) +velTip = vel + omega*armTip*[-sin(theta); cos(theta)]; % [m/s] + +force = -0.5*rho*surface*cd*forceCoefficient*norm(velTip)*velTip; + +end \ No newline at end of file