diff --git a/.gitignore b/.gitignore
index 49842d54e7de3dd0c4d4498e3406b47d735d96e0..d8fa1796da6d059c5f8a368017f5ec937f1614a3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,3 @@
 *.asv
-*.avi
-
+**/*.avi
+**/*.pdf
diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000000000000000000000000000000000000..c18998f5ff69dedd742dcbe6649a9ca5c7c3d60e
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,6 @@
+[submodule "RoccarasoFlight\\postprocessing\\commonFunctions\\general-utilities"]
+	path = RoccarasoFlight\\postprocessing\\commonFunctions\\general-utilities
+	url = https://git.skywarder.eu/skyward/general-utilities
+[submodule "general-utilities"]
+	path = general-utilities
+	url = https://git.skywarder.eu/skyward/general-utilities
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/GNC/configLaunch.m b/RoccarasoFlight/postprocessing/GNC/configLaunch.m
new file mode 100644
index 0000000000000000000000000000000000000000..c9f3393eb6206211d50a09e67bb006bb368e5c1a
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/GNC/configLaunch.m
@@ -0,0 +1,8 @@
+% configuration file for roccaraso:
+
+settings.launchRail.elevation = deg2rad(79);        % [rad] launchrail elevation angle
+settings.launchRail.azimuth   = deg2rad(178);       % [rad] launchrail azimuth angle
+settings.launchRail.lat0      = 41.8084579;         % [deg] launchrail latitude
+settings.launchRail.lon0      = 14.0546408;         % [deg] launchrail longitude
+
+settings.sensorBoard.board2rocket = deg2rad([0 0 -45]);  % [deg] euler rotation angles between the sensorboard and the rocket
\ No newline at end of file
diff --git a/RoccarasoFlight/postprocessing/GNC/mainPP.m b/RoccarasoFlight/postprocessing/GNC/mainPP.m
index bf0333e7d87e61d29ee62f6991ef1681bd6b3e69..0045d05ba3fdba84368e86f70a001434c99cc7ce 100644
--- a/RoccarasoFlight/postprocessing/GNC/mainPP.m
+++ b/RoccarasoFlight/postprocessing/GNC/mainPP.m
@@ -1,12 +1,15 @@
 clear; close all; clc;
 
 addpath(genpath('../commonFunctions/'))
+addpath('../RCS/Functions&Files/')
+checkSubmodules;
+addpath('../../../general-utilities/')
 
 %% import data from logs
 
 % import and parse data
 [main,motor,COTS,RIG,GS] = importData('main','cots','motor','rig','gs');
-
+configLaunch;
 % import algorithm configuration files
 ABK.referenceTrajectory = readmatrix('softwareConfigurationFiles\ABK_references_Gemini_Roccaraso_September_2023.csv');
 ABK.closedTrajectory = ABK.referenceTrajectory(:,[1,12]);
diff --git a/RoccarasoFlight/postprocessing/GNC/postProcess.m b/RoccarasoFlight/postprocessing/GNC/postProcess.m
index a0fa8e97a93587847f0406d2ef7e31b765d0fef7..8facc9a8d00086d3b0465bfaa072ffcf00cc943c 100644
--- a/RoccarasoFlight/postprocessing/GNC/postProcess.m
+++ b/RoccarasoFlight/postprocessing/GNC/postProcess.m
@@ -1,5 +1,12 @@
 % post process information
 
+
+figure
+stairs(motor.mainServoData(:,1),motor.mainServoData(:,2),'DisplayName','motore')
+hold on;
+stairs(RIG.mainValve(:,1),RIG.mainValve(:,2),'DisplayName','RIG')
+
+
 %% engine acceleration vs pressure in chamber
 figure
 sgtitle('Engine: Pressure and Acceleration')
@@ -35,72 +42,92 @@ xlabel('Time [s]')
 
 
 %% ADA or COTS expulsion?
+textDim = 0.49;
 figure
-sgtitle('Expulsion: ADA or COTS?')
-subplot(2,1,1)
+title('ADA vs COTS expulsion zoomed')
 yyaxis left
-plot(main.ADAData(:,1),main.ADAData(:,3),'DisplayName','Altitude')
+plot(main.ADAData(:,1),main.ADAData(:,3),'DisplayName','Altitude','color','#4658A9')
 ylabel('ADA Altitude')
 hold on;
-xline(main.ADAEvents(main.ADAEvents(:,2)>5,1),'b--','ADA detection','HandleVisibility','off')
-xline(COTS.MainData{COTS.MainApogeeIndex,4},'r--','COTS detection','HandleVisibility','off')
-xline(COTS.MainData{COTS.MainDeploymentIndex,4},'r--','COTS deploy','HandleVisibility','off')
-
 ylabel('ADA altitude')
+% ylim([1300,1360])
+exportStandardizedFigure(gcf, 'Expulsion far',textDim, 'addMarkers',false,'changeColors',false,'exportPDF',false)
+set(gca,"YColor",'#4658A9')
 yyaxis right
-xline(main.EVENTData{findEventIndex(main.EVENTData,2,'FLIGHT_DPL_ALT_DETECTED',main.EVENTNames),1},'b--','FLIGHT DPL ALT DETECTED',HandleVisibility='off')
-plot(main.RCS_PRESSURE(:,1),main.RCS_PRESSURE(:,2),'DisplayName','RCS bay pressure')
+xline(main.ADAEvents(main.ADAEvents(:,2)>5,1),'--','ADA detection','HandleVisibility','off','Color','#00A86B')
+hold on;
+xline(COTS.MainData{COTS.MainApogeeIndex,4},'--','COTS detection','HandleVisibility','off','Color','#FDAF00')
+xline(main.EVENTData{findEventIndex(main.EVENTData,2,'FLIGHT_DPL_ALT_DETECTED',main.EVENTNames),1},'--','FLIGHT DPL ALT DETECTED','HandleVisibility','off','Color','#00A86B')
+xline(COTS.MainData{COTS.MainDeploymentIndex,4},'--','COTS deploy','HandleVisibility','off','Color','#FDAF00')
+plot(main.RCS_PRESSURE(:,1),main.RCS_PRESSURE(:,2),'DisplayName','RCS bay pressure','Color','#D3212D')
 legend
+set(gca,"YColor",'#D3212D')
+% xlim([18,19])
 
-subplot(2,1,2)
+exportStandardizedFigure(gcf, 'Expulsion zoom',textDim, 'addMarkers',false,'changeColors',false,'exportPDF',false)
+
+figure
+title('ADA vs COTS expulsion')
 yyaxis left
 plot(main.ADAData(:,1),main.ADAData(:,3),'DisplayName','Altitude')
 xline(COTS.MainData{COTS.MainApogeeIndex,4},'r--','COTS detection','HandleVisibility','off')
 xline(COTS.MainData{COTS.MainDeploymentIndex,4},'r--','COTS deploy','HandleVisibility','off')
-
 ylabel('ADA Altitude')
 hold on;
 xline(main.ADAEvents(main.ADAEvents(:,2)>5,1),'b--','ADA detection','HandleVisibility','off')
 xline(main.EVENTData{findEventIndex(main.EVENTData,2,'FLIGHT_DPL_ALT_DETECTED',main.EVENTNames),1},'b--','FLIGHT DPL ALT DETECTED',HandleVisibility='off')
-
 ylabel('ADA altitude')
 yyaxis right
 plot(main.IMU(:,1),vecnorm(main.IMU(:,2:4),2,2),'DisplayName','Accelerometer')
 legend
+
+
 % % % commento: sembra proprio che abbiamo espulso noi siummmmm e anche il
 % % %           deployment noi, GODO
 
 
 
 %% TARS
+textDim = 0.49;
 figure
-sgtitle('TARS')
-subplot(2,1,1)
+% sgtitle('TARS')
 yyaxis left
-plot(RIG.TARSData(:,1),RIG.TARSData(:,3),'DisplayName','TARS Pressure')
-ylabel('Pressure [bar]')
+plot(RIG.TARSData(:,1),RIG.TARSData(:,3),'DisplayName','Tank Pressure','Color','#4658A9')
+ylabel('Tank Pressure [bar]')
+set(gca,"YColor",'#4658A9')
+exportStandardizedFigure(gcf, 'Expulsion zoom',textDim, 'addMarkers',false,'changeColors',false,'exportPDF',false)
 yyaxis right
-plot(RIG.TARSData(:,1),RIG.TARSData(:,4),'DisplayName','TARS Mass')
-xline(RIG.EVENTData(findEventIndex(RIG.EVENTData,2,'MOTOR_START_TARS',main.EVENTNames),1),'k--',"START TARS")
-xline(RIG.EVENTData(findEventIndex(RIG.EVENTData,2,'TARS_PRESSURE_STABILIZED',main.EVENTNames),1),'b--',"PRESSURE STABILIZED")
+plot(RIG.TARSData(:,1),RIG.TARSData(:,4),'DisplayName','TARS Mass','Color','#D3212D')
+xline(RIG.EVENTData(findEventIndex(RIG.EVENTData,2,'MOTOR_START_TARS',main.EVENTNames),1),'k--',"START TARS",'HandleVisibility','off')
+% xline(RIG.EVENTData(findEventIndex(RIG.EVENTData,2,'TARS_PRESSURE_STABILIZED',main.EVENTNames),1),'--',"PRESSURE STABILIZED",'HandleVisibility','off','Color','#00A86B')
 title('Pressure and Mass measures')
 xlabel('Time [s]')
 ylabel('Mass [kg]')
