diff --git a/RoccarasoFlight/postprocessing/AFD/descent_para_plot.mat b/RoccarasoFlight/postprocessing/AFD/descent_para_plot.mat
new file mode 100644
index 0000000000000000000000000000000000000000..8be38e73c75f0ab0db3f7ff80ff1c408ec242cab
Binary files /dev/null and b/RoccarasoFlight/postprocessing/AFD/descent_para_plot.mat differ
diff --git a/RoccarasoFlight/postprocessing/AFD/estimationCA.m b/RoccarasoFlight/postprocessing/AFD/estimationCA.m
index d7f08131d7d65be148bd90a8ed5d58513e4c6e31..6f8df89b604f07b9498a9b3ff32eb1a3a1c5099c 100644
--- a/RoccarasoFlight/postprocessing/AFD/estimationCA.m
+++ b/RoccarasoFlight/postprocessing/AFD/estimationCA.m
@@ -4,11 +4,6 @@ function results = estimationCA(main, alt0)
NPOINTS = 1000;
MSTRUCT = 25.2885; % [kg] Gemini structural mass
AREF = 0.15^2 * pi * 0.25;
- PROP_FINAL_MASS = 0.7799 + 0.4445; % [kg] initial propellant mass
- tIGN = 0.97500; % Ignition time in the thrust curve
- tCO = 5.32702; % Cutoff time in the thrsut curve
- mDot = 0.9660; % [kg/s] Averaged propellant mass flow rate
- DELTA_TIME = 0.4988;
load("CA_alpha0.mat");
% acceleration raw
@@ -18,49 +13,25 @@ function results = estimationCA(main, alt0)
accXB = movmean(accXB_raw, 30);
% Motor data
- load("../RPS/cleanData/ROC2/engineFlightData.mat");
- timeThrust = flightData.timeThrust;
- indexTimeThrust = and(timeThrust > tIGN, timeThrust < tCO);
- timeThrust = timeThrust(indexTimeThrust);
- timeThrust = timeThrust - timeThrust(1) - DELTA_TIME;
-
- Thrust = flightData.Thrust(indexTimeThrust);
-
- tNodal = [timeThrust(1), -0.2248, 3.33322, timeThrust(end)];
- indexTransIGN = and(timeThrust>=tNodal(1), timeThrust<tNodal(2));
- indexTransCO = and(timeThrust>tNodal(3), timeThrust<=tNodal(4));
-
- thrustTransIGN = Thrust(indexTransIGN);
- thrustTransCO = Thrust(indexTransCO);
-
- offsetIGN = thrustTransIGN(1);
- offsetCO = thrustTransCO(end);
-
- thrustTransIGN = thrustTransIGN - offsetIGN;
- thrustTransCO = thrustTransCO - offsetCO;
-
- thrustTransIGN = thrustTransIGN * (1 + offsetIGN/thrustTransIGN(end));
- thrustTransCO = thrustTransCO * (1 + offsetCO/thrustTransCO(1));
-
- Thrust = [thrustTransIGN; Thrust(and(timeThrust>=tNodal(2), timeThrust<=tNodal(3))); thrustTransCO];
- mMotor = linspace(PROP_FINAL_MASS + mDot*(tCO - tIGN), PROP_FINAL_MASS, sum(indexTimeThrust));
- timeThrust = timeThrust - timeThrust(1);
-
+ motor = load("roccEngine.mat");
+ timeThrust = motor.time;
+ Thrust = motor.Thrust;
+ mMotor = motor.mass;
+
% Velocity
- tNas = main.NASData(1:main.NASApogeeIndex, 1);
- altNas = -main.NASData(1:main.NASApogeeIndex, 4);
- velNas = vecnorm(main.NASData(1:main.NASApogeeIndex, 5:7), 2, 2);
-
- load sensorPitot.mat
- indNasPit = sensorTot_WithPitot.nas.time < main.tApogee;
+ tNasFlight = main.NASData(1:main.NASApogeeIndex, 1);
+ altNasFlight = -main.NASData(1:main.NASApogeeIndex, 4);
+ velNasFlight = vecnorm(main.NASData(1:main.NASApogeeIndex, 5:7), 2, 2);
- tNasPit = sensorTot_WithPitot.nas.time(indNasPit);
- velNasPit = vecnorm(sensorTot_WithPitot.nas.states(indNasPit, 4:6), 2,2);
- altNasPit = -sensorTot_WithPitot.nas.states(indNasPit, 3);
+ load nasSim.mat
+ indNasSim = sensorTot_NoPitot.nas.time < main.tApogee;
+ tNasSim = sensorTot_NoPitot.nas.time(indNasSim);
+ velNasSim = vecnorm(sensorTot_NoPitot.nas.states(indNasSim, 4:6), 2,2);
+ altNasSim = -sensorTot_NoPitot.nas.states(indNasSim, 3);
% Normalize data
- timeRef = linspace(max([tNas(1), tAcc(1), tNasPit(1)]), min([tNas(end), tAcc(end), tNasPit(end)]) , NPOINTS);
+ timeRef = linspace(max([tNasFlight(1), tAcc(1), tNasSim(1)]), min([tNasFlight(end), tAcc(end), tNasSim(end)]) , NPOINTS);
indexEndThrust = find(timeRef > timeThrust(end), 1, 'first') - 1;
accXBRef = interp1(tAcc, accXB, timeRef);
@@ -71,68 +42,67 @@ function results = estimationCA(main, alt0)
mRocket = (mMotor(end) + MSTRUCT) * ones(1, NPOINTS);
mRocket(1:indexEndThrust) = interp1(timeThrust, mMotor, timeRef(1:indexEndThrust)) + MSTRUCT;
- velRef = interp1(tNas, movmean(velNas, 10), timeRef);
+ velNasFlightNorm = interp1(tNasFlight, movmean(velNasFlight, 10), timeRef);
+ velNasSimNorm = interp1(tNasSim, velNasSim, timeRef);
- velNasPitRef = interp1(tNasPit, velNasPit, timeRef);
-
- altRef = interp1(tNasPit, altNasPit, timeRef);
+ velVertNasFlightNorm = interp1(tNasFlight, movmean(-main.NASData(1:main.NASApogeeIndex, 7), 10), timeRef);
+ velVertNasSimNorm = interp1(tNasSim, -sensorTot_NoPitot.nas.states(indNasSim, 6), timeRef);
+
+ altFlightNorm = interp1(tNasFlight, altNasFlight, timeRef);
+ altSimNorm = interp1(tNasSim, altNasSim, timeRef);
% Air density and sonic velocity
- [~, vSonRef, ~, rhoRef] = atmosisa(altRef + alt0);
- machRef = velRef./vSonRef;
- machPitRef = velNasPitRef./vSonRef;
+ [~, vSonSimNorm, ~, rhoSimNorm] = atmosisa(altSimNorm + alt0);
+ [~, vSonFlightNorm, ~, rhoFlightNorm] = atmosisa(altFlightNorm + alt0);
+
+ machFlightNorm = velNasFlightNorm./vSonFlightNorm;
+ machSimNorm = velNasSimNorm./vSonSimNorm;
% index opening
indexOpeningARB = find(timeRef > 5.12, 1, 'first');
% Compute CA
- CA = (ThrustRef - mRocket.*accXBRef) ./ (0.5 .* rhoRef .* velRef.^2 * AREF);
+ CA_NasFlight = (ThrustRef - mRocket.*accXBRef) ./ (0.5 .* rhoFlightNorm .* velNasFlightNorm.^2 * AREF);
- CA_NasPitot = (ThrustRef - mRocket.*accXBRef) ./ (0.5 .* rhoRef .* velNasPitRef.^2 * AREF);
+ CA_NasSim = (ThrustRef - mRocket.*accXBRef) ./ (0.5 .* rhoSimNorm .* velNasSimNorm.^2 * AREF);
- tCO = tCO -1.8519;
indexCO = find(ThrustRef == 0, 1, 'first');
- mRocketNoPwr = mRocket(indexCO:end);
- CANoPwr = CA(indexCO:end);
CA_100ARB = zeros(1, length(timeRef));
CA_0ARB = zeros(1, length(timeRef));
for i = 1:length(timeRef)
- [~, indexMach] = nearestValSorted(MachRef, machPitRef(i));
- [~, indexAlt] = nearestValSorted(AltRef, altRef(i));
+ [~, indexMach] = nearestValSorted(MachRef, machSimNorm(i));
+ [~, indexAlt] = nearestValSorted(AltRef, altSimNorm(i));
CA_100ARB(i) = CA_MD(indexMach, indexAlt, 3);
CA_0ARB(i) = CA_MD(indexMach, indexAlt, 1);
end
- ThrustComp = 0.5 .* rhoRef(1:indexCO) .* velRef(1:indexCO).^2 .* AREF .* CA_0ARB(1:indexCO) + mRocket(1:indexCO).*accXBRef(1:indexCO);
+ ThrustComp = 0.5 .* rhoSimNorm(1:indexCO) .* velNasFlightNorm(1:indexCO).^2 .* AREF .* CA_0ARB(1:indexCO) + mRocket(1:indexCO).*accXBRef(1:indexCO);
- figure
- plot(timeRef(indexCO:end),CANoPwr);
- hold on
- plot(timeRef(indexCO:end), CA_NasPitot(indexCO:end));
- plot(timeRef(indexCO:end), [CA_0ARB(indexCO:indexOpeningARB) CA_100ARB(indexOpeningARB+1:end)]);
- xlabel('Time [s]');
- ylabel('CA [-]');
- title('CA from telemetry');
- grid on
- legend('NAS', 'NAS with PITOT', 'MD', 'Location','southoutside', 'Orientation', 'horizontal');
-
%% SAVE RESULTS
results.time = timeRef;
results.Thrust = ThrustRef;
results.ThrustComputed = ThrustComp;
- results.velocity.vel_pitot = velNasPitRef;
- results.velocity.val_noPitot = velRef;
-
- results.alt = altRef;
+ results.velocity.vel_NasSim = velNasSimNorm;
+ results.velocity.vel_NasFlight = velNasFlightNorm;
+ results.velocity.velVert_NasSim = velVertNasSimNorm;
+ results.velocity.velVert_NasFlight = velVertNasFlightNorm;
+
+ results.altitude.alt_NasSim = altSimNorm;
+ results.altitude.alt_NasFlight = altFlightNorm;
+
+ results.mach.mach_NasSim = machSimNorm;
+ results.mach.mach_NasFlight = machFlightNorm;
+
+
results.CA.CA100 = CA_100ARB;
results.CA.CA0 = CA_0ARB;
- results.CA.CA_pitot = CA_NasPitot;
- results.CA.CA_noPitot = CA;
+ results.CA.CA_NasSim = CA_NasSim;
+ results.CA.CA_NasFlight = CA_NasFlight;
results.mass.propMass = mMotor;
results.mass.rocketMass = mRocket;
diff --git a/RoccarasoFlight/postprocessing/AFD/mainPostprocess.m b/RoccarasoFlight/postprocessing/AFD/mainPostprocess.m
index 3616d5a4580deab3333a7d604fcc6a5dae1556b5..c7b1fde450c038776b8e966e81df3b186d692a5f 100644
--- a/RoccarasoFlight/postprocessing/AFD/mainPostprocess.m
+++ b/RoccarasoFlight/postprocessing/AFD/mainPostprocess.m
@@ -3,13 +3,14 @@ clc;
close all;
%% Select options
-opt.flagTraj = false;
+opt.flagTraj = true;
opt.flagCA = true;
alt0 = 1414;
%% adding path and recall data
addpath("..\commonFunctions\");
-addpath(genpath("..\..\..\..\msa-toolkit\commonFunctions\"));
+addpath(genpath("..\..\..\..\msa-toolkit\commonFunctions\"));
+addpath(genpath("..\..\..\..\msa-toolkit\simulator\src\"))
main = importData('main');
%% Trajectory
@@ -21,18 +22,18 @@ end
%% CA Estimation
if opt.flagCA
+ clc
res = estimationCA(main,alt0);
-
+ indexCO = res.events.shutdonw;
+ indexOpeningARB = res.events.ARB_opening;
end
%% TEST Simulation
if opt.flagCA
-
- indexCO = res.events.shutdonw;
- indexOpeningARB = res.events.ARB_opening;
configSimPostp;
+ settings.Local = [1414, 25 + 273, 88369];
% CONTINUE
Thrust = res.Thrust(1:indexCO);
@@ -60,8 +61,61 @@ if opt.flagCA
settings.motor.expThrust = Thrust;
settings.tb = settings.motor.expTime(end);
settings.tControl = settings.tb;
+
+ settings.z0 = settings.z0 - 100;
[Tf, Yf, Ta, Ya, bound_value] = stdRun(settings);
+ load ascent_plot.mat;
+
+ %%
+ close all;
+
+ figure
+ plot(res.time, res.altitude.alt_NasFlight);
+ hold on
+ plot(res.time, res.altitude.alt_NasSim);
+ plot(Ta, -Ya(:, 3))
+ xlabel('Time [s]');
+ ylabel('Altitude [m]');
+ title('Altitude');
+ grid on
+ legend('NAS - flight', 'NAS - simulated', 'Simulation', 'Location','southoutside', 'Orientation', 'horizontal');
+
+ figure
+ plot(res.time, res.velocity.vel_NasFlight);
+ hold on
+ plot(res.time, res.velocity.vel_NasSim);
+ plot(Ta, vecnorm(Ya(:, 4:6), 2, 2));
+ xlabel('Time [s]');
+ ylabel('Velocity [m/s]');
+ title('Total Velocity');
+ grid on
+ legend('NAS - flight', 'NAS - simulated', 'Simulation', 'Location','southoutside', 'Orientation', 'horizontal');
+
+ figure
+ plot(res.time, res.velocity.velVert_NasFlight);
+ hold on
+ plot(res.time, res.velocity.velVert_NasSim);
+ plot(Ta, -data_ascent.velocities(3, :));
+ xlabel('Time [s]');
+ ylabel('Velocity [m/s]');
+ title('Vertical Velocity');
+ grid on
+ legend('NAS - flight', 'NAS - simulated', 'Simulation', 'Location','southoutside', 'Orientation', 'horizontal');
+
+
+ figure
+ plot(res.time(indexCO:end),res.CA.CA_NasFlight(indexCO:end));
+ hold on
+ plot(res.time(indexCO:end),res.CA.CA_NasSim(indexCO:end));
+ plot(res.time(indexCO:end), [res.CA.CA0(indexCO:indexOpeningARB) res.CA.CA100(indexOpeningARB+1:end)]);
+ plot(data_ascent.integration.t, data_ascent.coeff.CA);
+ xlabel('Time [s]');
+ ylabel('CA [-]');
+ title('CA from telemetry');
+ grid on
+ legend('NAS - flight', 'NAS - simulated', 'MD', 'Simulation', 'Location','southoutside', 'Orientation', 'horizontal');
+
end
diff --git a/RoccarasoFlight/postprocessing/AFD/nasSim.mat b/RoccarasoFlight/postprocessing/AFD/nasSim.mat
new file mode 100644
index 0000000000000000000000000000000000000000..60c73ce495cf34f198d334dbb3e23d0d9b30c4ab
Binary files /dev/null and b/RoccarasoFlight/postprocessing/AFD/nasSim.mat differ
diff --git a/RoccarasoFlight/postprocessing/AFD/roccEngine.mat b/RoccarasoFlight/postprocessing/AFD/roccEngine.mat
new file mode 100644
index 0000000000000000000000000000000000000000..d5c3231a657a2d4c57092071be9c4047465c4795
Binary files /dev/null and b/RoccarasoFlight/postprocessing/AFD/roccEngine.mat differ
diff --git a/RoccarasoFlight/postprocessing/AFD/src/roccEngine.mat b/RoccarasoFlight/postprocessing/AFD/src/roccEngine.mat
new file mode 100644
index 0000000000000000000000000000000000000000..d5c3231a657a2d4c57092071be9c4047465c4795
Binary files /dev/null and b/RoccarasoFlight/postprocessing/AFD/src/roccEngine.mat differ
diff --git a/RoccarasoFlight/postprocessing/AFD/src/thrustBricolage.m b/RoccarasoFlight/postprocessing/AFD/src/thrustBricolage.m
new file mode 100644
index 0000000000000000000000000000000000000000..59bde01708f551575b9a785f4d942ca1065a4f64
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/AFD/src/thrustBricolage.m
@@ -0,0 +1,56 @@
+clear;
+clc;
+close all;
+
+%% Load data
+
+load("..\..\RPS\msaEngineData.mat"); % curva nuova
+load("..\..\RPS\cleanData\ROC2\engineFlightDataOld.mat"); % curva originale
+addpath("..\..\..\..\..\msa-toolkit\commonFunctions\miscellaneous\");
+
+%% Isolate transients
+tNodal = [0.982008, 1.25, 4.79702, 5.10202];
+timeOld = flightData.timeThrust;
+ThrustOld = flightData.Thrust;
+
+indIGN = and(timeOld > tNodal(1), timeOld < tNodal(2));
+indCO = and(timeOld > tNodal(3), timeOld < tNodal(4));
+
+timeIGN = timeOld(indIGN);
+ThrustIGN = ThrustOld(indIGN);
+
+timeCO = timeOld(indCO);
+ThrustCO = ThrustOld(indCO);
+
+%% Add transient to computed thrust
+
+timeComp = msaEngineData.timeThrustReal;
+ThrustComp = msaEngineData.ThrustReal;
+
+timeIGN = timeIGN - timeIGN(1);
+timeComp = timeComp - timeComp(1) + timeIGN(end);
+timeCO = timeCO - timeCO(1) + timeComp(end);
+
+%% Reshape transient
+
+T0 = ThrustComp(1);
+TEnd = ThrustComp(end);
+
+ThrustIGN = ThrustIGN - ThrustIGN(1); % parte da 0
+ThrustCO = ThrustCO - ThrustCO(end); % finisce a 0
+
+ThrustIGN = ThrustIGN * T0/ThrustIGN(end);
+ThrustCO = ThrustCO * TEnd/ThrustCO(1);
+
+%% Plot results
+
+time = [timeIGN(1:end-1); timeComp; timeCO(2:end)];
+Thrust = [ThrustIGN(1:end-1); ThrustComp; ThrustCO(2:end)];
+mass = reshapeVector(time, msaEngineData.time, msaEngineData.Mp);
+
+plot(time, Thrust);
+
+save('roccEngine.mat', 'time', 'Thrust', 'mass');
+
+
+
diff --git a/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m b/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m
index 49f6d72be7e681a71d6b9028e05ea7951db7cf84..c86f3356f311c2a972f5860ec61b769af0cbc109 100644
--- a/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m
+++ b/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m
@@ -56,6 +56,8 @@ if settings.saveResults
msaEngineData.mp = engineSim.mfr.mp;
msaEngineData.time = engineSim.time;
msaEngineData.Thrust = engineSim.performances.T;
+ msaEngineData.ThrustReal = engineFly.T;
+ msaEngineData.timeThrustReal = engineFly.time;
save('msaEngineData', 'msaEngineData');
end
diff --git a/RoccarasoFlight/postprocessing/RPS/msaEngineData.mat b/RoccarasoFlight/postprocessing/RPS/msaEngineData.mat
index 0cdec558af9edf9ae92a240128b897ec6ee850c5..7ad64beb50dd65246eaf943eb67c0b7183694f8f 100644
Binary files a/RoccarasoFlight/postprocessing/RPS/msaEngineData.mat and b/RoccarasoFlight/postprocessing/RPS/msaEngineData.mat differ
diff --git a/RoccarasoFlight/postprocessing/commonFunctions/atmosphereAlt.m b/RoccarasoFlight/postprocessing/commonFunctions/atmosphereAlt.m
new file mode 100644
index 0000000000000000000000000000000000000000..52764d6cb71e24875a543be72113a95d559ba391
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/commonFunctions/atmosphereAlt.m
@@ -0,0 +1,35 @@
+function absoluteAltitude = atmosphereAlt(pressure, local)
+
+
+
+atmCostants = [0.0065, 401.87434, 1.86584515, 1.225, ...
+ 101325, 288.15, 287.05]; % atmosisa constants: [a, gamma*R, R*a, rho0, p0, T0, R]
+g = 9.81;
+P0 = atmCostants(5);
+T0 = atmCostants(6);
+
+if length(local) == 4 % if T, P, rho (at local altitude) and local altitude are given as inputs, T0, P0 and rho0 are computed in order to have a more accurate atmosphere model
+ altLocal = local(1);
+ TLocal = local(2);
+ PLocal = local(3);
+ T0 = TLocal + atmCostants(1)*altLocal;
+ theta0 = TLocal/T0;
+ P0 = PLocal/(theta0^(g/atmCostants(3)));
+elseif length(local) == 3 % if T, P (at local altitude) and local altitude are given as inputs, T0, P0 and rho0 are computed in order to have a more accurate atmosphere model
+ altLocal = local(1);
+ TLocal = local(2);
+ PLocal = local(3);
+ T0 = TLocal + atmCostants(1)*altLocal;
+ theta0 = TLocal/T0;
+ P0 = PLocal/(theta0^(g/atmCostants(3)));
+elseif length(local) == 2 % if T (at a local altitude) and local altitude are given as inputs, T0, P0 and rho0 are computed in order to have a more accurate atmosphere model
+ altLocal = local(1);
+ TLocal = local(2);
+ T0 = TLocal + atmCostants(1)*altLocal;
+end
+
+theta = (pressure./P0).^(atmCostants(3)./g);
+T = T0.*theta;
+absoluteAltitude = (T0 - T)./atmCostants(1);
+
+end