diff --git a/RoccarasoFlight/postprocessing/ASD/plots.m b/RoccarasoFlight/postprocessing/ASD/plots.m
index 79288cb13410a1548c8807644c4d605792b8f449..5a4c522c4b61be43b7eb70c66b740405b07e1a6d 100644
--- a/RoccarasoFlight/postprocessing/ASD/plots.m
+++ b/RoccarasoFlight/postprocessing/ASD/plots.m
@@ -33,7 +33,7 @@ for i = 1:nBays
' - ']), 'FontSize', 12);
end
- xlim('Tight');
+ xlim('tight');
xlabel('Time [s]', 'FontSize', 12);
ylabel('Force [N], Moment [Nm]', 'FontSize', 12);
legend('Measured N', 'Simulated N', 'Measured S',...
diff --git a/RoccarasoFlight/postprocessing/ASD/postProcessDrogue.m b/RoccarasoFlight/postprocessing/ASD/postProcessDrogue.m
index 18b4c42182c276966ba5564a43d9c81f4d24e736..e3c5265e7f388a7f29df3d2004adfae4e214598a 100644
--- a/RoccarasoFlight/postprocessing/ASD/postProcessDrogue.m
+++ b/RoccarasoFlight/postprocessing/ASD/postProcessDrogue.m
@@ -2,12 +2,12 @@
% 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.
+% IMU time is used for interpolation.
% 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.
+% frame. Body frame can be 2D or 3D (IMU frame)
+% 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
%
@@ -26,11 +26,24 @@ load("data\geminiRocketRoccaraso.mat");
load("data\geminiFlightRoccaraso.mat");
%% Set parameters for the post process
+duration = 2.5; % Time window duration of the analysis [s]
+
% Pressure distribution type, only 'uniform' supported for now
aeroDistribution = 'uniform';
-indexOpen = 968; % Drogue opening esimtated index in quaternions time
-gNED = [0; 0; 9.81]; % Gravity acceleration in NED [m/s^2]
+% Choose on which time span interpolate, 'NAS' or 'IMU'. Latter one with
+% higher sampling frequency
+interpOn = 'IMU';
+
+switch interpOn
+ case 'NAS'
+ indexOpen = 968; % Drogue opening esimtated index in NAS time
+ time = NAS.time'; % NAS time is used
+ case 'IMU'
+ indexOpen = 8526; % Drogue opening esimtated index in IMU time
+ time = IMU.accTime'; % IMU time is used
+end
+
imu2bay = 0.1854; % IMU longitudinal distance from top of its bay (ELC) [m]
% IMU y and z position from longitudinal axis, sign according to IMU
@@ -39,22 +52,24 @@ yImu = 0; % [m]
zImu = -0.075; % diameter/2 [m]
% Options for non-linear system
-guess = [200*ones(2,1); ones(2,1)]; % Initial guess
+guess = [200*ones(2,1); ones(2,1)]; % Initial guess, (see eqResidual.m)
nonlinOptions = optimoptions('fsolve');
nonlinOptions.Algorithm = 'levenberg-marquardt';
+gNED = [0; 0; 9.81]; % Gravity acceleration in NED [m/s^2]
+
%% 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.time = time(indexOpen) + [0 duration]; % Simulation time [s]
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);
+indexFirstBay = find(namesBays == 'recovery');
bays = bays(indexFirstBay:end);
namesBays = namesBays(indexFirstBay:end);
nBays = length(bays);
@@ -67,15 +82,15 @@ diameter = rocket.diameter; % [m]
xCg = rocket.xCg(end); % [m]
armImu = xCg - (rocket.absolutePositions('electronics') +...
imu2bay); % Accelerometer position from xCg [m]
+armTip = xCg - absolutePositions(1); % xCg distance from rocket's tip [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);
+indexEnd = indexOpen + ceil(duration/dt);
% Quaternions [-]
quat = zeros(4,1,nTime);
@@ -150,7 +165,6 @@ 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; yImu; zImu]) - cross(omegaBody(:,i),...
cross(omegaBody(:,i), [armImu; yImu; zImu]));
@@ -164,8 +178,8 @@ for i = 1:nTime
% Solve non-linear system to impose global equilibrium through
% aerodynamic forces
- sol = fsolve(@(var) eqResidual(var, dStateBody(:,i), iRnb*gNED,...
- aeroDistribution, rocket, settings),...
+ sol = fsolve(@(var) eqResidual(var, dStateBody(:,i), armTip,...
+ iRnb*gNED, aeroDistribution, rocket),...
guess, nonlinOptions);
aeroAmplitude = [0; sol(3:4)];
@@ -210,7 +224,9 @@ for i = 1:nTime
end
-% Free interface is not of interest
+% Free interface is not of interest, neither the rear bay
+nBays = nBays - 1;
+namesBays = namesBays(1:nBays);
measuredForces = measuredForces(:,1:nBays,:);
measuredMoments = measuredMoments(:,1:nBays,:);
measuredLoads = measuredLoads(:,1:nBays,:);
@@ -225,5 +241,8 @@ state0 = [alt0; velX0; velY0; AoA0];
[~, simulatedLoads, simulationTime] =...
internalLoads(state0, rocket, settings);
+% Don't consider rear bay
+simulatedLoads = simulatedLoads(:,1:nBays,:);
+
%% Plot results
plots;
\ No newline at end of file
diff --git a/RoccarasoFlight/postprocessing/ASD/postProcessMain.m b/RoccarasoFlight/postprocessing/ASD/postProcessMain.m
index 37f9b971bdaef1cd2ff3c685f1f122a61f5ceab3..5e7c30c3798653f58898d1a05eb31e9418c3b36d 100644
--- a/RoccarasoFlight/postprocessing/ASD/postProcessMain.m
+++ b/RoccarasoFlight/postprocessing/ASD/postProcessMain.m
@@ -2,7 +2,7 @@
% 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.
+% IMU time is used for interpolation.
% No specification of reference frame means 2D inertial
% frame. Body frame can be 2D or 3D (IMU frame) depending
% on the context. Anyway, frames are only used to rotate
@@ -26,11 +26,24 @@ load("data\geminiRocketRoccaraso.mat");
load("data\geminiFlightRoccaraso.mat");
%% Set parameters for the post process
+duration = 2.5; % Time window duration of the analysis [s]
+
% 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]
+% Choose on which time span to interpolate, 'NAS' or 'IMU'. Latter one with
+% higher sampling frequency
+interpOn = 'IMU';
+
+switch interpOn
+ case 'NAS'
+ indexOpen = 3585; % Main opening esimtated index in NAS time
+ time = NAS.time'; % NAS time is used
+ case 'IMU'
+ indexOpen = 31586; % Main opening esimtated index in IMU time
+ time = IMU.accTime'; % IMU time is used
+end
+
imu2bay = 0.1854; % IMU longitudinal distance from top of its bay (ELC) [m]
% IMU y and z position from longitudinal axis, sign according to IMU
@@ -39,22 +52,24 @@ yImu = 0; % [m]
zImu = -0.075; % diameter/2 [m]
% Options for non-linear system
-guess = [200*ones(2,1); ones(2,1)]; % Initial guess
+guess = [200*ones(2,1); ones(2,1)]; % Initial guess, (see eqResidual.m)
nonlinOptions = optimoptions('fsolve');
nonlinOptions.Algorithm = 'levenberg-marquardt';
+gNED = [0; 0; 9.81]; % Gravity acceleration in NED [m/s^2]
+
%% 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.time = time(indexOpen) + [0 duration]; % Simulation time [s]
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);
+indexFirstBay = find(namesBays == 'recovery');
bays = bays(indexFirstBay:end);
namesBays = namesBays(indexFirstBay:end);
nBays = length(bays);
@@ -67,15 +82,15 @@ diameter = rocket.diameter; % [m]
xCg = rocket.xCg(end); % [m]
armImu = xCg - (rocket.absolutePositions('electronics') +...
imu2bay); % Accelerometer x position from xCg [m]
+armTip = xCg - absolutePositions(1); % xCg distance from rocket's tip [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);
+indexEnd = indexOpen + ceil(duration/dt);
% Quaternions [-]
quat = zeros(4,1,nTime);
@@ -163,8 +178,8 @@ for i = 1:nTime
% Solve non-linear system to impose global equilibrium through
% aerodynamic forces
- sol = fsolve(@(var) eqResidual(var, dStateBody(:,i), iRnb*gNED,...
- aeroDistribution, rocket, settings),...
+ sol = fsolve(@(var) eqResidual(var, dStateBody(:,i), armTip,...
+ iRnb*gNED, aeroDistribution, rocket),...
guess, nonlinOptions);
aeroAmplitude = [0; sol(3:4)];
@@ -209,7 +224,9 @@ for i = 1:nTime
end
-% Free interface is not of interest
+% Free interface is not of interest, neither the rear bay
+nBays = nBays - 1;
+namesBays = namesBays(1:nBays);
measuredForces = measuredForces(:,1:nBays,:);
measuredMoments = measuredMoments(:,1:nBays,:);
measuredLoads = measuredLoads(:,1:nBays,:);
@@ -224,5 +241,8 @@ state0 = [alt0; velX0; velY0; AoA0];
[~, simulatedLoads, simulationTime] =...
internalLoads(state0, rocket, settings);
+% Don't consider rear bay
+simulatedLoads = simulatedLoads(:,1:nBays,:);
+
%% Plot results
plots;
\ No newline at end of file
diff --git a/RoccarasoFlight/postprocessing/ASD/src/eqResidual.m b/RoccarasoFlight/postprocessing/ASD/src/eqResidual.m
index 1c819c0f55fee8b9eaee2dc439e55cb4330063dd..94fa4c6850bfaf6be6660aa818ab6ae39f5f46e6 100644
--- a/RoccarasoFlight/postprocessing/ASD/src/eqResidual.m
+++ b/RoccarasoFlight/postprocessing/ASD/src/eqResidual.m
@@ -1,4 +1,5 @@
-function res = eqResidual(var, dStateBody, gBody, aeroDistribution, rocket, settings)
+function res =...
+ eqResidual(var, dStateBody, armTip, gBody, aeroDistribution, rocket)
% eqResidual - Computes the residual in the 3D equations of equilibrium of
% the whole rocket, along y and z body axes
@@ -12,11 +13,13 @@ function res = eqResidual(var, dStateBody, gBody, aeroDistribution, rocket, sett
% CG acceleration, double [3, 1], [m/s^2];
% angular velocity, double [3, 1], [rad/s];
% angular acceleration, double [3, 1], [rad/s^2];
+% - armTip, double [1], xCg distance from rocket's
+% tip in parachute's deployment
+% configuration [m];
% - 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
@@ -26,16 +29,7 @@ function res = eqResidual(var, dStateBody, gBody, aeroDistribution, rocket, sett
% - 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]