-xlim([12900,13100])
+set(gca,"YColor",'#D3212D')
+xlim([-1160,-1000])
+legend
+exportStandardizedFigure(gcf, 'TARS pressure mass',textDim, 'addMarkers',false,'changeColors',false,'overwriteFigure',true,'exportPDF',false)
 
-subplot(2,1,2)
+figure
 yyaxis left
-plot(RIG.TARScheckMass(2:end,1),abs(diff(RIG.TARScheckMass(:,2))),'DisplayName','M tank')
+plot(RIG.TARScheckMass(2:end,1),abs(diff(RIG.TARScheckMass(:,2))),'DisplayName','OX mass difference','Color','#4658A9')
+set(gca,"YColor",'#4658A9')
 yline(0.2,'k--','TARGET OFFSET','HandleVisibility','off')
-ylabel('Mass [kg]')
+ylabel('Mass Offset [kg]')
 hold on; grid on;
+exportStandardizedFigure(gcf, 'Expulsion zoom',textDim, 'addMarkers',false,'changeColors',false,'exportPDF',false)
 yyaxis right
-plot(RIG.ventValve(:,1),RIG.ventValve(:,2),'DisplayName','Venting Valve Position')
-ylabel('Position [%]')
+stairs(RIG.ventValve(:,1),RIG.ventValve(:,2)*100,'DisplayName','Venting Valve Position','Color','#D3212D')
+xline(RIG.EVENTData(findEventIndex(RIG.EVENTData,2,'MOTOR_START_TARS',main.EVENTNames),1),'k--',"START TARS",'HandleVisibility','off')
+% xline(RIG.EVENTData(findEventIndex(RIG.EVENTData,2,'TARS_PRESSURE_STABILIZED',main.EVENTNames),1),'--',"PRESSURE STABILIZED",'HandleVisibility','off','Color','#00A86B')
+ylabel('Valve Position [\%]')
 xlabel('Time [s]')
 legend
-xlim([12900,13100])
-title('Mass finite difference 1HZ vs venting valve')
+set(gca,"YColor",'#D3212D')
+xlim([-1160,-1000])
+title('Mass finite difference vs venting valve')
+exportStandardizedFigure(gcf, 'TARS mass offset',textDim, 'addMarkers',false,'changeColors',false,'overwriteFigure',true,'exportPDF',false)
+
 % % % commento: non abbiamo lasciato che TARS finisse, ma abbiamo stoppato
 % % %           noi perché abbiamo visto fase liquida. Probabilmente si
 % % %           sarebbe stoppato a vedere l'andamento della massa prima del venting
@@ -109,7 +136,9 @@ title('Mass finite difference 1HZ vs venting valve')
 %% NAS obsw vs NAS risimulato
 % build sensors
 sensorTot.nas.time = main.NASData(:,1);
-sensorTot.nas.states = [main.NASData(1,2:end)];
+q_init = eul2quat([settings.launchRail.azimuth,settings.launchRail.elevation,0]);
+bias_init = zeros(1,3);
+sensorTot.nas.states = [main.NASData(1,2:7), q_init ,bias_init];% pos vel quat bias
 sensorTot.nas.n_old = 2;
 sensorTot.nas.P   = eye(12);
 sensorTot.nas.P(3,3) = 10;
@@ -139,6 +168,10 @@ sensorTot.pitot.speed = main.PITOT(:,3);
 %%% magnetometer calibration:
 center_offset = [+0.3, -0.3, 0.46];
 sensorTot.imu.magnetometerMeasuresCalib = sensorTot.imu.magnetometerMeasures - center_offset;
+R_b2r = quat2dcm(eul2quat(settings.sensorBoard.board2rocket));
+% R_calib = findRotationMatrix(sensorTot.imu.magnetometerMeasures(1,:),[1 0 0]);
+% 
+sensorTot.imu.magnetometerMeasuresCalib = (R_b2r*sensorTot.imu.magnetometerMeasuresCalib')';
 
 % set configuration parameters ( equal to OBSW )
 nas.dt_k          =   0.02;                                    % [s]        nas time step
@@ -260,7 +293,8 @@ sgtitle('NAS OBSW vs NAS SIMULATORE - Velocity norm')
 % velocity vs pitot
 plot(main.NASData(:,1),vecnorm(main.NASData(:,5:7),2,2),'DisplayName','OBSW')
 hold on;
-plot(sensorTot_NoPitot.nas.time,vecnorm(sensorTot_NoPitot.nas.states(:,4:6),2,2),'DisplayName','Simulatore')
+plot(sensorTot_NoPitot.nas.time,vecnorm(sensorTot_NoPitot.nas.states(:,4:6),2,2),'DisplayName','Simulator No Pitot')
+plot(sensorTot_WithPitot.nas.time,vecnorm(sensorTot_WithPitot.nas.states(:,4:6),2,2),'DisplayName','Simulator With Pitot')
 plot(main.PITOT(:,1),main.PITOT(:,3),'DisplayName','Pitot OBSW')
 legend
 xline(main.NASData(main.NASApogeeIndex,1),'--','APOGEE','HandleVisibility','off')
@@ -350,18 +384,19 @@ xline(main.NASData(main.NASApogeeIndex,1),'--','APOGEE','HandleVisibility','off'
 % % %           Mi devo inventare un modo per sincronizzarle
 
 %% MAGNETOMETER CALIBRATION
-figure
-scatter3(main.MAGNETOMETER(:,2),main.MAGNETOMETER(:,3),main.MAGNETOMETER(:,4),'.','DisplayName','Raw')
+fig_magnetometerCalibration = figure;
+% scatter3(main.MAGNETOMETER(:,2),main.MAGNETOMETER(:,3),main.MAGNETOMETER(:,4),'.','DisplayName','Raw')
 hold on;
 scatter3(main.IMU_CORRECTED(:,10),main.IMU_CORRECTED(:,11),main.IMU_CORRECTED(:,12),'.','DisplayName','Rotated')
 scatter3(sensorTot.imu.magnetometerMeasuresCalib(:,1),sensorTot.imu.magnetometerMeasuresCalib(:,2),sensorTot.imu.magnetometerMeasuresCalib(:,3),'.','DisplayName','Calibrated')
-
 axis equal
 legend
 title('Magnetometer calibration')
 xlabel('X')
 ylabel('Y')
 zlabel('Z')
+view(-60,40)
+% exportgraphics(fig_magnetometerCalibration,'MagnetometerCalibration.pdf')
 
 % % % commento: i dati raw ok, però anche i dati non raw non sono calibrati
 % % %           purtroppo, la sfera non è centrata in zero, a meno che non
@@ -485,6 +520,8 @@ title('Pressure')
 legend
 xlabel('Time [s]')
 ylabel('Pressure [Pa]')
+xline(motor.mainServoData(find(motor.mainServoData(:,2)>0.965,1),1),'k--','OPEN MAIN','HandleVisibility','off')
+xline(motor.mainServoData(find(motor.mainServoData(:,2)<0.001,1),1),'k--','CLOSE MAIN','HandleVisibility','off')
 % apogee
 subplot(3,1,3)
 plot(main.MEAControllerStatus(main.MEAControllerStatus(:,1)<t_end,1),main.MEAControllerStatus(main.MEAControllerStatus(:,1)<t_end,3),'DisplayName','OBSW')
@@ -494,36 +531,65 @@ plot(sensorTot_NoPitot.mea.time,sensorTot_NoPitot.mea.predictedApogee,'DisplayNa
 plot(sensorTot_WithPitot.mea.time,sensorTot_WithPitot.mea.predictedApogee,'DisplayName','Simulator WithPitot' )
 plot(main.NASData(main.NASData(:,1)<t_end,1),-main.NASData(main.NASData(:,1)<t_end,4),'DisplayName','NAS Altitude')
 yline(main.MEATarget,'--','SHUTDOWN TARGET','HandleVisibility','off')
+xline(motor.mainServoData(find(motor.mainServoData(:,2)>0.965,1),1),'k--','OPEN MAIN','HandleVisibility','off')
+xline(motor.mainServoData(find(motor.mainServoData(:,2)<0.001,1),1),'k--','CLOSE MAIN','HandleVisibility','off')
 legend
 title('Apogee prediction')
 xlabel('Time [s]')
 ylabel('Altitude [m]')
 
+%% report images for MEA
+figure
+plot(main.MEAControllerStatus(main.MEAControllerStatus(:,1)<t_end,1),main.MEAControllerStatus(main.MEAControllerStatus(:,1)<t_end,2),'DisplayName','OBSW')
+hold on;
+plot(sensorTot_OBSW.mea.time,sensorTot_OBSW.mea.estimatedMass,'DisplayName','Last HIL test')
+legend
+title('Estimated Mass')
+xlabel('Time [s]')
+ylabel('Mass [kg]')
+
+figure
+plot(main.MEAControllerStatus(main.MEAControllerStatus(:,1)<t_end,1),main.MEAControllerStatus(main.MEAControllerStatus(:,1)<t_end,3),'DisplayName','OBSW')
+hold on;
+plot(sensorTot_OBSW.mea.time,sensorTot_OBSW.mea.predictedApogee,'DisplayName','Last HIL test')
+legend
+title('Estimated Mass')
+xlabel('Time [s]')
+ylabel('Mass [kg]')
+
+
+
 %% ABK wrt references
+textDim = 0.6;
 main.ABKServoData(:,5) = zeros(length(main.ABKServoData(:,1)),1);
 for i = 1:length(main.ABKServoData(:,1))
     main.ABKServoData(i,5) = -main.NASData(find((main.NASData(:,1)-main.ABKServoData(i,1))>0,1,'first'),4);
 end
 figure
 yyaxis left
-plot(-main.NASData(1:main.NASApogeeIndex,4),-main.NASData(1:main.NASApogeeIndex,7),'b','DisplayName','NAS OBSW')
+plot(-main.NASData(1:main.NASApogeeIndex,4),-main.NASData(1:main.NASApogeeIndex,7),'DisplayName','NAS OBSW','Color','#4658A9')
 hold on;
-plot(-sensorTot_NoPitot.nas.states(1:main.NASApogeeIndex,3),-sensorTot_NoPitot.nas.states(1:main.NASApogeeIndex,6),'b','DisplayName','NAS Simulatore NoPitot')
-plot(-sensorTot_WithPitot.nas.states(1:main.NASApogeeIndex,3),-sensorTot_WithPitot.nas.states(1:main.NASApogeeIndex,6),'b','DisplayName','NAS Simulatore Pitot')
+% plot(-sensorTot_NoPitot.nas.states(1:main.NASApogeeIndex,3),-sensorTot_NoPitot.nas.states(1:main.NASApogeeIndex,6),'b','DisplayName','NAS Simulatore NoPitot')
+% plot(-sensorTot_WithPitot.nas.states(1:main.NASApogeeIndex,3),-sensorTot_WithPitot.nas.states(1:main.NASApogeeIndex,6),'b','DisplayName','NAS Simulatore Pitot')
 
 plot(ABK.closedTrajectory(:,1),ABK.closedTrajectory(:,2),'k--','DisplayName','Closed Trajectory')
 plot(ABK.openTrajectory(:,1),ABK.openTrajectory(:,2),'ks-','DisplayName','Open Trajectory')
 ylabel('Vertical Velocity [m/s]')
 h = gca;
-h.YColor = 'b';
+h.YColor = '#4658A9';
+exportStandardizedFigure(gcf, 'Expulsion far',textDim, 'addMarkers',false,'forcedMarkers',5,'changeColors',false,'exportPDF',false)
+
 yyaxis right
-stairs(main.ABKServoData(:,5),main.ABKServoData(:,4),'r','DisplayName','ABK')
+stairs(main.ABKServoData(:,5),main.ABKServoData(:,4),'DisplayName','ABK','Color','#D3212D')
 h = gca;
-h.YColor = 'r';
+h.YColor = '#D3212D';
 ylabel('ABK position')
 xlabel('Altitude [m]')
 legend
 title('ABK vs References')
+ylim([0,1.1])
+exportStandardizedFigure(gcf, 'ABK wrt references',textDim, 'addMarkers',false,'forcedMarkers',5,'changeColors',false,'exportPDF',false)
+
 
 % % % commento: la shadowmode probabilmente si è mangiata il fatto che il primo
 % % %           comando dovrebbe essere senza filtro, per permettere agli aerofreni di
@@ -532,4 +598,29 @@ title('ABK vs References')
 % % %           Il pitot anticipa tanto come si è visto anche nel volo di
 % % %           euroc di pyxis, il che potrebbe portare a un comportamento
 % % %           del controllore inaspettato. Sicuramente però, se è tarato è la migliore
-% % %           stima della velocità xbody che abbiamo
\ No newline at end of file
+% % %           stima della velocità xbody che abbiamo
+
+%% NAS consistency check
+
+v_int_PITOT = 0;
+v_int_NOpitot = 0;
+for i = 2:length(sensorTot_WithPitot.nas.states(:,3))
+    v_int_PITOT(i) = v_int_PITOT(i-1) + sum(sensorTot_WithPitot.nas.states([i-1,i],6))*nas.dt_k/2;
+    v_int_NOpitot(i) = v_int_NOpitot(i-1) + sum(sensorTot_NoPitot.nas.states([i-1,i],6))*nas.dt_k/2;
+end
+figure
+sgtitle('NAS consistency check')
+plot(sensorTot_WithPitot.nas.time,v_int_PITOT,'DisplayName','Int vel PITOT')
+hold on;
+plot(sensorTot_WithPitot.nas.time,sensorTot_WithPitot.nas.states(:,3),'DisplayName','Altitude PITOT')
+plot(sensorTot_NoPitot.nas.time,v_int_NOpitot,'DisplayName','Int vel NO PITOT')
+plot(sensorTot_NoPitot.nas.time,sensorTot_NoPitot.nas.states(:,3),'DisplayName','Altitude NOPITOT')
+legend
+
+
+
+%% rocket animation
+for i = 1:length(main.NASData(:,1))
+    dcm(:,:,i) = quat2dcm(sensorTot_NoPitot.nas.states(i,[10,7:9]));
+end
+% animateOrientation(dcm, sensorTot_NoPitot.nas.time, 'prova1')
\ No newline at end of file
diff --git a/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m b/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m
index 081f1fdd1e5451b9665d1ced8464a8e4c3c9ef40..c86f3356f311c2a972f5860ec61b769af0cbc109 100644
--- a/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m
+++ b/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m
@@ -39,6 +39,11 @@ engineSim = combustionSimulation(engine, settings);
 % Computing actual performances
 engineFly = performanceReconstruction(engine, engineSim);
 
+% Loading MSA simulated motor for the flight
+load("HREmotors.mat");
+engineMSA = HREmotors(59);
+clear HREmotors;
+
 %% PLOTTING
 fprintf('DATA ANALYSIS CONCLUDED: \n\n')
 fprintf('\t - Total Impulse simulated: %.2f\n', engineSim.performances.Itot)
@@ -51,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
 
@@ -78,8 +85,9 @@ if settings.plots
     plot(engineSim.time, engineSim.performances.T)
     hold on; grid minor
     plot(engineFly.time, engineFly.T)
+    plot(engineMSA.t, engineMSA.T)
     ylabel('Thrust [N]'); xlabel('Time [s]'); title('Thrust compare')
-    legend('Simulated', 'Computed')
+    legend('Simulated', 'Computed', 'MSA-SFT06')
     axis padded
 end
 
diff --git a/RoccarasoFlight/postprocessing/RPS/msaData/HREmotors.mat b/RoccarasoFlight/postprocessing/RPS/msaData/HREmotors.mat
new file mode 100644
index 0000000000000000000000000000000000000000..b7b3c10914d77aee38c5c66ed75e6f72bdf9926f
Binary files /dev/null and b/RoccarasoFlight/postprocessing/RPS/msaData/HREmotors.mat differ
diff --git a/RoccarasoFlight/postprocessing/RPS/msaEngineData.mat b/RoccarasoFlight/postprocessing/RPS/msaEngineData.mat
index 1603214baff3dd313ed2ff0b4e98cc1f8bbf7a99..7ad64beb50dd65246eaf943eb67c0b7183694f8f 100644
Binary files a/RoccarasoFlight/postprocessing/RPS/msaEngineData.mat and b/RoccarasoFlight/postprocessing/RPS/msaEngineData.mat differ
diff --git a/RoccarasoFlight/postprocessing/RPS/src/combustionSimulation.m b/RoccarasoFlight/postprocessing/RPS/src/combustionSimulation.m
index bdacbf726dbde03c0ea8d8ea92a7cf194c9893e7..1bfaeaf8c5f67537f3d6fdc6068dfdb133e9878b 100644
--- a/RoccarasoFlight/postprocessing/RPS/src/combustionSimulation.m
+++ b/RoccarasoFlight/postprocessing/RPS/src/combustionSimulation.m
@@ -59,6 +59,8 @@ g0          = settings.g0;
 %% COMBUSTION SIMULATION
 % Defining the time
 engineSim.time = engine.acquisition.cutData.timePTop - engine.acquisition.cutData.timePTop(1);
+
+% Saving the oxidizer mass flow in the out-struct
 engineSim.mfr.mox = engine.performances.mox;
 
 % Performing Marxman simulation
@@ -83,6 +85,7 @@ OF = engineSim.mfr.mox./engineSim.mfr.mf;
 OF(1) = OF(2);
 engineSim.performances.OF = OF;
 
+% Fixing convergence problems in the first transitory phase (performances are not relevant in this phase)
 ii = find(OF < 2);
 OF(ii) = 2;
 
@@ -114,7 +117,7 @@ end
 Pe = real(Pe');
 
 % Thrust coefficient
-cT = sqrt(2*y.^2./(y - 1).*(2./(y + 1)).^((y + 1)./(y - 1))).*sqrt(1 - (Pe./Pc).^((y - 1)./y)) * lambda + (Pe - Pamb)./Pc.*epsilon;
+cT = sqrt(2*y.^2./(y - 1).*(2./(y + 1)).^((y + 1)./(y - 1))).*sqrt(1 - (Pe./Pc).^((y - 1)./y)) * lambda;
 
 % Thrust
 T  = cT.*Pc*1e5*At*etacT;
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
diff --git a/RoccarasoFlight/postprocessing/commonFunctions/graphics/BaseZoom.m b/RoccarasoFlight/postprocessing/commonFunctions/graphics/BaseZoom.m
new file mode 100644
index 0000000000000000000000000000000000000000..78a039cf56f8079b68609610b51752b01a4f0d43
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/commonFunctions/graphics/BaseZoom.m
@@ -0,0 +1,885 @@
+classdef BaseZoom < handle
+    %{
+
+        Interactive Magnification of Customized Regions.
+
+        Email: iqiukp@outlook.com
+    
+        -------------------------------------------------------------
+  
+        Version 1.3.1, 24-JAN-2022
+            -- Fixed bugs when applied to logarithmic-scale coordinates. 
+
+        Version 1.3, 17-JAN-2022
+            -- Fixed minor bugs.
+            -- Added support for image class.
+
+        Version 1.2, 4-OCT-2021
+            -- Added support for interaction
+
+        Version 1.1, 1-SEP-2021
+            -- Fixed minor bugs.
+            -- Added description of parameters.   
+
+        Version 1.0, 10-JUN-2021
+            -- Magnification of Customized Regions.
+
+        -------------------------------------------------------------
+        
+        BSD 3-Clause License
+        Copyright (c) 2022, Kepeng Qiu
+        All rights reserved.
+        
+    %}
+
+    % main  properties
+    properties
+        % main handles
+        mainFigure
+        subFigure
+        mainAxes
+        subAxes
+        roi
+        rectangleZoomedZone
+
+        % image-related properties
+        uPixels
+        vPixels
+        vuRatio
+        CData_
+        Colormap_
+        imageDim
+        imagePosition = [0.1, 0.1, 0.8, 0.6]
+        imageRectangleEdgePosition
+        imageArrow
+
+        % figure-related properties
+        mappingParams
+        figureRectangleEdgePosition
+        lineDirection
+        axesPosition
+        figureArrow
+
+        % others
+        drawFunc
+        axesClass
+        axesDone = 'off'
+        rectangleDone = 'off'
+        pauseTime = 0.2
+        textDisplay = 'on'
+        axesScale
+    end
+
+    % theme of inserted axes (sub-axes)
+    properties
+        subAxesBox = 'on'
+        subAxesinsertedLineWidth = 1.2
+        subAxesTickDirection = 'in'
+        subAxesBackgroundColor = 'w'
+    end
+
+    % theme of the zoomed zone (figures)
+    properties
+        rectangleColor = 'k'
+        rectangleFaceColor = 'none'
+        rectangleFaceAlpha = 0
+        rectangleLineStyle = '-'
+        rectangleLineWidth = 1.2
+        rectangleInteractionsAllowed = 'none'
+    end
+
+    % theme of the zoomed zone (images)
+    properties
+        imageRectangleColor = [0, 114, 189]/255
+        imageRectangleFaceColor = [0, 114, 189]/255
+        imageRectangleFaceAlpha = 0.2
+        imageRectangleLineStyle = '-'
+        imageRectangleLineWidth = 2
+        imageRectangleInteractionsAllowed = 'none'
+    end
+
+    % theme of the connected lines (images)
+    properties
+        % setting of lines between arrows
+        imageConnectedLineStyle = ':'
+        imageConnectedLineColor = 'r'
+        imageConnectedLineWidth = 1.2
+        % setting of start arrow
+        imageConnectedLineStartHeadStyle = 'ellipse' % shape of start arrow
+        imageConnectedLineStartHeadLength = 3
+        imageConnectedLineStartHeadWidth = 3
+        % setting of end arrow
+        imageConnectedLineEndHeadStyle = 'cback2' % shape of ending arrow
+        imageConnectedLineEndHeadLength = 7
+        imageConnectedLineEndHeadWidth = 7
+
+    end
+
+    % theme of the connected lines (figures)
+    properties
+        % setting of lines between arrows
+        figureConnectedLineStyle = ':'
+        figureConnectedLineColor = 'k'
+        figureConnectedLineWidth = 1.2
+        % setting of start arrow
+        figureConnectedLineStartHeadStyle = 'ellipse' % shape of start arrow
+        figureConnectedLineStartHeadLength = 3
+        figureConnectedLineStartHeadWidth = 3
+        % setting of end arrow
+        figureConnectedLineEndHeadStyle = 'cback2' % shape of ending arrow
+        figureConnectedLineEndHeadLength = 7
+        figureConnectedLineEndHeadWidth = 7
+    end
+
+    % theme of the dynamic rectangle
+    properties(Constant)
+        dynamicRectFaceColor = [0, 114, 189]/255
+        dynamicRectFaceAspect = 0.2
+        dynamicRectFacAngleMarker = 's'
+        dynamicRectFacAngleMarkerSize = 12
+        dynamicRectLineWidth = 2.5
+        dynamicRectLineColor = [0, 114, 189]/255
+    end
+
+    % dynamic properties
+    properties(Dependent)
+        XLimNew
+        YLimNew
+        affinePosition
+        dynamicPosition
+        newCData_
+        newCData
+        newCMap
+        newCMap_
+    end
+
+    methods
+        function plot(obj)
+            % main steps
+            obj.checkVersion;
+            obj.mainAxes = gca;
+            obj.axesScale.XScale = obj.mainAxes.XScale;
+            obj.axesScale.YScale = obj.mainAxes.YScale;            
+            obj.mainFigure = gcf;
+
+            if size(imhandles(obj.mainAxes),1) ~= 0
+                obj.axesClass = 'image';
+                % information about the image
+                obj.CData_ = get(obj.mainAxes.Children, 'CData');
+                %
+                obj.Colormap_ = colormap(gca);
+                if size(obj.Colormap_, 1) == 64
+                    obj.Colormap_ = colormap(gcf);
+                end
+                [obj.vPixels, obj.uPixels, ~] = size(obj.CData_);
+                obj.vuRatio = obj.vPixels/obj.uPixels;
+                obj.imageDim = length(size(obj.CData_));
+            else
+                obj.axesClass = 'figure';
+            end
+
+            switch obj.axesClass
+                case 'image'
+                    obj.insertSubAxes;
+                    obj.axesDone = 'off';
+                    fprintf('Use the left mouse button to draw a rectangle\n')
+                    fprintf('for the magnification zone...\n')
+                    obj.insertRectangle;
+                    obj.rectangleDone = 'off';
+
+                case 'figure'
+                    fprintf('Use the left mouse botton to draw a rectangle\n')
+                    fprintf('for the sub-coordinate system...\n')
+                    obj.insertSubAxes;
+                    obj.axesDone = 'off';
+                    fprintf('Use the left mouse button to draw a rectangle\n')
+                    fprintf('for the magnification zone...\n')
+                    obj.insertRectangle;
+                    obj.rectangleDone = 'off';
+            end
+        end
+
+        function checkVersion(obj)
+            % check the MATLAB version
+            version_ = version('-release');
+            year_ = str2double(version_(1:4));
+            if year_ < 2014 || (year_ == 2014 && version_(5) == 'a')
+                error('ZoomPlot V1.2 is not compatible with the versions lower than R2014b.')
+            end
+
+            if year_ >= 2017
+                set(findobj(gcf, 'type', 'Legend'), 'AutoUpdate', 'off');
+            end
+
+            if year_ > 2018 || (year_ == 2018 && version_(5) == 'b')
+                obj.drawFunc = 'drawrectangle';
+            else
+                obj.drawFunc = 'imrect';
+            end
+        end
+
+        function insertSubAxes(obj)
+            % insert an axes
+            switch obj.axesClass
+                case 'image'  %
+                    % close(obj.mainFigure)
+                    obj.subFigure = figure;
+                    obj.imagePosition(4) = obj.imagePosition(3)*obj.vuRatio;
+                    set(obj.subFigure, 'Units', 'Normalized', 'OuterPosition', obj.imagePosition);
+                    %
+                    subplot(1, 2, 1, 'Parent', obj.subFigure);
+                    image(obj.CData_);
+                    obj.mainAxes = gca;
+                    if obj.imageDim == 2
+                        colormap(obj.mainAxes, obj.Colormap_);
+                    end
+                    axis off
+                    subplot(1, 2, 2, 'Parent', obj.subFigure);
+                    image((ones(obj.vPixels, obj.uPixels)));
+                    obj.subAxes = gca;
+                    colormap(obj.subAxes, [240, 240, 240]/255);
+                    axis off
+
+                case 'figure' %
+                    switch obj.drawFunc
+                        case 'drawrectangle'
+                            obj.roi = drawrectangle(obj.mainAxes);
+                            obj.setTheme;
+                            obj.creatSubAxes;
+                            addlistener(obj.roi, 'MovingROI', @obj.allEventsForSubAxesNew);
+                            addlistener(obj.roi, 'ROIMoved', @obj.allEventsForSubAxesNew);
+                            set(gcf, 'WindowButtonDownFcn', @obj.clickEventsForSubAxes)
+                            while strcmp(obj.axesDone, 'off')
+                                pause(obj.pauseTime)
+                            end
+
+                        case 'imrect'
+                            obj.roi = imrect;
+                            obj.setTheme;
+                            func_ = makeConstrainToRectFcn('imrect',...
+                                get(obj.mainAxes, 'XLim'), get(obj.mainAxes, 'YLim'));
+                            setPositionConstraintFcn(obj.roi, func_);
+                            obj.creatSubAxes;
+                            addNewPositionCallback(obj.roi, @obj.allEventsForSubAxesOld);
+                            set(gcf, 'WindowButtonDownFcn', @obj.clickEventsForSubAxes);
+                            wait(obj.roi);
+                            while strcmp(obj.axesDone, 'off')
+                                pause(obj.pauseTime)
+                            end
+                    end
+            end
+
+        end
+
+        function insertRectangle(obj)
+            % insert an rectangle
+            switch obj.axesClass
+                case 'image' %
+                    switch obj.drawFunc
+                        case 'drawrectangle'
+                            obj.roi = drawrectangle(obj.mainAxes);
+                            obj.setTheme;
+                            obj.creatSubAxes;
+                            if strcmp(obj.subAxesBox, 'on')
+                                obj.connectAxesAndBox;
+                            end
+                            set(gcf, 'WindowButtonDownFcn', @obj.clickEventForRectangle);
+                            addlistener(obj.roi, 'MovingROI', @obj.allEventsForRectangleNew);
+                            addlistener(obj.roi, 'ROIMoved', @obj.allEventsForRectangleNew);
+                            while strcmp(obj.rectangleDone, 'off')
+                                pause(obj.pauseTime)
+                            end
+
+                        case 'imrect'
+                            obj.roi = imrect(obj.mainAxes);
+                            obj.setTheme;
+                            obj.creatSubAxes;
+                            func_ = makeConstrainToRectFcn('imrect',...
+                                get(obj.mainAxes, 'XLim'), get(obj.mainAxes, 'YLim'));
+                            setPositionConstraintFcn(obj.roi, func_);
+                            if strcmp(obj.subAxesBox, 'on')
+                                obj.connectAxesAndBox;
+                            end
+                            addNewPositionCallback(obj.roi, @obj.allEventsForRectangleOld);
+                            set(gcf, 'WindowButtonDownFcn', @obj.clickEventForRectangle);
+                            wait(obj.roi);
+                            while strcmp(obj.rectangleDone, 'off')
+                                pause(obj.pauseTime)
+                            end
+                    end
+                    %
+                    for iArrow = 1:length(obj.imageArrow)
+                        obj.imageArrow{iArrow}.Tag = 'ZoomPlot';
+                    end
+
+                case 'figure' %
+                    switch obj.drawFunc
+                        case 'drawrectangle'
+                            obj.roi = drawrectangle(obj.mainAxes);
+                            obj.setTheme;
+                            if strcmp(obj.subAxesBox, 'on')
+                                obj.connectAxesAndBox;
+                            end
+                            set(obj.subAxes, 'XLim', obj.XLimNew, 'YLim', obj.YLimNew);
+                            set(gcf, 'WindowButtonDownFcn', @obj.clickEventForRectangle);
+                            addlistener(obj.roi, 'MovingROI', @obj.allEventsForRectangleNew);
+                            addlistener(obj.roi, 'ROIMoved', @obj.allEventsForRectangleNew);
+                            while strcmp(obj.rectangleDone, 'off')
+                                pause(obj.pauseTime)
+                            end
+
+                        case 'imrect'
+                            obj.roi = imrect;
+                            obj.setTheme;
+                            func_ = makeConstrainToRectFcn('imrect',...
+                                get(obj.mainAxes, 'XLim'), get(obj.mainAxes, 'YLim'));
+                            setPositionConstraintFcn(obj.roi, func_);
+                            if strcmp(obj.subAxesBox, 'on')
+                                obj.connectAxesAndBox;
+                            end
+                            set(obj.subAxes, 'XLim', obj.XLimNew, 'YLim', obj.YLimNew);
+                            addNewPositionCallback(obj.roi, @obj.allEventsForRectangleOld);
+                            set(gcf, 'WindowButtonDownFcn', @obj.clickEventForRectangle);
+                            wait(obj.roi);
+                            while strcmp(obj.rectangleDone, 'off')
+                                pause(obj.pauseTime)
+                            end
+                    end
+                    %
+                    for iArrow = 1:length(obj.figureArrow)
+                        obj.figureArrow{iArrow}.Tag = 'ZoomPlot';
+                    end
+
+            end
+        end
+
+        function allEventsForSubAxesOld(obj, ~)
+            % callback funcion for inserted subAxes when using 'imrect'
+            if strcmp(obj.textDisplay, 'on')
+                fprintf('adjust the inserted subAxes...\n')
+            end
+            delete(obj.subAxes);
+            obj.creatSubAxes;
+            obj.subAxes.Color = obj.subAxesBackgroundColor;
+        end
+
+        function allEventsForSubAxesNew(obj, ~, evt)
+            % callback funcion for inserted subAxes when using 'drawrectangle'
+            eventName = evt.EventName;
+            if ismember(eventName, {'MovingROI', 'ROIMoved'})
+                if strcmp(obj.textDisplay, 'on')
+                    fprintf('adjust the inserted subAxes...\n')
+                end
+                delete(obj.subAxes);
+                obj.creatSubAxes;
+                obj.subAxes.Color = obj.subAxesBackgroundColor;
+            end
+        end
+
+        function clickEventsForSubAxes(obj, ~, ~)
+            % callback funcion for inserted subAxes
+            switch get(gcf, 'SelectionType')
+                % right-click
+                case 'alt'
+                    obj.axesDone = 'on';
+                    set(obj.subAxes, 'Visible', 'on');
+                    set(gcf, 'WindowButtonDownFcn', []);
+                    if strcmp(obj.textDisplay, 'on')
+                        fprintf('Inserted subAxes adjustment is done.\n\n')
+                    end
+                    delete(obj.roi);
+                    obj.subAxes.Color = obj.subAxesBackgroundColor;
+                    % left-click
+
+                case 'normal'
+                    obj.axesDone = 'off';
+                    if strcmp(obj.textDisplay, 'on')
+                        fprintf('Right-click to stop adjusting.\n')
+                    end
+                    obj.subAxes.Color = obj.subAxesBackgroundColor;
+
+                otherwise
+                    obj.axesDone = 'off';
+                    if strcmp(obj.textDisplay, 'on')
+                        fprintf('Right-click to stop adjusting.\n')
+                    end
+                    obj.subAxes.Color = obj.subAxesBackgroundColor;
+            end
+        end
+
+        function allEventsForRectangleOld(obj, ~)
+            % callback funcion for inserted rectangle when using 'imrect'
+            fprintf('adjust the inserted rectangle...\n')
+            delete(findall(gcf, 'Tag', 'ZoomPlot_'))
+            if strcmp(obj.subAxesBox, 'on')
+                obj.connectAxesAndBox;
+            end
+
+            switch obj.axesClass
+                case 'image'
+                    obj.creatSubAxes;
+                case 'figure'
+                    set(obj.subAxes, 'XLim', obj.XLimNew, 'YLim', obj.YLimNew);
+            end
+
+        end
+
+        function allEventsForRectangleNew(obj, ~, evt)
+            % callback funcion for inserted rectangle when using 'drawrectangle'
+            eventName = evt.EventName;
+            switch obj.axesClass
+                case 'image' %
+                    obj.creatSubAxes;
+                    if ismember(eventName, {'MovingROI', 'ROIMoved'})
+                        if strcmp(obj.textDisplay, 'on')
+                            fprintf('adjust the zoomed zone...\n')
+                        end
+                        delete(findall(gcf, 'Tag', 'ZoomPlot_'))
+                        if strcmp(obj.subAxesBox, 'on')
+                            obj.connectAxesAndBox;
+                        end
+                    end
+
+                case 'figure' %
+                    if ismember(eventName, {'MovingROI', 'ROIMoved'})
+                        if strcmp(obj.textDisplay, 'on')
+                            fprintf('adjust the zoomed zone...\n')
+                        end
+                        delete(findall(gcf, 'Tag', 'ZoomPlot_'))
+                        if strcmp(obj.subAxesBox, 'on')
+                            obj.connectAxesAndBox;
+                        end
+                        set(obj.subAxes,'XLim', obj.XLimNew, 'YLim', obj.YLimNew);
+                    end
+            end
+        end
+
+        function clickEventForRectangle(obj, ~, ~)
+            % callback funcion for inserted rectangle
+            switch get(gcf, 'SelectionType')
+                % right-click
+                case 'alt'
+                    obj.rectangleDone = 'on';
+                    obj.creatRectangle;
+                    set(gcf, 'WindowButtonDownFcn', []);
+                    delete(obj.roi);
+                    if strcmp(obj.textDisplay, 'on')
+                        fprintf('Inserted rectangle adjustment is done.\n\n')
+                    end
+
+                    % left-click
+                case 'normal'
+                    obj.rectangleDone = 'off';
+                    if strcmp(obj.textDisplay, 'on')
+                        fprintf('Right-click to stop adjusting.\n')
+                    end
+
+                otherwise
+                    obj.rectangleDone = 'off';
+                    if strcmp(obj.textDisplay, 'on')
+                        fprintf('Right-click to stop adjusting.\n')
+                    end
+            end
+        end
+
+        function creatSubAxes(obj)
+            % creat sub-axes
+            switch obj.axesClass
+                case 'image' %
+                    set(obj.subAxes.Children, 'CData', obj.newCData);
+                    if obj.imageDim == 2
+                        colormap(obj.subAxes, obj.newCMap)
+                    end
+
+                case 'figure'
+                    obj.subAxes = axes('Position', obj.affinePosition,...
+                                       'XScale', obj.axesScale.XScale,...
+                                       'YScale', obj.axesScale.YScale);
+                    children_ = get(obj.mainAxes, 'children');
+                    numChildren_ = 1:length(children_);
+                    for ii = 1:length(children_)
+                        if strcmp(children_(ii, 1).Type, 'images.roi.rectangle') ||...
+                                strcmp(children_(ii, 1).Type, 'hggroup')
+                            numChildren_(ii) = [];
+                        end
+                    end
+                    copyobj(children_(numChildren_), obj.subAxes);
+                    hold(obj.subAxes, 'on');
+                    set(obj.subAxes, 'LineWidth', obj.subAxesinsertedLineWidth,...
+                        'TickDir', obj.subAxesTickDirection,...
+                        'Box', obj.subAxesBox,...
+                        'Color', obj.subAxesBackgroundColor,...
+                        'XLim', get(obj.mainAxes, 'XLim'),...
+                        'YLim', get(obj.mainAxes, 'YLim'),...
+                        'XScale', obj.axesScale.XScale,...
+                        'YScale', obj.axesScale.YScale);
+                    set(obj.subAxes, 'Visible', 'off');
+            end
+        end
+
+        function creatRectangle(obj)
+            % creat rectangle
+            switch obj.axesClass
+                case 'image' %
+                    obj.rectangleZoomedZone = annotation( ...
+                        'rectangle', obj.imageRectangleEdgePosition, ...
+                        'LineWidth', obj.imageRectangleLineWidth,...
+                        'LineStyle', obj.imageRectangleLineStyle,...
+                        'FaceAlpha', obj.imageRectangleFaceAlpha,...
+                        'FaceColor', obj.imageRectangleFaceColor,...
+                        'Color', obj.imageRectangleColor);
+
+                case 'figure'
+                    obj.rectangleZoomedZone = annotation( ...
+                        'rectangle', obj.affinePosition, ...
+                        'LineWidth', obj.rectangleLineWidth,...
+                        'LineStyle', obj.rectangleLineStyle,...
+                        'FaceAlpha', obj.rectangleFaceAlpha,...
+                        'FaceColor', obj.rectangleFaceColor,...
+                        'Color', obj.rectangleColor);
+            end
+        end
+
+        function mappingParams = computeMappingParams(obj)
+            % compute the mapping parameters
+
+            switch obj.axesScale.XScale
+                case 'linear'
+                    rangeXLim = obj.mainAxes.XLim(1, 2)-obj.mainAxes.XLim(1, 1);
+                case 'log'
+                    rangeXLim = log10(obj.mainAxes.XLim(1, 2))-log10(obj.mainAxes.XLim(1, 1));
+            end
+            map_k_x = rangeXLim/obj.mainAxes.Position(3);
+
+            switch obj.axesScale.YScale
+                case 'linear'
+                    rangeYLim = obj.mainAxes.YLim(1, 2)-obj.mainAxes.YLim(1, 1);
+                case 'log'
+                    rangeYLim = log10(obj.mainAxes.YLim(1, 2))-log10(obj.mainAxes.YLim(1, 1));
+            end
+            map_k_y = rangeYLim/obj.mainAxes.Position(4);
+
+            switch obj.axesScale.XScale
+                case 'linear'
+                    map_b_x = obj.mainAxes.XLim(1)-obj.mainAxes.Position(1)*map_k_x;
+                case 'log'
+                    map_b_x = log10(obj.mainAxes.XLim(1))-obj.mainAxes.Position(1)*map_k_x;
+            end
+
+            switch obj.axesScale.YScale
+                case 'linear'
+                    map_b_y = obj.mainAxes.YLim(1)-obj.mainAxes.Position(2)*map_k_y;
+                case 'log'
+                    map_b_y = log10(obj.mainAxes.YLim(1))-obj.mainAxes.Position(2)*map_k_y;
+            end
+            mappingParams = [map_k_x, map_b_x; map_k_y, map_b_y];
+        end
+
+        function connectAxesAndBox(obj)
+            % insert lines between the inserted axes and rectangle
+
+            %   Rectangle        subAxes
+            %    2----1          2----1
+            %    3----4          3----4
+
+            switch obj.axesClass
+                case 'image' %
+                    uPixelsAll = obj.uPixels/obj.mainAxes.Position(3);
+                    vPixelsAll = obj.vPixels/obj.mainAxes.Position(4);
+                    switch obj.drawFunc
+                        case 'drawrectangle'
+                            Position_ = obj.roi.Position;
+
+                        case 'imrect'
+                            Position_ = getPosition(obj.roi);
+                    end
+
+                    obj.imageRectangleEdgePosition(1) = Position_(1)/uPixelsAll+obj.mainAxes.Position(1);
+                    obj.imageRectangleEdgePosition(2) = (obj.vPixels-Position_(2)-Position_(4))/...
+                        vPixelsAll+obj.subAxes.Position(2);
+                    obj.imageRectangleEdgePosition(3) = Position_(3)/uPixelsAll;
+                    obj.imageRectangleEdgePosition(4) = Position_(4)/vPixelsAll;
+
+                    % annotation position 1
+                    annotationPosX_1(1) = obj.imageRectangleEdgePosition(1)+obj.imageRectangleEdgePosition(3);
+                    annotationPosX_1(2) = obj.subAxes.Position(1);
+                    annotationPosY_1(1) = obj.imageRectangleEdgePosition(2);
+                    annotationPosY_1(2) = obj.subAxes.Position(2);
+                    obj.imageArrow{1} = annotation(gcf, 'doublearrow',...
+                        annotationPosX_1, annotationPosY_1,...
+                            'Color', obj.imageConnectedLineColor,...
+                            'LineWidth', obj.imageConnectedLineWidth,...
+                            'LineStyle', obj.imageConnectedLineStyle,...
+                            'Head1Style', obj.imageConnectedLineStartHeadStyle,...
+                            'Head1Length', obj.imageConnectedLineStartHeadLength,...
+                            'Head1Width', obj.imageConnectedLineStartHeadWidth,...
+                            'Head2Style', obj.imageConnectedLineEndHeadStyle,...
+                            'Head2Length', obj.imageConnectedLineEndHeadLength,...
+                            'Head2Width', obj.imageConnectedLineEndHeadWidth,...
+                            'Tag', 'ZoomPlot_');
+
+                    % annotation position 2
+                    annotationPosX_2(1) = obj.imageRectangleEdgePosition(1)+obj.imageRectangleEdgePosition(3);
+                    annotationPosX_2(2) = obj.subAxes.Position(1);
+                    annotationPosY_2(1) = obj.imageRectangleEdgePosition(2)+obj.imageRectangleEdgePosition(4);
+                    annotationPosY_2(2) = obj.subAxes.Position(2)+obj.subAxes.Position(4);
+                    obj.imageArrow{2} = annotation(gcf, 'doublearrow',...
+                        annotationPosX_2, annotationPosY_2,...
+                            'Color', obj.imageConnectedLineColor,...
+                            'LineWidth', obj.imageConnectedLineWidth,...
+                            'LineStyle', obj.imageConnectedLineStyle,...
+                            'Head1Style', obj.imageConnectedLineStartHeadStyle,...
+                            'Head1Length', obj.imageConnectedLineStartHeadLength,...
+                            'Head1Width', obj.imageConnectedLineStartHeadWidth,...
+                            'Head2Style', obj.imageConnectedLineEndHeadStyle,...
+                            'Head2Length', obj.imageConnectedLineEndHeadLength,...
+                            'Head2Width', obj.imageConnectedLineEndHeadWidth,...
+                            'Tag', 'ZoomPlot_');
+
+                case 'figure'
+                    % real coordinates of the inserted rectangle and axes
+                    obj.getAxesAndBoxPosition;
+                    % get the line direction
+                    obj.getLineDirection;
+                    % insert lines
+                    numLine = size(obj.lineDirection, 1);
+                    for i = 1:numLine
+                        tmp1 = [obj.figureRectangleEdgePosition(obj.lineDirection(i, 1), 1),...
+                            obj.figureRectangleEdgePosition(obj.lineDirection(i, 1), 2)];
+                        tmp2 = [obj.axesPosition(obj.lineDirection(i, 2), 1),...
+                            obj.axesPosition(obj.lineDirection(i, 2), 2)];
+                        pos1 = obj.transformCoordinate(tmp1, 'a2n');
+                        pos2 = obj.transformCoordinate(tmp2, 'a2n');
+                        obj.figureArrow{i} = annotation(gcf, 'doublearrow',...
+                            [pos1(1, 1), pos2(1, 1)], [pos1(1, 2), pos2(1, 2)],...
+                            'Color', obj.figureConnectedLineColor,...
+                            'LineWidth', obj.figureConnectedLineWidth,...
+                            'LineStyle', obj.figureConnectedLineStyle,...
+                            'Head1Style', obj.figureConnectedLineStartHeadStyle,...
+                            'Head1Length', obj.figureConnectedLineStartHeadLength,...
+                            'Head1Width', obj.figureConnectedLineStartHeadWidth,...
+                            'Head2Style', obj.figureConnectedLineEndHeadStyle,...
+                            'Head2Length', obj.figureConnectedLineEndHeadLength,...
+                            'Head2Width', obj.figureConnectedLineEndHeadWidth,...
+                            'Tag', 'ZoomPlot_');
+                    end
+            end
+        end
+
+        function getAxesAndBoxPosition(obj)
+            % real coordinates of the inserted rectangle
+            box1_1 = [obj.XLimNew(1, 2), obj.YLimNew(1, 2)];
+            box1_2 = [obj.XLimNew(1, 1), obj.YLimNew(1, 2)];
+            box1_3 = [obj.XLimNew(1, 1), obj.YLimNew(1, 1)];
+            box1_4 = [obj.XLimNew(1, 2), obj.YLimNew(1, 1)];
+            box1 = [box1_1; box1_2; box1_3; box1_4];
+            % real coordinates of the inserted axes
+            tmp1 = [obj.subAxes.Position(1)+obj.subAxes.Position(3),...
+                obj.subAxes.Position(2)+obj.subAxes.Position(4)];
+            box2_1 = obj.transformCoordinate(tmp1, 'n2a');
+            tmp2 = [obj.subAxes.Position(1),...
+                obj.subAxes.Position(2)+obj.subAxes.Position(4)];
+            box2_2 = obj.transformCoordinate(tmp2, 'n2a');
+            tmp3 = [obj.subAxes.Position(1), obj.subAxes.Position(2)];
+            box2_3 = obj.transformCoordinate(tmp3, 'n2a');
+            tmp4 = [obj.subAxes.Position(1)+obj.subAxes.Position(3),...
+                obj.subAxes.Position(2)];
+            box2_4 = obj.transformCoordinate(tmp4, 'n2a');
+            box2 = [box2_1; box2_2; box2_3; box2_4];
+            obj.figureRectangleEdgePosition = box1;
+            obj.axesPosition = box2;
+        end
+
+        function getLineDirection(obj)
+            % get the line direction
+            % left-upper
+            if (obj.figureRectangleEdgePosition(4, 1) < obj.axesPosition(1, 1) &&...
+                    obj.figureRectangleEdgePosition(4, 2) > obj.axesPosition(2, 2))
+                obj.lineDirection = [3, 3; 1, 1];
+            end
+            % middle-upper
+            if (obj.figureRectangleEdgePosition(4, 1) > obj.axesPosition(2, 1) &&...
+                    obj.figureRectangleEdgePosition(4, 2) > obj.axesPosition(2, 2)) &&...
+                    obj.figureRectangleEdgePosition(3, 1) < obj.axesPosition(1, 1)
+                obj.lineDirection = [4, 1; 3, 2];
+            end
+            % right-upper
+            if (obj.figureRectangleEdgePosition(3, 1) > obj.axesPosition(1, 1) &&...
+                    obj.figureRectangleEdgePosition(3, 2) > obj.axesPosition(1, 2))
+                obj.lineDirection = [2, 2; 4, 4];
+            end
+            % right-middle
+            if (obj.figureRectangleEdgePosition(3, 1) > obj.axesPosition(1, 1) &&...
+                    obj.figureRectangleEdgePosition(3, 2) < obj.axesPosition(1, 2)) &&...
+                    obj.figureRectangleEdgePosition(2, 2) > obj.axesPosition(4, 2)
+                obj.lineDirection = [2, 1; 3, 4];
+            end
+            % right-down
+            if (obj.figureRectangleEdgePosition(2, 1) > obj.axesPosition(4, 1) &&...
+                    obj.figureRectangleEdgePosition(2, 2) < obj.axesPosition(4, 2))
+                obj.lineDirection = [1, 1; 3, 3];
+            end
+            % down-middle
+            if (obj.figureRectangleEdgePosition(1, 1) > obj.axesPosition(3, 1) &&...
+                    obj.figureRectangleEdgePosition(1, 2) < obj.axesPosition(3, 2) &&...
+                    obj.figureRectangleEdgePosition(2, 1) < obj.axesPosition(4, 1))
+                obj.lineDirection = [2, 3; 1, 4];
+            end
+            % left-down
+            if (obj.figureRectangleEdgePosition(1, 1) < obj.axesPosition(3, 1) &&...
+                    obj.figureRectangleEdgePosition(1, 2) < obj.axesPosition(3, 2))
+                obj.lineDirection = [2, 2; 4, 4];
+            end
+            % left-middle
+            if (obj.figureRectangleEdgePosition(4, 1) <obj.axesPosition(2, 1) &&...
+                    obj.figureRectangleEdgePosition(4, 2) < obj.axesPosition(2, 2)) &&...
+                    obj.figureRectangleEdgePosition(1, 2) > obj.axesPosition(3, 2)
+                obj.lineDirection = [1, 2; 4, 3];
+            end
+        end
+
+        function setTheme(obj)
+            % set the theme of the dynamic rectangle
+            switch obj.drawFunc
+                case 'drawrectangle'
+                    try
+                        obj.roi.MarkerSize = obj.dynamicRectFacAngleMarkerSize;
+                    catch
+
+                    end
+                    obj.roi.Color = obj.dynamicRectFaceColor;
+                    obj.roi.FaceAlpha = obj.dynamicRectFaceAspect;
+                    obj.roi.LineWidth = obj.dynamicRectLineWidth;
+                case 'imrect'
+                    children_ = get(findobj(gca, 'type', 'hggroup'), 'children');
+                    % 8 angles
+                    for i = [1:4, 6:2:12]
+                        children_(i).LineWidth = obj.dynamicRectLineWidth*0.6;
+                        children_(i).Color = obj.dynamicRectLineColor;
+                        children_(i).MarkerSize = obj.dynamicRectFacAngleMarkerSize;
+                        children_(i).Marker = obj.dynamicRectFacAngleMarker;
+                        children_(i).MarkerEdgeColor = 'k';
+                        children_(i).MarkerFaceColor = obj.dynamicRectFaceColor;
+                    end
+                    % 4 lines
+                    for i = 5:2:11
+                        children_(i).Color = obj.dynamicRectFaceColor;
+                        children_(i).LineWidth = obj.dynamicRectLineWidth;
+                        children_(i).Marker = 'none';
+                    end
+                    % dynamic rectangle
+                    children_(13).FaceAlpha = obj.dynamicRectFaceAspect;
+                    children_(13).FaceColor = obj.dynamicRectFaceColor;
+            end
+        end
+
+        function coordinate = transformCoordinate(obj, coordinate, type)
+            % coordinate transformation
+            switch type
+                % absolute coordinates to normalized coordinates
+                case 'a2n'
+                    switch obj.axesScale.XScale
+                        case 'linear'
+                            coordinate(1, 1) = (coordinate(1, 1)-obj.mappingParams(1, 2))...
+                                /obj.mappingParams(1, 1);
+                        case 'log'
+                            coordinate(1, 1) = (log10(coordinate(1, 1))-obj.mappingParams(1, 2))...
+                                /obj.mappingParams(1, 1);
+                    end
+
+                    switch obj.axesScale.YScale
+                        case 'linear'
+                            coordinate(1, 2) = (coordinate(1, 2)-obj.mappingParams(2, 2))...
+                                /obj.mappingParams(2, 1);
+                        case 'log'
+                            coordinate(1, 2) = (log10(coordinate(1, 2))-obj.mappingParams(2, 2))...
+                                /obj.mappingParams(2, 1);
+                    end
+
+                % normalized coordinates to absolute coordinates
+                case 'n2a'
+                    switch obj.axesScale.XScale
+                        case 'linear'
+                            coordinate(1, 1) = coordinate(1, 1)*obj.mappingParams(1, 1)...
+                                +obj.mappingParams(1, 2);
+                        case 'log'
+                            coordinate(1, 1) = 10^(coordinate(1, 1)*obj.mappingParams(1, 1)...
+                                +obj.mappingParams(1, 2));
+                    end
+
+                    switch obj.axesScale.YScale
+                        case 'linear'
+                            coordinate(1, 2) = coordinate(1, 2)*obj.mappingParams(2, 1)...
+                                +obj.mappingParams(2, 2);
+                        case 'log'
+                            coordinate(1, 2) = 10^(coordinate(1, 2)*obj.mappingParams(2, 1)...
+                                +obj.mappingParams(2, 2));
+                    end
+            end
+        end
+
+        % dependent properties
+        function dynamicPosition = get.dynamicPosition(obj)
+            switch obj.drawFunc
+                case 'drawrectangle'
+                    dynamicPosition = obj.roi.Position;
+                case 'imrect'
+                    dynamicPosition = getPosition(obj.roi);
+            end
+        end
+
+        % dependent properties
+        function XLimNew = get.XLimNew(obj)
+            XLimNew = [obj.dynamicPosition(1), obj.dynamicPosition(1)+obj.dynamicPosition(3)];
+        end
+
+        % dependent properties
+        function YLimNew = get.YLimNew(obj)
+            YLimNew = [obj.dynamicPosition(2), obj.dynamicPosition(2)+obj.dynamicPosition(4)];
+        end
+
+        % dependent properties
+        function affinePosition = get.affinePosition(obj)
+            obj.mappingParams = obj.computeMappingParams;
+            tmp1 = obj.transformCoordinate([obj.XLimNew(1, 1), obj.YLimNew(1, 1)], 'a2n');
+            tmp2 = obj.transformCoordinate([obj.XLimNew(1, 2), obj.YLimNew(1, 2)], 'a2n');
+            affinePosition(1, 1) = tmp1(1, 1);
+            affinePosition(1, 2) = tmp1(1, 2);
+            affinePosition(1, 3) = tmp2(1, 1)-tmp1(1, 1);
+            affinePosition(1, 4) = tmp2(1, 2)-tmp1(1, 2);
+        end
+
+        % dependent properties
+        function newCData_ = get.newCData_(obj)
+            switch obj.drawFunc
+                case 'drawrectangle'
+                    Position_ = obj.roi.Position;
+                case 'imrect'
+                    Position_ = getPosition(obj.roi);
+            end
+            newCData_ = imcrop(obj.CData_,obj.Colormap_, Position_);
+        end
+
+        % dependent properties
+        function newCData = get.newCData(obj)
+            switch obj.imageDim
+                case 2
+                    [newCData, ~] = imresize(obj.newCData_, obj.Colormap_, [obj.vPixels, obj.uPixels]);
+                    %  [~, newCMap] = imresize(obj.newCData_, obj.newCMap_, [obj.vPixels, obj.uPixels]);
+                case 3
+                    newCData = imresize(obj.newCData_, [obj.vPixels, obj.uPixels]);
+            end
+        end
+
+        % dependent properties
+        function newCMap = get.newCMap(obj)
+            switch obj.imageDim
+                case 2
+                    [~, newCMap] = imresize(obj.newCData_, obj.Colormap_, [obj.vPixels, obj.uPixels]);
+                case 3
+                    newCMap=[];
+            end
+        end
+    end
+end
+
diff --git a/RoccarasoFlight/postprocessing/commonFunctions/miscellaneous/checkSubmodules.m b/RoccarasoFlight/postprocessing/commonFunctions/miscellaneous/checkSubmodules.m
new file mode 100644
index 0000000000000000000000000000000000000000..f7d075e9ca96251b95c277ccb3a79cf7efe35e9a
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/commonFunctions/miscellaneous/checkSubmodules.m
@@ -0,0 +1,25 @@
+currentPath = pwd;
+localRepoPath_generalUtilities = '../../../general-utilities';
+if ~exist([localRepoPath_generalUtilities,'/exportStandardizedFigure'],'file')
+    % clone repo
+    answer = input('WARNING! You don''have the exportStandardizedFigure installed. Do you want to install it? (recommended) (y/n)','s');
+    if answer == "y" || answer == "yes"
+        
+        cloneType = input('Which type of clone do you want to use? (SSH/HTTP) ','s');
+        if cloneType =="SSH" || cloneType == "ssh"
+            mkdir(localRepoPath_generalUtilities)    
+            system(['git submodule add "git@git.skywarder.eu:skyward/general-utilities.git" ',localRepoPath_generalUtilities])
+        elseif cloneType == "HTTP" || cloneType == "http"
+            mkdir(localRepoPath_generalUtilities)
+            system(['git submodule add "https://git.skywarder.eu/skyward/general-utilities.git" ',localRepoPath_generalUtilities])
+        else
+            fprintf('\n WARNING! Input not valid: aborting')
+        return
+        end
+    end
+else 
+    answer_pull = input('Do you want to update the graphics submodule? (y/n)','s');
+    if answer_pull == "yes" || answer_pull == "y"
+        system('git submodule update --recursive')
+    end
+end
diff --git a/RoccarasoFlight/postprocessing/commonFunctions/miscellaneous/findRotationMatrix.m b/RoccarasoFlight/postprocessing/commonFunctions/miscellaneous/findRotationMatrix.m
new file mode 100644
index 0000000000000000000000000000000000000000..93e196757e7f3a1cfa0423773047e5f9ee5b94a9
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/commonFunctions/miscellaneous/findRotationMatrix.m
@@ -0,0 +1,12 @@
+function R = findRotationMatrix(a,b)
+
+% finds the rotation matrix to apply to rotate a into b
+v = cross(a,b);
+s = norm(v);
+c = dot(a,b);
+
+if size(v,2)> size(v,1)
+    v = v';
+end
+skew = cross([v,v,v],eye(length(v)));
+R = eye(length(a)) + skew + skew^2 * (1-c)/s^2;
\ No newline at end of file
diff --git a/RoccarasoFlight/postprocessing/commonFunctions/parseData.m b/RoccarasoFlight/postprocessing/commonFunctions/parseData.m
index 620cdd499246542c1339eb794d903ac3cf30cf70..e7f558f78d6e5a0dc8bf15789fdd88a58a07525f 100644
--- a/RoccarasoFlight/postprocessing/commonFunctions/parseData.m
+++ b/RoccarasoFlight/postprocessing/commonFunctions/parseData.m
@@ -160,11 +160,18 @@ if ~isempty(fieldnames(motor))
 end
 %% RIG
 if ~isempty(fieldnames(RIG))
+    RIG.t_offset_main = 14094.1-0.281972;
     RIG.actuatorData(:,1) = RIG.actuatorData(:,1)/1e6;
     RIG.LoadCellData(:,1) = RIG.LoadCellData(:,1)/1e6;
     RIG.can_pressure(:,1) = RIG.can_pressure(:,1)/1e6;
     RIG.TARSData(:,1)     = RIG.TARSData(:,1)/1e6;
     RIG.EVENTData(:,1)    = RIG.EVENTData(:,1)/1e6;
+    % readjust timestamps to main valve opening
+    RIG.actuatorData(:,1) = RIG.actuatorData(:,1)-RIG.t_offset_main;
+    RIG.LoadCellData(:,1) = RIG.LoadCellData(:,1)-RIG.t_offset_main;
+    RIG.can_pressure(:,1) = RIG.can_pressure(:,1)-RIG.t_offset_main;
+    RIG.TARSData(:,1)     = RIG.TARSData(:,1)-RIG.t_offset_main;
+    RIG.EVENTData(:,1)    = RIG.EVENTData(:,1)-RIG.t_offset_main;
     % extract actuators from ADC
     RIG.mainValve         = RIG.actuatorData(RIG.actuatorData(:,2) == 5,[1,3]);
     RIG.ventValve         = RIG.actuatorData(RIG.actuatorData(:,2) == 6,[1,3]);
diff --git a/general-utilities b/general-utilities
new file mode 160000
index 0000000000000000000000000000000000000000..e9d4e4d276c44e2d5cc39608f46f37e84b58eee2
--- /dev/null
+++ b/general-utilities
@@ -0,0 +1 @@
+Subproject commit e9d4e4d276c44e2d5cc39608f46f37e84b58eee2