Skip to content
Snippets Groups Projects
Commit 1ad1a2ff authored by gabrigghia's avatar gabrigghia Committed by Marco Luigi Gaibotti
Browse files

first changes, to be reviewed 

pusho solo cosi da non perdere niente nell'aggiornamento delle branch
parent af3955c3
Branches
No related tags found
1 merge request!11Merge code-refactoring into master
Show whitespace changes
Inline Side-by-side
stabilityAnalysis/gabrleggi.m 0 → 100644
+ 10
− 0
View file @ 1ad1a2ff
% % % % % togliere qualunque cosa correlata con openrocket. togliere qualunque cosa
% % % % % con or alla fine che fa riferimento a or
% % % %
% guardare ballistic per vedere come mettere la descrizione inputs ouputs
% %
% % andare nel drive afd utility matlab code format
%
% % fai copia incolla delle ose nel file, mettere su matlab in favorites new faovrite, cambiare label in align comments
% % guardare anche standard header for functions
% non mettere varargin vrarargout ma mettere come arguments in ballistic
\ No newline at end of file
stabilityAnalysis/mainStability.m 0 → 100644
+ 345
− 0
View file @ 1ad1a2ff
function mainStability(rocket, wind, environment)
arguments
rocket Rocket = Rocket.empty
wind WindCustom = WindCustom.empty
environment Environment = Environment.empty
end
%{
mainStability - Main script to compute the rocket stability at
launchpad exit given a set of wind intensity and
azimuth
REVISIONS:
- 0 12/11/2022, release Riccardo Cadamuro
- 1 03/12/2022, update Maria Teresa Cazzola, Riccardo Cadamuro
SM computed using also alpha total and phif
- 2 06/04/2023, update Riccardo Cadamuro
added wind magnitude dependence
Copyright © 2021, Skyward Experimental Rocketry, AFD department
All rights reserved
SPDX-License-Identifier: GPL-3.0-or-later
%}
%% PATH
currentPath = fileparts(mfilename("fullpath"));
commonPath = trimPath(fullfile(currentPath, '..', '..', 'common'));
if ~contains(path, currentPath), addpath(genpath(currentPath)); end
if ~contains(path, commonPath), addpath(genpath(commonPath)); end
filePath = fileparts(mfilename('fullpath'));
currentPath = pwd;
if not(strcmp(filePath, currentPath))
cd (filePath);
currentPath = filePath;
end
addpath(genpath(currentPath));
%% LOAD DATA / SETTINGS
mission = Mission('2024_Lyra_Portugal_October');
if isempty(rocket), rocket = Rocket(mission); end
if isempty(environment), environment = Environment(mission, rocket.motor); end
if isempty(wind), wind = WindCustom(mission); end
if vars.externalSource && ~vars.XCPreg
settings.stability.xcg = Stability.rocket.xcg(1);
settings.stability.Lcenter = Stability.rocket.Lcenter;
settings.stability.C = Stability.rocket.C;
% fins
settings.Chord1 = Stability.rocket.fins.chord1;
settings.Chord2 = Stability.rocket.fins.chord2;
settings.height = Stability.rocket.fins.height;
settings.deltaXLE = Stability.rocket.fins.deltaXLE;
% ogive
settings.OgType = 'MHAACK';
settings.cMod = Stability.rocket.ogive.cMod;
settings.pMod = Stability.rocket.ogive.pMod;
settings.Lnose = Stability.rocket.ogive.Lnose;
% boat
settings.boatL = Stability.rocket.boatL;
settings.boatD = Stability.rocket.boatD;
Q0 = angleToQuat(settings.PHI, settings.OMEGA, 180*pi/180)';
Ypad = 0;
input = createDissileInput(settings, dissileVars);
input.refq.XCG = settings.xcg;
else
%% LAUNCHPAD DYNAMICS
% Selecting the most constraining condition (max wind and azimuth tents
% to 180 deg)
settings.wind.altitudes = 0;
settings.wind.magVec = max(vars.windMag);
[~, iAz] = min(abs(vars.windAz - pi));
settings.wind.azVec = vars.windAz(iAz);
[uw, vw, ww] = computeWindVels(settings.wind, 0);
settings.constWind = [uw, vw, ww];
% states to compute the exit pad velocity
T = (288.15 - 0.0065*settings.z0); % temperature
a = sqrt(T*1.4*287.0531); % sound speed @launchpad
X0 = [0 0 0]';
V0 = [0 0 0]';
W0 = [0 0 0]';
Q0 = angleToQuat(settings.PHI, settings.OMEGA, 180*pi/180)';
%%% Initial State
Y0 = [X0; V0; W0; Q0];
%%% ode run
[Tpad, Ypad] = ode113(@ascent, [0, 10], Y0, settings.ode.optionspad,settings);
%%% saving data
vExit = Ypad(end, 4);
machExit = vExit/a;
input = createDissileInput(settings, dissileVars);
input.fltcon.MACH = machExit; % Mach @launchpad
input.fltcon.ALT = settings.z0; % Altitude @launchpad
input.refq.XCG = interpLinear(settings.motor.expTime, settings.xcg, Tpad(end));
end
%% COMPUTING THE LAUNCHPAD STABILITY DATA
tic
data_stability = computeStability(vars, input, Ypad, Q0);
time = toc;
%% REMOVE NON PHISICAL POINTS
indexNan = or(data_stability.XCPmat < vars.SMlim(1), data_stability.XCPmat > vars.SMlim(2));
data_stability.XCPmat(indexNan) = nan;
%% RETRIVE DATA
XCPlon = squeeze(data_stability.XCPmat(1, :, :));
XCPlat = squeeze(data_stability.XCPmat(2, :, :));
XCPtotSTD = squeeze(data_stability.XCPmat(3, :, :));
XCPtotMOD = squeeze(data_stability.XCPmat(4, :, :));
XCPlonATOT = squeeze(data_stability.XCPmat(5, :, :));
XCPlatATOT = squeeze(data_stability.XCPmat(6, :, :));
XCPor = squeeze(data_stability.XCPmat(7, :, :));
%%
if vars.prints || vars.plots
flag = true;
if length(vars.windMag) == 1
refVec = vars.windAz.*180/pi;
nameVec = {'Wind Azimuth','[deg]'};
const = vars.windMag;
nameConst = {'Wind Magnitude','[m/s]'};
elseif length(vars.windAz) == 1
refVec = vars.windMag;
nameVec = {'Wind Magnitude', '[m/s]'};
const = vars.windAz.*180/pi;
nameConst = {'Wind Azimuth', '[deg]'};
else
flag = false;
warning('Data printing not allowed');
end
end
%% PRINT
if vars.prints
if flag
nEl = length(refVec);
for i = 1:nEl
fprintf('--------------------------- Case %d ---------------------------\n',i);
fprintf(strcat(nameVec{1},': %5.2f ',nameVec{2}, ' | ',nameConst{1},': %5.2f ', nameConst{2}, '\n'), refVec(i), const);
fprintf('[ v-IN = %5.1f | v-AER = %5.1f m/s | M = %5.1f]\n',vExit, norm(data_stability.vRelMat(:, i)), norm(data_stability.vRelMat(:, i))/a);
fprintf('[alpha = %5.1f° | beta = %5.1f°]\n', data_stability.alphaExit(i), data_stability.betaExit(i));
fprintf('[alpha total = %5.1f° | phi = %5.1f°]\n', data_stability.alphaTotExit(i), data_stability.phiExit(i));
fprintf('\tSM longitudinal = %5.3f\n', XCPlon(i));
fprintf('\tSM lateral = %5.3f\n', XCPlat(i));
fprintf('\tSM total STD = %5.3f\n', XCPtotSTD(i));
fprintf('\tSM total MOD = %5.3f\n', XCPtotMOD(i));
fprintf('\tSM longitudinal ATOT = %5.3f\n', XCPlonATOT(i));
fprintf('\tSM lateral ATOT = %5.3f\n', XCPlatATOT(i)')
if vars.checkOpenRocket
fprintf('SM OpenRocket = %5.3f\n', XCPor(i));
end
fprintf('\n');
fprintf('\tAngle between Faero and Maero: %5.3f [deg]\n', data_stability.qIndexMat(1,i)*180/pi);
fprintf('\tAngle between Fyz and Myz: %5.3f [deg]\n', data_stability.qIndexMat(2,i)*180/pi);
fprintf('\tAngle between Fyz and vRyz: %5.3f [deg]\n', data_stability.qIndexMat(3,i)*180*pi);
fprintf('\tT - index: %5.3f [deg]\n', data_stability.qIndexMat(4,i));
fprintf('\tNorm of Fpar: %5.3f\n', data_stability.qIndexMat(5,i));
fprintf('\n\n');
end
end
fprintf('Computational Effort = %5.3f [s]\n', time);
end
%% PLOTS
if vars.plots
if flag
figure('Name','Stability margin');
hold on;
if (vars.plotXCPlon) plot(refVec, XCPlon, '-o', 'DisplayName','XCPlon'); end
if (vars.plotXCPlat) plot(refVec, XCPlat, '-o', 'DisplayName','XCPlat'); end
if (vars.plotXCPtotSTD) plot(refVec, XCPtotSTD, '-o', 'DisplayName', 'XCPtot - STD'); end
if (vars.plotXCPtotMOD) plot(refVec, XCPtotMOD, '-o', 'DisplayName', 'XCPtot - MOD'); end
if (vars.plotXCPlonATOT) plot(refVec, XCPlonATOT, '-o', 'DisplayName', 'XCPlon - ATOT'); end
if (vars.plotXCPlatATOT) plot(refVec, XCPlatATOT, '-o', 'DisplayName', 'XCPlat - ATOT'); end
if (vars.checkOpenRocket) plot(refVec, XCPor, '-o', 'DisplayName', 'XCP Open Rocket'); end
grid on;
xlabel(strcat(nameVec{1},{' '}, nameVec{2}));
ylabel('SM [-]');
title('STABILITY MARGIN');
subtitle(strcat(nameConst{1},':',{' '},num2str(const), {' '}, nameConst{2}));
legend;
if vars.plotQINDEX
figure('Name','Quality indices (angular)')
names = {'Angle between Faero and Maero', 'Angle between Fyz and Myz', 'Angle between Fyz and vRyz'};
hold on;
for i = 1:3
plot(refVec, data_stability.qIndexMat(i, :).*180/pi,'-o','DisplayName',names{i});
end
grid on;
xlabel(strcat(nameVec{1},{' '}, nameVec{2}));
ylabel('Angle [deg]');
title('Angular indices');
subtitle(strcat(nameConst{1},':',{' '},num2str(const), {' '}, nameConst{2}));
legend;
figure('Name','Tindex')
plot(refVec, data_stability.qIndexMat(4, :),'-o', 'DisplayName','T - index');
grid on;
xlabel(strcat(nameVec{1},{' '}, nameVec{2}));
ylabel('[-]');
title('T - index')
subtitle(strcat(nameConst{1},':',{' '},num2str(const), {' '}, nameConst{2}));
legend;
figure('Name','Quality index (force)')
plot(refVec, data_stability.qIndexMat(5, :),'-o', 'DisplayName','Norm of Fpar');
grid on;
xlabel(strcat(nameVec{1},{' '}, nameVec{2}));
ylabel('[-]');
title('Force index');
end
if vars.plotCoeffs
names = {'Cl', 'Cm', 'Cn', 'CA', 'CY', 'CN'};
ind = [9 3 7 10 8 4];
figure('Name','Coefficients');
hold on;
for j = 1:6
coeff = data_stability.stabilityCoeffs(ind(j), :);
coeff = coeff./max(abs(coeff));
if j <=3
mk = '-o';
else
mk = '-^';
end
plot(refVec, coeff, mk , 'DisplayName', names{j});
end
grid on
legend
ylabel('[-]');
xlabel(strcat(nameVec{1},{' '}, nameVec{2}));
title('Coefficients normalized');
end
else
names = {'XCPlon', 'XCPlat', 'XCPtot-STD', 'XCPtot-STD','XCPlon-ATOT','XCPlat-ATOT','XCPor'};
plotXCPflags =[vars.plotXCPlon, vars.plotXCPlat, vars.plotXCPtotSTD, vars.plotXCPtotMOD, vars.plotXCPlonATOT, vars.plotXCPlatATOT, vars.checkOpenRocket];
for i = 1:7
if plotXCPflags(i)
tit = strcat('Stability margin [', names{i}, ']');
figure('Name',tit);
SM = squeeze(data_stability.XCPmat(i, :, :));
surf(vars.windMag, vars.windAz*180/pi,SM','FaceColor','interp');
colormap jet;
cb = colorbar;
cb.Label.String = 'SM [-]';
xlabel('Wind Magnitude [m/s]');
ylabel('Wind Azimuth [deg]');
zlabel('SM [-]');
title(names{i});
end
end
if vars.plotQINDEX
names = {'Angle between Faero and Maero', 'Angle between Fyz and Myz', 'Angle between Fyz and vRyz', 'T - index', 'Norm of Fpar'};
for i = 1:5
tit = strcat('Quality index num:', num2str(i));
figure('Name',tit);
qInd = squeeze(data_stability.qIndexMat(i, :, :));
if i <= 3
qInd = qInd.*180/pi;
end
surf(vars.windMag, vars.windAz*180/pi,qInd','FaceColor','interp');
colormap jet;
cb = colorbar;
xlabel('Wind Magnitude [m/s]');
ylabel('Wind Azimuth [deg]');
if i <= 4
cb.Label.String = 'qIndex [deg]';
zlabel('qIndex [deg]');
else
cb.Label.String = 'qIndex [-]';
zlabel('qIndex [-]');
end
title(names{i})
end
end
end
end
%% SAVING
if vars.saveRes
alphaExit = data_stability.alphaExit;
betaExit = data_stability.betaExit;
phiExit = data_stability.phiExit;
vRelMat = data_stability.vRelMat;
vRelCFD = vRelMat;
vRelCFD(1, :) = -vRelMat(1, :);
if vars.useAlphaTot
[alphaExit, betaExit] = getAlphaBeta(alphaExit*pi/180, phiExit*pi/180);
alphaExit = alphaExit * 180/pi;
betaExit = betaExit * 180/pi;
end
stabilityRes.fltcon.alpha = alphaExit;
stabilityRes.fltcon.beta = betaExit;
stabilityRes.fltcon.mach = machExit;
stabilityRes.fltcon.alt = settings.z0;
stabilityRes.fltcon.vRelSIM = vRelMat;
stabilityRes.fltcon.vRelCFD = vRelCFD;
stabilityRes.xcg = input.refq.XCG;
stabilityRes.XCPlon = XCPlon;
stabilityRes.XCPlat = XCPlat;
stabilityRes.XCPtotSTD = XCPtotSTD;
stabilityRes.XCPtotMOD = XCPtotMOD;
stabilityRes.XCPlonATOT = XCPlonATOT;
stabilityRes.XCPlatATOT = XCPlatATOT;
stabilityRes.XCPor = XCPor;
stabilityRes.coeff = data_stability.stabilityCoeffs;
save('stabilityRes.mat', 'stabilityRes');
fprintf('RESULT SAVED!\n')
end
if not(vars.externalSource)
clearvars -except data_stability XCPlon XCPlat XCPtotSTD XCPtotMOD XCPlonATOT XCPlatATOT XCPor vars flag settings
end
%% LOAD DATA / SETTINGS
mission = Mission('2024_Lyra_Portugal_October');
if isempty(rocket), rocket = Rocket(mission); end
if isempty(environment), environment = Environment(mission, rocket.motor); end
if isempty(wind), wind = WindCustom(mission); end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment