diff --git a/.gitattributes b/.gitattributes
index a36b49d94ecd239c2cfc9fbfab3914b2806d9826..4d5d84092bfe590c4a69bb1ea1b4515f78bbea97 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -1 +1 @@
-.mat filter=lfs diff=lfs merge=lfs -text
+*.mat filter=lfs diff=lfs merge=lfs -text
diff --git a/README.md b/functions/README.md
similarity index 100%
rename from README.md
rename to functions/README.md
diff --git a/angleTransformation/angleToQuat.m b/functions/angleTransformation/angleToQuat.m
similarity index 100%
rename from angleTransformation/angleToQuat.m
rename to functions/angleTransformation/angleToQuat.m
diff --git a/angleTransformation/dcmToAngle.m b/functions/angleTransformation/dcmToAngle.m
similarity index 100%
rename from angleTransformation/dcmToAngle.m
rename to functions/angleTransformation/dcmToAngle.m
diff --git a/angleTransformation/dcmToQuat.m b/functions/angleTransformation/dcmToQuat.m
similarity index 100%
rename from angleTransformation/dcmToQuat.m
rename to functions/angleTransformation/dcmToQuat.m
diff --git a/angleTransformation/quatToDcm.m b/functions/angleTransformation/quatToDcm.m
similarity index 100%
rename from angleTransformation/quatToDcm.m
rename to functions/angleTransformation/quatToDcm.m
diff --git a/Datcom/createFor006.m b/functions/datcom/createFor006.m
similarity index 100%
rename from Datcom/createFor006.m
rename to functions/datcom/createFor006.m
diff --git a/Datcom/datcomParser.m b/functions/datcom/datcomParser.m
similarity index 100%
rename from Datcom/datcomParser.m
rename to functions/datcom/datcomParser.m
diff --git a/eventFunctions/eventAB.m b/functions/eventFunctions/eventAB.m
similarity index 100%
rename from eventFunctions/eventAB.m
rename to functions/eventFunctions/eventAB.m
diff --git a/eventFunctions/eventApogee.m b/functions/eventFunctions/eventApogee.m
similarity index 100%
rename from eventFunctions/eventApogee.m
rename to functions/eventFunctions/eventApogee.m
diff --git a/eventFunctions/eventLanding.m b/functions/eventFunctions/eventLanding.m
similarity index 100%
rename from eventFunctions/eventLanding.m
rename to functions/eventFunctions/eventLanding.m
diff --git a/eventFunctions/eventPad.m b/functions/eventFunctions/eventPad.m
similarity index 100%
rename from eventFunctions/eventPad.m
rename to functions/eventFunctions/eventPad.m
diff --git a/eventFunctions/eventParaCut.m b/functions/eventFunctions/eventParaCut.m
similarity index 100%
rename from eventFunctions/eventParaCut.m
rename to functions/eventFunctions/eventParaCut.m
diff --git a/eventFunctions/parfor_progress.txt b/functions/eventFunctions/parfor_progress.txt
similarity index 100%
rename from eventFunctions/parfor_progress.txt
rename to functions/eventFunctions/parfor_progress.txt
diff --git a/miscellaneous/PoolWaitbar.m b/functions/miscellaneous/PoolWaitbar.m
similarity index 100%
rename from miscellaneous/PoolWaitbar.m
rename to functions/miscellaneous/PoolWaitbar.m
diff --git a/miscellaneous/addTransient.m b/functions/miscellaneous/addTransient.m
similarity index 100%
rename from miscellaneous/addTransient.m
rename to functions/miscellaneous/addTransient.m
diff --git a/miscellaneous/addTransientNew.m b/functions/miscellaneous/addTransientNew.m
similarity index 100%
rename from miscellaneous/addTransientNew.m
rename to functions/miscellaneous/addTransientNew.m
diff --git a/miscellaneous/atmosphereData.m b/functions/miscellaneous/atmosphereData.m
similarity index 100%
rename from miscellaneous/atmosphereData.m
rename to functions/miscellaneous/atmosphereData.m
diff --git a/miscellaneous/chainStruct.m b/functions/miscellaneous/chainStruct.m
similarity index 100%
rename from miscellaneous/chainStruct.m
rename to functions/miscellaneous/chainStruct.m
diff --git a/miscellaneous/checkGeometry.m b/functions/miscellaneous/checkGeometry.m
similarity index 100%
rename from miscellaneous/checkGeometry.m
rename to functions/miscellaneous/checkGeometry.m
diff --git a/miscellaneous/computeTangentBoatPoints.m b/functions/miscellaneous/computeTangentBoatPoints.m
similarity index 100%
rename from miscellaneous/computeTangentBoatPoints.m
rename to functions/miscellaneous/computeTangentBoatPoints.m
diff --git a/miscellaneous/createDissileInput.m b/functions/miscellaneous/createDissileInput.m
similarity index 100%
rename from miscellaneous/createDissileInput.m
rename to functions/miscellaneous/createDissileInput.m
diff --git a/miscellaneous/delayControl.m b/functions/miscellaneous/delayControl.m
similarity index 100%
rename from miscellaneous/delayControl.m
rename to functions/miscellaneous/delayControl.m
diff --git a/miscellaneous/exitVelocity.m b/functions/miscellaneous/exitVelocity.m
similarity index 100%
rename from miscellaneous/exitVelocity.m
rename to functions/miscellaneous/exitVelocity.m
diff --git a/miscellaneous/gaplotshape.m b/functions/miscellaneous/gaplotshape.m
similarity index 100%
rename from miscellaneous/gaplotshape.m
rename to functions/miscellaneous/gaplotshape.m
diff --git a/miscellaneous/getAlphaBeta.m b/functions/miscellaneous/getAlphaBeta.m
similarity index 100%
rename from miscellaneous/getAlphaBeta.m
rename to functions/miscellaneous/getAlphaBeta.m
diff --git a/miscellaneous/getAlphaPhi.m b/functions/miscellaneous/getAlphaPhi.m
similarity index 100%
rename from miscellaneous/getAlphaPhi.m
rename to functions/miscellaneous/getAlphaPhi.m
diff --git a/miscellaneous/interpCoeffs.m b/functions/miscellaneous/interpCoeffs.m
similarity index 100%
rename from miscellaneous/interpCoeffs.m
rename to functions/miscellaneous/interpCoeffs.m
diff --git a/miscellaneous/interpCoeffs_old.m b/functions/miscellaneous/interpCoeffs_old.m
similarity index 100%
rename from miscellaneous/interpCoeffs_old.m
rename to functions/miscellaneous/interpCoeffs_old.m
diff --git a/miscellaneous/interpLinear.m b/functions/miscellaneous/interpLinear.m
similarity index 100%
rename from miscellaneous/interpLinear.m
rename to functions/miscellaneous/interpLinear.m
diff --git a/miscellaneous/interpN.m b/functions/miscellaneous/interpN.m
similarity index 100%
rename from miscellaneous/interpN.m
rename to functions/miscellaneous/interpN.m
diff --git a/miscellaneous/nearestValSorted.m b/functions/miscellaneous/nearestValSorted.m
similarity index 100%
rename from miscellaneous/nearestValSorted.m
rename to functions/miscellaneous/nearestValSorted.m
diff --git a/miscellaneous/parfor_progress.m b/functions/miscellaneous/parfor_progress.m
similarity index 100%
rename from miscellaneous/parfor_progress.m
rename to functions/miscellaneous/parfor_progress.m
diff --git a/miscellaneous/recallOdeFcn.m b/functions/miscellaneous/recallOdeFcn.m
similarity index 100%
rename from miscellaneous/recallOdeFcn.m
rename to functions/miscellaneous/recallOdeFcn.m
diff --git a/miscellaneous/reshapeVector.m b/functions/miscellaneous/reshapeVector.m
similarity index 100%
rename from miscellaneous/reshapeVector.m
rename to functions/miscellaneous/reshapeVector.m
diff --git a/miscellaneous/safeEllipse.m b/functions/miscellaneous/safeEllipse.m
similarity index 100%
rename from miscellaneous/safeEllipse.m
rename to functions/miscellaneous/safeEllipse.m
diff --git a/miscellaneous/settingsEngineCut.m b/functions/miscellaneous/settingsEngineCut.m
similarity index 100%
rename from miscellaneous/settingsEngineCut.m
rename to functions/miscellaneous/settingsEngineCut.m
diff --git a/miscellaneous/singleGArun.m b/functions/miscellaneous/singleGArun.m
similarity index 100%
rename from miscellaneous/singleGArun.m
rename to functions/miscellaneous/singleGArun.m
diff --git a/odeFunctions/ascent.m b/functions/odeFunctions/ascent.m
similarity index 100%
rename from odeFunctions/ascent.m
rename to functions/odeFunctions/ascent.m
diff --git a/odeFunctions/ascentMultipleAB.m b/functions/odeFunctions/ascentMultipleAB.m
similarity index 100%
rename from odeFunctions/ascentMultipleAB.m
rename to functions/odeFunctions/ascentMultipleAB.m
diff --git a/odeFunctions/descentBal.m b/functions/odeFunctions/descentBal.m
similarity index 100%
rename from odeFunctions/descentBal.m
rename to functions/odeFunctions/descentBal.m
diff --git a/odeFunctions/descentParachute.m b/functions/odeFunctions/descentParachute.m
similarity index 100%
rename from odeFunctions/descentParachute.m
rename to functions/odeFunctions/descentParachute.m
diff --git a/odeFunctions/descentParafoil.m b/functions/odeFunctions/descentParafoil.m
similarity index 100%
rename from odeFunctions/descentParafoil.m
rename to functions/odeFunctions/descentParafoil.m
diff --git a/odeFunctions/launchPadFreeDyn.m b/functions/odeFunctions/launchPadFreeDyn.m
similarity index 100%
rename from odeFunctions/launchPadFreeDyn.m
rename to functions/odeFunctions/launchPadFreeDyn.m
diff --git a/roccarasoTerrainData/funZGen.m b/functions/roccarasoTerrainData/funZGen.m
similarity index 100%
rename from roccarasoTerrainData/funZGen.m
rename to functions/roccarasoTerrainData/funZGen.m
diff --git a/functions/roccarasoTerrainData/zdata.mat b/functions/roccarasoTerrainData/zdata.mat
new file mode 100644
index 0000000000000000000000000000000000000000..2c148ef7fe9ed43f473f54fc7ef2c8e5cdad7e46
--- /dev/null
+++ b/functions/roccarasoTerrainData/zdata.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c95b9364010ce84d32cb202e266e1a83bef03c78ba9fa3a315f5ce3d270324d3
+size 3441807
diff --git a/simulations/quickApogeeOnly.m b/functions/simulations/quickApogeeOnly.m
similarity index 100%
rename from simulations/quickApogeeOnly.m
rename to functions/simulations/quickApogeeOnly.m
diff --git a/simulations/stochAscentLoop.m b/functions/simulations/stochAscentLoop.m
similarity index 100%
rename from simulations/stochAscentLoop.m
rename to functions/simulations/stochAscentLoop.m
diff --git a/simulations/stochDescentBalLoop.m b/functions/simulations/stochDescentBalLoop.m
similarity index 100%
rename from simulations/stochDescentBalLoop.m
rename to functions/simulations/stochDescentBalLoop.m
diff --git a/simulations/stochDescentLoop.m b/functions/simulations/stochDescentLoop.m
similarity index 100%
rename from simulations/stochDescentLoop.m
rename to functions/simulations/stochDescentLoop.m
diff --git a/simulations/stochPostProcess.m b/functions/simulations/stochPostProcess.m
similarity index 100%
rename from simulations/stochPostProcess.m
rename to functions/simulations/stochPostProcess.m
diff --git a/simulations/stochStabilityLoop.m b/functions/simulations/stochStabilityLoop.m
similarity index 100%
rename from simulations/stochStabilityLoop.m
rename to functions/simulations/stochStabilityLoop.m
diff --git a/simulations/stochWindData.m b/functions/simulations/stochWindData.m
similarity index 100%
rename from simulations/stochWindData.m
rename to functions/simulations/stochWindData.m
diff --git a/windModels/windConstGenerator.m b/functions/windModels/windConstGenerator.m
similarity index 100%
rename from windModels/windConstGenerator.m
rename to functions/windModels/windConstGenerator.m
diff --git a/windModels/windInputGenerator.m b/functions/windModels/windInputGenerator.m
similarity index 100%
rename from windModels/windInputGenerator.m
rename to functions/windModels/windInputGenerator.m
diff --git a/windModels/windMatlabGenerator.m b/functions/windModels/windMatlabGenerator.m
similarity index 100%
rename from windModels/windMatlabGenerator.m
rename to functions/windModels/windMatlabGenerator.m
diff --git a/windModels/windVariableGenerator.m b/functions/windModels/windVariableGenerator.m
similarity index 100%
rename from windModels/windVariableGenerator.m
rename to functions/windModels/windVariableGenerator.m
diff --git a/windModels/windVertGenerator.m b/functions/windModels/windVertGenerator.m
similarity index 100%
rename from windModels/windVertGenerator.m
rename to functions/windModels/windVertGenerator.m
diff --git a/missions/2014_R1X_Roccaraso_July/empty.mat b/missions/2014_R1X_Roccaraso_July/empty.mat
new file mode 100644
index 0000000000000000000000000000000000000000..9a6a40edbe3d6f0d93e2fa57e55cd500f9c87637
--- /dev/null
+++ b/missions/2014_R1X_Roccaraso_July/empty.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6d69c0b06cdce76b064cf9de44a77fc8c5d1c43bc7306fa28887448dca5948ef
+size 1897653
diff --git a/missions/2014_R1X_Roccaraso_July/full.mat b/missions/2014_R1X_Roccaraso_July/full.mat
new file mode 100644
index 0000000000000000000000000000000000000000..67aab28c8fb3c837cb19580e46781592a9a1a4ac
--- /dev/null
+++ b/missions/2014_R1X_Roccaraso_July/full.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:65eb151c817782ef2213ed3467cab155652457f0aca5c6d3ff334fe8ee453835
+size 1899555
diff --git a/missions/2014_R1X_Roccaraso_July/simulationsData.m b/missions/2014_R1X_Roccaraso_July/simulationsData.m
new file mode 100644
index 0000000000000000000000000000000000000000..0450e125c9068adff0f1404d2124de8d7e8d2f1e
--- /dev/null
+++ b/missions/2014_R1X_Roccaraso_July/simulationsData.m
@@ -0,0 +1,256 @@
+%{
+
+ROCKET DATA FILE
+
+Project name:                          R1X
+Year of launch:                       2014
+Location of Launch:        Roccaraso (ITA)
+Date of launch:                 26/05/2014
+
+Last update:                    20/05/2021
+
+Copyright © 2021, Skyward Experimental Rocketry, AFD department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+% launchpad - roccaraso
+settings.z0 = 1414;                                     % [m] Launchpad Altitude
+settings.lpin = 0.5;                                    % [m] Distance from base of second pin
+settings.lrampa = 5.9 - settings.lpin;                  % [m] LaunchPad route (total available route)
+settings.lat0 = 41.8086605;                             % Launchpad latitude
+settings.lon0 = 14.0543387;                             % Launchpad longitude
+settings.g0 = gravitywgs84(settings.z0, settings.lat0); % Gravity costant at launch latitude and altitude
+
+settings.Local = [];                                    % vector [altLocal, TLocal, PLocal, rhoLocal], up to 4 components, the first one is mandatory
+
+%% ENGINE DETAILS
+engineName = 'R1X'; 
+tSample = linspace(0, 6, 20);
+Fthrust = @(t) -2.52*t.^4 + 33.64*t.^3 - 185.33*t.^2 + 438.92*t;
+settings.motor.expTime = tSample;
+settings.motor.expThrust = Fthrust(tSample);
+settings.motor.expM = linspace(1.980, 0, 20);
+settings.mp = 1.980;                                    % [kg]   Propellant Mass                                                
+settings.tb = tSample(end);                             % [s]    Burning time
+mc = 1.1;
+settings.mNoMot = 12.34;
+settings.ms = settings.mNoMot + mc;                     % [kg]   Structural Mass
+settings.m0 = settings.ms + settings.mp;                % [kg]   Total Mass
+settings.mnc = 0;                                       % [kg]   Nosecone Mass
+
+%% GEOMETRY/DATCOM DETAILS
+% This parameters should be the same parameters set up in MISSILE DATCOM
+% simulation.
+
+%%% Boat-tail
+% The boat-tail is simulated as a trunc of cone with the major diameter
+% equal to the rocket one (settings.C)
+settings.boatL = 0;                                        % [m] Boat-tail length
+settings.boatD = 0;                                        % [m] Boat-tail base diameter
+settings.caseInBoat = 0;                                   % [m] Length of motor case inside the boat-tail
+
+%%% Rocket ogive and central body
+settings.C = 0.045;                                        % [m] Caliber (Fuselage Diameter)
+settings.S = pi*settings.C^2/4;                            % [m^2] Cross-sectional Surface
+settings.xcg = [1.1, 1];                                   % [m] CG postion [full, empty] from nosecone base
+settings.Lnose = 0.35;                                     % [m] Nosecone Length
+settings.OgType = 'KARMAN';
+settings.rocketLength = 2.233;                             % [m] Rocket Length
+settings.Lcenter = settings.rocketLength -...
+    settings.Lnose - settings.caseInBoat;                  % [m] fuselage length
+
+%%% finset
+settings.Chord1 = 0.372;                                   % [m] attached chord length
+settings.Chord2 = 0.112;                                   % [m] free chord length
+settings.Height = 0.11;                                    % [m] fin heigth     
+settings.deltaXLE = 0.26;                                  % [m] start of Chord 2 measured from start of Chord 1
+settings.Npanel = 4;                                       % [m] number of fins
+settings.Ler = 0.003;                                      % [deg] Leading edge radius
+settings.d = 0.008;                                        % [m] rocket tip-fin distance
+settings.zupRaw = 0.0015;                                  % [m] fin semi-thickness 
+settings.LmaxuRaw = 0.0015;                                % [m] Fraction of chord from leading edge to max thickness
+
+%%% protub data
+settings.xprot = settings.Lcenter + settings.Lnose - 0.85; % axial position 
+settings.nloc = 3;                                         % number of brakes
+settings.lprot = 0.005;                                    % brakes thickness
+settings.wprot = 0.088;                                    % brakes width
+
+%% AIRBRAKES CONTROL SETTINGS
+settings.tControl = NaN;                                   % [s] time after which the airbrakes can be used
+settings.MachControl = NaN;                                % [-] Maximum Mach at which airbrakes can be used
+settings.vServo = NaN;                                     % [rad/s] Servo-motor angular velocity
+% First entry, must be always 0 !
+settings.hprot = 0;                                        % [m] Block airbrakes opening coordinate
+
+%% MASS GEOMERTY DETAILS
+% x-axis: along the fuselage
+% y-axis: right wing
+% z-axis: downward
+
+% inertias for full configuration (with all the propellant embarqued) obtained with CAD's
+settings.Ixxf = 0.27;          % [kg*m^2] Inertia to x-axis
+settings.Iyyf = 84.42;         % [kg*m^2] Inertia to y-axis
+settings.Izzf = 84.42;         % [kg*m^2] Inertia to z-axis
+
+% inertias for empty configuration (all the propellant consumed) obtained with CAD’s
+settings.Ixxe = 0.21;          % [kg*m^2] Inertia to x-axis
+settings.Iyye = 65.77;         % [kg*m^2] Inertia to y-axis
+settings.Izze = 65.77;         % [kg*m^2] Inertia to z-axis
+
+%% AERODYNAMICS DETAILS
+% These coefficients are obtained using MISSILE DATCOM
+% (after parsing with the tool datcom_parser.py)
+% The files are stored in the ../data folder with the following rule:
+% rocket_name_full.mat | rocket_name_empty.mat
+% e.g. R2a_full.mat    | R2a_empty.mat
+% Relative Path of the data files (default: ../data/). Remember the trailing slash!!
+
+% Coeffs is a 4D matrix given by Datcom that contains the aerodynamics
+% coefficient computed for the input parameters (AoA,Betas,Altitudes,Machs)
+% Note: All the parameters (AoA,Betas,Altitudes,Machs) must be the same for
+% empty and full configuration
+
+% DATA_PATH = '../data/';
+
+filename_coeffs = strcat('SRM_', engineName, '.mat');
+
+if not(isfield(settings,'autoMatProt_flag')) || settings.autoMatProt_flag == false
+    % Coefficients
+    Coeffs = load(filename_coeffs);
+    settings.Coeffs = Coeffs.CoeffsTot;
+
+    % Geometry
+    settings.Geometry = Coeffs.Geometry;
+    settings.Geometry.xcg = Coeffs.Geometry.xcg';
+
+    % State
+    settings.State.Alphas = Coeffs.State.Alphas';
+    settings.State.Betas = Coeffs.State.Betas';
+    settings.State.Altitudes = Coeffs.State.Altitudes';
+    settings.State.Machs = Coeffs.State.Machs';
+    settings.State.Controls = Coeffs.State.hprot';
+    settings.State.xcgTime = Coeffs.State.xcgTime';
+    clear('Coeffs');
+end
+%% PARACHUTES DETAILS
+settings.NdescentStages = 1; 
+settings.Npara = 3; 
+
+% parachute 1
+settings.para(1, 1).name = 'Rocket DROGUE'; 
+settings.para(1, 1).S = 1.55;      % [m^2]   Surface
+settings.para(1, 1).mass = 0.25;   % [kg]   Parachute Mass
+settings.para(1, 1).CD = 0.8;      % [/] Parachute Drag Coefficient
+settings.para(1, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 1).delay = 1;     % [s] drogue opening delay
+settings.para(1, 1).z_cut = 5000;  % [m] Final altitude of the parachute
+
+% parachute 2
+settings.para(2, 1).name = 'Rocket CHUTE 1'; 
+settings.para(2, 1).S = 17.5;      % [m^2]   Surface
+settings.para(2, 1).mass = 1.14;   % [kg]   Parachute Mass
+settings.para(2, 1).CD = 0.59;     % [/] Parachute Drag Coefficient
+settings.para(2, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(2, 1).z_cut = 2000;  % [m] Final altitude of the parachute
+
+% parachute 3
+settings.para(3, 1).name = 'Rocket CHUTE 2'; 
+settings.para(3, 1).S = 15;        % [m^2]   Surface
+settings.para(3, 1).mass = 1.466;  % [kg]   Parachute Mass
+settings.para(3, 1).CD = 0.4;      % [/] Parachute Drag Coefficient
+settings.para(3, 1).CL = 0.8;      % [/] Parachute Lift Coefficient
+settings.para(3, 1).z_cut = 0;     % [m] Final altitude of the parachute
+
+% parachute 1
+settings.para(1, 1).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 1).L = 1.5;       % [m] Shock Chord Length
+settings.para(1, 1).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 1).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 1).Vexit = 5;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 1).CX = 1.2;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(2, 1).L = 6;         % [m] Shock Chord Length
+settings.para(2, 1).K = 3000;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(2, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(2, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(2, 1).nf = 8.7;      % [/] Adimensional Opening Time
+
+% parachute 2
+settings.para(3, 1).CX = 1.2;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(3, 1).L = 6;         % [m] Shock Chord Length
+settings.para(3, 1).K = 3000;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(3, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(3, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(3, 1).nf = 8.7;      % [/] Adimensional Opening Time
+
+%% INTEGRATION OPTIONS
+settings.ode.finalTime =  2000;    % [s] Final integration time
+
+% create an option structure for the integrations:
+
+% - AbsTol is the threshold below which the value of the solution becomes unimportant
+% - RelTol is the tolerance betweeen two consecutive values
+% - Events is the event function that defines when the integration must be
+% - stopped (it has to be created)
+% - InitialStep is the highest value tried by the solver
+
+settings.ode.optionsasc1 = odeset('Events', @eventApogee, 'InitialStep', 1);       %ODE options for ascend
+settings.ode.optionsasc2 = odeset('InitialStep', 1);                                %ODE options for due to the opening delay of the parachute
+settings.ode.optionspara = odeset('Events', @eventParaCut);                       %ODE options for the parachutes
+settings.ode.optionsdesc = odeset('Events', @eventLanding);                        %ODE options for ballistic descent
+settings.ode.optionspad = odeset('Events', @eventPad);                              %ODE options for the launchpad phase
+
+% Settings for descent 6dof simulation
+settings.ode.optionsDrogue6DOF = odeset('Events', @eventParaCut, 'AbsTol', 1e-6,'RelTol', 1e-6);         %ODE options for due to cutting of the drogue chute
+settings.ode.optionsMainExt6DOF = odeset('Events', @eventMainExit, 'AbsTol', 1e-6,'RelTol', 1e-6);       %ODE options for due to the extraction of the main chute
+settings.ode.optionsMain6DOF = odeset('Events', @eventLanding, 'AbsTol', 1e-6,'RelTol', 1e-6);            %ODE options to terminate descent phase
+
+
+%% STOCHASTIC DETAILS
+%%% launch probability details
+settings.stoch.prob.xLim = 2e3;                    % Max ovest displacement [m]
+settings.stoch.prob.VLim = 50;                     % Max drogue velocity [Pa]
+settings.stoch.prob.XCPLim = 1.5;                  % Min XCP
+
+%%% Safe Ellipse - roccaraso
+settings.prob.SafeEllipse.a = 1100;
+settings.prob.SafeEllipse.b = 2800;
+settings.prob.SafeEllipse.x0  = 0;
+settings.prob.SafeEllipse.y0 = -300;
+settings.prob.SafeEllipse.alpha = 10;
+
+
+%% LANDING POINTS
+% satellite maps of the landing zone 
+% delta limit on the coordinates (for the landing map)
+settings.limLat = 0.04; 
+settings.limLon = 0.025;
+
+
+%% HRE - compatibility 
+settings.HREmot = false;
+
+%%% efflux data
+settings.motor.Ae = 0; 
+settings.motor.Pe = zeros(size(settings.motor.expTime)); 
+
+%%% intertias 
+settings.Ixx = (settings.motor.expTime./settings.tb).*(settings.Ixxe - settings.Ixxf) + settings.Ixxf; 
+settings.Iyy = (settings.motor.expTime./settings.tb).*(settings.Iyye - settings.Iyyf) + settings.Iyyf;
+settings.Izz = (settings.motor.expTime./settings.tb).*(settings.Izze - settings.Izzf) + settings.Izzf;
+settings.I = [settings.Ixx; settings.Iyy; settings.Izz]; 
+settings.Idot = ones(size(settings.I)).*(settings.I(:, end) - settings.I(:, 1))./settings.tb; 
+
+%%% mass & center of gravity
+settings.mTotalTime = settings.ms + settings.motor.expM; 
+settings.xcg = (settings.motor.expTime./settings.tb).*(settings.xcg(2) - settings.xcg(1)) + settings.xcg(1);
+
+%%% times 
+settings.timeEngineCut = inf; 
\ No newline at end of file
diff --git a/missions/2017_R2A_Sardinia_July/empty.mat b/missions/2017_R2A_Sardinia_July/empty.mat
new file mode 100644
index 0000000000000000000000000000000000000000..4546f4b9de7962c1df54769516130d83312bd845
--- /dev/null
+++ b/missions/2017_R2A_Sardinia_July/empty.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7e856bc225a9dd07646a095754a97fa7420d6fd0682621efc8e877f774ad3462
+size 1593391
diff --git a/missions/2017_R2A_Sardinia_July/full.mat b/missions/2017_R2A_Sardinia_July/full.mat
new file mode 100644
index 0000000000000000000000000000000000000000..045d118c79d2b6dd7a1dea22de26974f1e17cb50
--- /dev/null
+++ b/missions/2017_R2A_Sardinia_July/full.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4a455d04ccdfddef2bf2e86ef8311257fae68aa6afc172cd31678935c638645a
+size 1589290
diff --git a/missions/2017_R2A_Sardinia_July/simulationsData.m b/missions/2017_R2A_Sardinia_July/simulationsData.m
new file mode 100644
index 0000000000000000000000000000000000000000..6ba53a68d46993e15823fbdfe3b2213d203a6fe5
--- /dev/null
+++ b/missions/2017_R2A_Sardinia_July/simulationsData.m
@@ -0,0 +1,266 @@
+%{
+
+ROCKET DATA FILE
+
+Project name:                          R2A
+Year of launch:                       2017
+Location of Launch:         Sardinia (ITA)
+Date of launch:                        TBD
+
+Last update:                    20/11/2021
+
+Copyright © 2021, Skyward Experimental Rocketry, AFD department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+% launchpad - Sardinia
+settings.z0 = 5;                                        % [m] Launchpad Altitude
+settings.lpin = 0.5;                                    % [m] Distance from base of second pin
+settings.lrampa = 5.9 - settings.lpin;                  % [m] LaunchPad route (total available route)
+settings.lat0 = 39.552709;                              % Launchpad latitude
+settings.lon0 = 9.652400;                               % Launchpad longitude
+settings.g0 = gravitywgs84(settings.z0, settings.lat0); % Gravity costant at launch latitude and altitude
+
+settings.Local = [];                                    % vector [altLocal, TLocal, PLocal, rhoLocal], up to 4 components, the first one is mandatory
+
+%% ENGINE DETAILS
+% load motors data 
+filename = '../Motors.mat';
+load(filename)
+
+engineName = 'O8000-WT';
+crossCase = false;   % true if you want a Cesaroni case with an Aerotech reload, and viceversa
+
+iMotor = [Motors.MotorName] == engineName;
+settings.motor.expTime = Motors(iMotor).t;
+settings.motor.expThrust = Motors(iMotor).T;
+settings.motor.expM = Motors(iMotor).m;
+settings.mp = Motors(iMotor).mp;                        % [kg]   Propellant Mass                                                
+settings.tb = Motors(iMotor).t(end) ;                   % [s]    Burning time
+if not(crossCase)
+    mc = Motors(iMotor).mc;                             % [kg]   Case mass
+else
+    mc = Motors(iMotor).mcCross;
+end
+settings.mNoMot = 35.0890;
+settings.ms = settings.mNoMot + mc;                     % [kg]   Structural Mass
+settings.m0 = settings.ms + settings.mp;                % [kg]   Total Mass
+settings.mnc = 0;                                       % [kg]   Nosecone Mass
+
+%% GEOMETRY/DATCOM DETAILS
+% This parameters should be the same parameters set up in MISSILE DATCOM
+% simulation.
+
+%%% Boat-tail
+% The boat-tail is simulated as a trunc of cone with the major diameter
+% equal to the rocket one (settings.C)
+settings.boatL = 0;                                        % [m] Boat-tail length
+settings.boatD = 0;                                        % [m] Boat-tail base diameter
+settings.caseInBoat = 0;                                   % [m] Length of motor case inside the boat-tail
+
+%%% Rocket ogive and central body
+settings.C = 0.174;                                        % [m] Caliber (Fuselage Diameter)
+settings.S = pi*settings.C^2/4;                            % [m^2] Cross-sectional Surface
+settings.xcg = [2.194, 1.868];                             % [m] CG postion [full, empty] from nosecone base
+settings.Lnose = 0.9;                                      % [m] Nosecone Length
+settings.OgType = 'KARMAN';
+settings.rocketLength = 4.43619;                           % [m] Rocket Length
+settings.Lcenter = settings.rocketLength -...
+    settings.Lnose - settings.caseInBoat;                  % [m] fuselage length
+
+%%% finset
+settings.Chord1 = 0.17;                                    % [m] attached chord length
+settings.Chord2 = 0.08;                                    % [m] free chord length
+settings.Height = 0.08;                                    % [m] fin heigth     
+settings.deltaXLE = 0.09;                                  % [m] start of Chord 2 measured from start of Chord 1
+settings.Npanel = 4;                                       % [m] number of fins
+settings.Ler = 0.003;                                      % [deg] Leading edge radius
+settings.d = 0.008;                                        % [m] rocket tip-fin distance
+settings.zupRaw = 0.0015;                                  % [m] fin semi-thickness 
+settings.LmaxuRaw = 0.0015;                                % [m] Fraction of chord from leading edge to max thickness
+
+%%% protub data
+settings.xprot = settings.Lcenter + settings.Lnose - 0.85; % axial position 
+settings.nloc = 3;                                         % number of brakes
+settings.lprot = 0.005;                                    % brakes thickness
+settings.wprot = 0.088;                                    % brakes width
+
+%% AIRBRAKES CONTROL SETTINGS
+settings.tControl = NaN;                                   % [s] time after which the airbrakes can be used
+settings.MachControl = NaN;                                % [-] Maximum Mach at which airbrakes can be used
+settings.vServo = NaN;                                     % [rad/s] Servo-motor angular velocity
+% First entry, must be always 0 !
+settings.hprot = 0;                                        % [m] Block airbrakes opening coordinate
+
+%% MASS GEOMERTY DETAILS
+% x-axis: along the fuselage
+% y-axis: right wing
+% z-axis: downward
+
+% inertias for full configuration (with all the propellant embarqued) obtained with CAD's
+settings.Ixxf = 0.27;          % [kg*m^2] Inertia to x-axis
+settings.Iyyf = 84.42;         % [kg*m^2] Inertia to y-axis
+settings.Izzf = 84.42;         % [kg*m^2] Inertia to z-axis
+
+% inertias for empty configuration (all the propellant consumed) obtained with CAD’s
+settings.Ixxe = 0.21;          % [kg*m^2] Inertia to x-axis
+settings.Iyye = 65.77;         % [kg*m^2] Inertia to y-axis
+settings.Izze = 65.77;         % [kg*m^2] Inertia to z-axis
+
+%% AERODYNAMICS DETAILS
+% These coefficients are obtained using MISSILE DATCOM
+% (after parsing with the tool datcom_parser.py)
+% The files are stored in the ../data folder with the following rule:
+% rocket_name_full.mat | rocket_name_empty.mat
+% e.g. R2a_full.mat    | R2a_empty.mat
+% Relative Path of the data files (default: ../data/). Remember the trailing slash!!
+
+% Coeffs is a 4D matrix given by Datcom that contains the aerodynamics
+% coefficient computed for the input parameters (AoA,Betas,Altitudes,Machs)
+% Note: All the parameters (AoA,Betas,Altitudes,Machs) must be the same for
+% empty and full configuration
+
+% DATA_PATH = '../data/';
+
+filename_coeffs = strcat('SRM_', engineName, '.mat');
+
+if not(isfield(settings,'autoMatProt_flag')) || settings.autoMatProt_flag == false
+    % Coefficients
+    Coeffs = load(filename_coeffs);
+    settings.Coeffs = Coeffs.CoeffsTot;
+
+    % Geometry
+    settings.Geometry = Coeffs.Geometry;
+    settings.Geometry.xcg = Coeffs.Geometry.xcg';
+
+    % State
+    settings.State.Alphas = Coeffs.State.Alphas';
+    settings.State.Betas = Coeffs.State.Betas';
+    settings.State.Altitudes = Coeffs.State.Altitudes';
+    settings.State.Machs = Coeffs.State.Machs';
+    settings.State.Controls = Coeffs.State.hprot';
+    settings.State.xcgTime = Coeffs.State.xcgTime';
+    clear('Coeffs');
+end
+%% PARACHUTES DETAILS
+settings.NdescentStages = 1; 
+settings.Npara = 3; 
+
+% parachute 1
+settings.para(1, 1).name = 'Rocket DROGUE 1'; 
+settings.para(1, 1).S = 1.55;      % [m^2]   Surface
+settings.para(1, 1).mass = 0.25;   % [kg]   Parachute Mass
+settings.para(1, 1).CD = 0.8;      % [/] Parachute Drag Coefficient
+settings.para(1, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 1).delay = 1;     % [s] drogue opening delay
+settings.para(1, 1).z_cut = 5000;  % [m] Final altitude of the parachute
+
+% parachute 2
+settings.para(2, 1).name = 'Rocket CHUTE 1'; 
+settings.para(2, 1).S = 17.5;      % [m^2]   Surface
+settings.para(2, 1).mass = 1.14;   % [kg]   Parachute Mass
+settings.para(2, 1).CD = 0.59;     % [/] Parachute Drag Coefficient
+settings.para(2, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(2, 1).z_cut = 2000;  % [m] Final altitude of the parachute
+
+% parachute 3
+settings.para(3, 1).name = 'Rocket CHUTE 2'; 
+settings.para(3, 1).S = 15;        % [m^2]   Surface
+settings.para(3, 1).mass = 1.466;  % [kg]   Parachute Mass
+settings.para(3, 1).CD = 0.4;      % [/] Parachute Drag Coefficient
+settings.para(3, 1).CL = 0.8;      % [/] Parachute Lift Coefficient
+settings.para(3, 1).z_cut = 0;     % [m] Final altitude of the parachute
+
+% parachute 1
+settings.para(1, 1).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 1).L = 1.5;       % [m] Shock Chord Length
+settings.para(1, 1).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 1).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 1).Vexit = 5;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 1).CX = 1.2;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(2, 1).L = 6;         % [m] Shock Chord Length
+settings.para(2, 1).K = 3000;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(2, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(2, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(2, 1).nf = 8.7;      % [/] Adimensional Opening Time
+
+% parachute 2
+settings.para(3, 1).CX = 1.2;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(3, 1).L = 6;         % [m] Shock Chord Length
+settings.para(3, 1).K = 3000;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(3, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(3, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(3, 1).nf = 8.7;      % [/] Adimensional Opening Time
+
+%% INTEGRATION OPTIONS
+settings.ode.finalTime =  2000;    % [s] Final integration time
+
+% create an option structure for the integrations:
+
+% - AbsTol is the threshold below which the value of the solution becomes unimportant
+% - RelTol is the tolerance betweeen two consecutive values
+% - Events is the event function that defines when the integration must be
+% - stopped (it has to be created)
+% - InitialStep is the highest value tried by the solver
+
+settings.ode.optionsasc1 = odeset('Events', @eventApogee, 'InitialStep', 1);       %ODE options for ascend
+settings.ode.optionsasc2 = odeset('InitialStep', 1);                                %ODE options for due to the opening delay of the parachute
+settings.ode.optionspara = odeset('Events', @eventParaCut);                       %ODE options for the parachutes
+settings.ode.optionsdesc = odeset('Events', @eventLanding);                        %ODE options for ballistic descent
+settings.ode.optionspad = odeset('Events', @eventPad);                              %ODE options for the launchpad phase
+
+% Settings for descent 6dof simulation
+settings.ode.optionsDrogue6DOF = odeset('Events', @eventParaCut, 'AbsTol', 1e-6,'RelTol', 1e-6);         %ODE options for due to cutting of the drogue chute
+settings.ode.optionsMainExt6DOF = odeset('Events', @eventMainExit, 'AbsTol', 1e-6,'RelTol', 1e-6);       %ODE options for due to the extraction of the main chute
+settings.ode.optionsMain6DOF = odeset('Events', @eventLanding, 'AbsTol', 1e-6,'RelTol', 1e-6);            %ODE options to terminate descent phase
+
+
+%% STOCHASTIC DETAILS
+%%% launch probability details
+settings.stoch.prob.xLim = 2e3;                    % Max ovest displacement [m]
+settings.stoch.prob.VLim = 50;                     % Max drogue velocity [Pa]
+settings.stoch.prob.XCPLim = 1.5;                  % Min XCP
+
+%%% Safe Ellipse - roccaraso
+settings.prob.SafeEllipse.a = 1100;
+settings.prob.SafeEllipse.b = 2800;
+settings.prob.SafeEllipse.x0  = 0;
+settings.prob.SafeEllipse.y0 = -300;
+settings.prob.SafeEllipse.alpha = 10;
+
+
+%% LANDING POINTS
+% satellite maps of the landing zone 
+% delta limit on the coordinates (for the landing map)
+settings.limLat = 0.04; 
+settings.limLon = 0.025;
+
+
+%% HRE - compatibility 
+settings.HREmot = false;
+
+%%% efflux data
+settings.motor.Ae = 0; 
+settings.motor.Pe = zeros(size(settings.motor.expTime)); 
+
+%%% intertias 
+settings.Ixx = (settings.motor.expTime./settings.tb).*(settings.Ixxe - settings.Ixxf) + settings.Ixxf; 
+settings.Iyy = (settings.motor.expTime./settings.tb).*(settings.Iyye - settings.Iyyf) + settings.Iyyf;
+settings.Izz = (settings.motor.expTime./settings.tb).*(settings.Izze - settings.Izzf) + settings.Izzf;
+settings.I = [settings.Ixx; settings.Iyy; settings.Izz]; 
+settings.Idot = ones(size(settings.I)).*(settings.I(:, end) - settings.I(:, 1))./settings.tb; 
+
+%%% mass & center of gravity
+settings.mTotalTime = settings.ms + settings.motor.expM; 
+settings.xcg = (settings.motor.expTime./settings.tb).*(settings.xcg(2) - settings.xcg(1)) + settings.xcg(1);
+
+%%% times 
+settings.timeEngineCut = inf; 
+
diff --git a/missions/2019_HermesV1_Roccaraso_November/empty.mat b/missions/2019_HermesV1_Roccaraso_November/empty.mat
new file mode 100644
index 0000000000000000000000000000000000000000..e48b6eb4ea860adcedb47c37eb1a83525562a0e1
--- /dev/null
+++ b/missions/2019_HermesV1_Roccaraso_November/empty.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b46d28ab8be627d1e91ed0d10f30c339b997cb4e3aefcce5675aa4b53d53d26b
+size 1997398
diff --git a/missions/2019_HermesV1_Roccaraso_November/full.mat b/missions/2019_HermesV1_Roccaraso_November/full.mat
new file mode 100644
index 0000000000000000000000000000000000000000..e48b6eb4ea860adcedb47c37eb1a83525562a0e1
--- /dev/null
+++ b/missions/2019_HermesV1_Roccaraso_November/full.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b46d28ab8be627d1e91ed0d10f30c339b997cb4e3aefcce5675aa4b53d53d26b
+size 1997398
diff --git a/missions/2019_HermesV1_Roccaraso_November/simulationsData.m b/missions/2019_HermesV1_Roccaraso_November/simulationsData.m
new file mode 100644
index 0000000000000000000000000000000000000000..4783e9792544402426a2de8bfd038990b3a369cb
--- /dev/null
+++ b/missions/2019_HermesV1_Roccaraso_November/simulationsData.m
@@ -0,0 +1,251 @@
+%{
+
+ROCKET DATA FILE
+
+Project name:                    Hermes V1
+Year of launch:                       2019
+Location of Launch:        Roccaraso (ITA)
+Date of launch:                 15/11/2019
+
+Last update:                    16/11/2019
+
+Copyright © 2021, Skyward Experimental Rocketry, AFD department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+% launchpad - roccaraso
+settings.z0 = 1416;                                     % [m] Launchpad Altitude
+settings.lpin = 1.150;                                  % [m] Distance from base of second pin
+settings.lrampa = 5.9 - settings.lpin;                  % [m] LaunchPad route (total available route)
+settings.lat0 = 41.810093;                              % Launchpad latitude
+settings.lon0 = 14.052546;                              % Launchpad longitude
+
+settings.g0 = gravitywgs84(settings.z0, settings.lat0); % Gravity costant at launch latitude and altitude
+
+%% ENGINE DETAILS
+% load motors data 
+filename = '../Motors.mat';
+load(filename)
+
+% name = 'K828FJ';
+engineName = 'K550W';
+crossCase = false;   % true if you want a Cesaroni case with an Aerotech reload, and viceversa
+
+iMotor = [Motors.MotorName] == engineName;
+settings.motor.expTime = Motors(iMotor).t;
+settings.motor.expThrust = Motors(iMotor).T;
+settings.motor.expM = Motors(iMotor).m;
+settings.mp = Motors(iMotor).mp;                        % [kg]   Propellant Mass                                                
+settings.tb = Motors(iMotor).t(end) ;                   % [s]    Burning time
+if not(crossCase)
+    mc = Motors(iMotor).mc;                             % [kg]   Case mass
+else
+    mc = Motors(iMotor).mcCross;
+end
+settings.mNoMot = 7.0423;
+settings.ms = settings.mNoMot + mc;                     % [kg]   Structural Mass
+settings.m0 = settings.ms + settings.mp;                % [kg]   Total Mass
+settings.mnc = 0;                                       % [kg]   Nosecone Mass
+
+%% GEOMETRY/DATCOM DETAILS
+% This parameters should be the same parameters set up in MISSILE DATCOM
+% simulation.
+
+%%% Boat-tail
+% The boat-tail is simulated as a trunc of cone with the major diameter
+% equal to the rocket one (settings.C)
+settings.boatL = 0;                                        % [m] Boat-tail length
+settings.boatD = 0;                                        % [m] Boat-tail base diameter
+settings.caseInBoat = 0;                                   % [m] Length of motor case inside the boat-tail
+
+%%% Rocket ogive and central body
+settings.C = 0.09;                                         % [m] Caliber (Fuselage Diameter)
+settings.S = pi*settings.C^2/4;                            % [m^2] Cross-sectional Surface
+settings.xcg = [1.2, 1.2];                                 % [m] CG postion [full, empty] from nosecone base
+settings.Lnose = 0.30;                                     % [m] Nosecone Length
+settings.OgType = 'KARMAN';
+settings.rocketLength = 1.895;                             % [m] Rocket Length
+settings.Lcenter = settings.rocketLength -...
+    settings.Lnose - settings.caseInBoat;                  % [m] fuselage length
+
+%%% finset
+settings.Chord1 = 0.17;                                    % [m] attached chord length
+settings.Chord2 = 0.08;                                    % [m] free chord length
+settings.Height = 0.08;                                    % [m] fin heigth     
+settings.deltaXLE = 0.09;                                  % [m] start of Chord 2 measured from start of Chord 1
+settings.Npanel = 3;                                       % [m] number of fins
+settings.Ler = 0.003;                                      % [deg] Leading edge radius
+settings.d = 0.008;                                        % [m] rocket tip-fin distance
+settings.zupRaw = 0.0015;                                  % [m] fin semi-thickness 
+settings.LmaxuRaw = 0.0015;                                % [m] Fraction of chord from leading edge to max thickness
+
+%%% protub data
+settings.xprot = settings.Lcenter + settings.Lnose - 0.85; % axial position 
+settings.nloc = 3;                                         % number of brakes
+settings.lprot = 0.005;                                    % brakes thickness
+settings.wprot = 0.088;                                    % brakes width
+
+%% AIRBRAKES CONTROL SETTINGS
+settings.tControl = NaN;                                   % [s] time after which the airbrakes can be used
+settings.MachControl = NaN;                                % [-] Maximum Mach at which airbrakes can be used
+settings.vServo = NaN;                                     % [rad/s] Servo-motor angular velocity
+% First entry, must be always 0 !
+settings.hprot = 0;                                        % [m] Block airbrakes opening coordinate
+
+%% MASS GEOMERTY DETAILS
+% x-axis: along the fuselage
+% y-axis: right wing
+% z-axis: downward
+
+% inertias for full configuration (with all the propellant embarqued) obtained with CAD's
+settings.Ixxf = 0.004612382;          % [kg*m^2] Inertia to x-axis
+settings.Iyyf = 1.194050037;          % [kg*m^2] Inertia to y-axis
+settings.Izzf = 1.194116615;          % [kg*m^2] Inertia to z-axis
+
+% inertias for empty configuration (all the propellant consumed) obtained with CAD’s
+settings.Ixxe = 0.004293335;          % [kg*m^2] Inertia to x-axis
+settings.Iyye = 0.931311998;          % [kg*m^2] Inertia to y-axis
+settings.Izze = 0.931376985;          % [kg*m^2] Inertia to z-axis
+
+%% AERODYNAMICS DETAILS
+% These coefficients are obtained using MISSILE DATCOM
+% (after parsing with the tool datcom_parser.py)
+% The files are stored in the ../data folder with the following rule:
+% rocket_name_full.mat | rocket_name_empty.mat
+% e.g. R2a_full.mat    | R2a_empty.mat
+% Relative Path of the data files (default: ../data/). Remember the trailing slash!!
+
+% Coeffs is a 4D matrix given by Datcom that contains the aerodynamics
+% coefficient computed for the input parameters (AoA,Betas,Altitudes,Machs)
+% Note: All the parameters (AoA,Betas,Altitudes,Machs) must be the same for
+% empty and full configuration
+
+% DATA_PATH = '../data/';
+
+filename_coeffs = strcat('SRM_', engineName, '.mat');
+
+if not(isfield(settings,'autoMatProt_flag')) || settings.autoMatProt_flag == false
+    % Coefficients
+    Coeffs = load(filename_coeffs);
+    settings.Coeffs = Coeffs.CoeffsTot;
+
+    % Geometry
+    settings.Geometry = Coeffs.Geometry;
+    settings.Geometry.xcg = Coeffs.Geometry.xcg';
+
+    % State
+    settings.State.Alphas = Coeffs.State.Alphas';
+    settings.State.Betas = Coeffs.State.Betas';
+    settings.State.Altitudes = Coeffs.State.Altitudes';
+    settings.State.Machs = Coeffs.State.Machs';
+    settings.State.Controls = Coeffs.State.hprot';
+    settings.State.xcgTime = Coeffs.State.xcgTime';
+    clear('Coeffs');
+end
+
+
+%% PARACHUTES DETAILS
+settings.NdescentStages = 1; 
+settings.Npara = 2; 
+
+% parachute 1
+settings.para(1, 1).name = "Rocket DROGUE";
+settings.para(1, 1).S = 0.4;       % [m^2]   Surface
+settings.para(1, 1).mass = 0.2;    % [kg]   Parachute Mass
+settings.para(1, 1).CD = 0.78;     % [/] Parachute Drag Coefficient
+settings.para(1, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 1).delay = 1;     % [s] drogue opening delay
+settings.para(1, 1).z_cut = 450;   % [m] Final altitude of the parachute
+
+% parachute 2
+settings.para(2, 1).name = "Rocket MAIN chute";
+settings.para(2, 1).S = 10.5;      % [m^2]   Surface
+settings.para(2, 1).mass = 1.5;    % [kg]   Parachute Mass
+settings.para(2, 1).CD = 0.7;      % [/] Parachute Drag Coefficient
+settings.para(2, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(2, 1).z_cut = 0;     % [m] Final altitude of the parachute
+
+% parachute 1
+settings.para(1, 1).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 1).L = 1.5;         % [m] Shock Chord Length
+settings.para(1, 1).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 1).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 1).Vexit = 5;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 1).CX = 1.2;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(2, 1).L = 6;         % [m] Shock Chord Length
+settings.para(2, 1).K = 3000;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(2, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(2, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(2, 1).nf = 8.7;      % [/] Adimensional Opening Time
+
+%% INTEGRATION OPTIONS
+settings.ode.finalTime =  2000;    % [s] Final integration time
+
+% create an option structure for the integrations:
+
+% - AbsTol is the threshold below which the value of the solution becomes unimportant
+% - RelTol is the tolerance betweeen two consecutive values
+% - Events is the event function that defines when the integration must be
+% - stopped (it has to be created)
+% - InitialStep is the highest value tried by the solver
+
+settings.ode.optionsasc1 = odeset('Events', @eventApogee, 'InitialStep', 1);       %ODE options for ascend
+settings.ode.optionsasc2 = odeset('InitialStep', 1);                                %ODE options for due to the opening delay of the parachute
+settings.ode.optionspara = odeset('Events', @eventParaCut);                       %ODE options for the parachutes
+settings.ode.optionsdesc = odeset('Events', @eventLanding);                        %ODE options for ballistic descent
+settings.ode.optionspad = odeset('Events', @eventPad);                              %ODE options for the launchpad phase
+
+% Settings for descent 6dof simulation
+settings.ode.optionsDrogue6DOF = odeset('Events', @eventParaCut, 'AbsTol', 1e-6,'RelTol', 1e-6);         %ODE options for due to cutting of the drogue chute
+settings.ode.optionsMainExt6DOF = odeset('Events', @eventMainExit, 'AbsTol', 1e-6,'RelTol', 1e-6);       %ODE options for due to the extraction of the main chute
+settings.ode.optionsMain6DOF = odeset('Events', @eventLanding, 'AbsTol', 1e-6,'RelTol', 1e-6);            %ODE options to terminate descent phase
+
+
+%% STOCHASTIC DETAILS
+%%% launch probability details
+settings.stoch.prob.xLim = 2e3;                    % Max ovest displacement [m]
+settings.stoch.prob.VLim = 50;                     % Max drogue velocity [Pa]
+settings.stoch.prob.XCPLim = 1.5;                  % Min XCP
+
+%%% Safe Ellipse - roccaraso
+settings.prob.SafeEllipse.a = 1100;
+settings.prob.SafeEllipse.b = 2800;
+settings.prob.SafeEllipse.x0  = 0;
+settings.prob.SafeEllipse.y0 = -300;
+settings.prob.SafeEllipse.alpha = 10;
+
+
+%% LANDING POINTS
+% satellite maps of the landing zone 
+% delta limit on the coordinates (for the landing map)
+settings.limLat = 0.04; 
+settings.limLon = 0.025;
+
+
+%% HRE - compatibility 
+settings.HREmot = false;
+
+%%% efflux data
+settings.motor.Ae = 0; 
+settings.motor.Pe = zeros(size(settings.motor.expTime)); 
+
+%%% intertias 
+settings.Ixx = (settings.motor.expTime./settings.tb).*(settings.Ixxe - settings.Ixxf) + settings.Ixxf; 
+settings.Iyy = (settings.motor.expTime./settings.tb).*(settings.Iyye - settings.Iyyf) + settings.Iyyf;
+settings.Izz = (settings.motor.expTime./settings.tb).*(settings.Izze - settings.Izzf) + settings.Izzf;
+settings.I = [settings.Ixx; settings.Iyy; settings.Izz]; 
+settings.Idot = ones(size(settings.I)).*(settings.I(:, end) - settings.I(:, 1))./settings.tb; 
+
+%%% mass & center of gravity
+settings.mTotalTime = settings.ms + settings.motor.expM; 
+settings.xcg = (settings.motor.expTime./settings.tb).*(settings.xcg(2) - settings.xcg(1)) + settings.xcg(1);
+
+%%% times 
+settings.timeEngineCut = inf; 
\ No newline at end of file
diff --git a/missions/2021_Lynx_Portugal_October/CAinterpCoeffs.mat b/missions/2021_Lynx_Portugal_October/CAinterpCoeffs.mat
new file mode 100644
index 0000000000000000000000000000000000000000..39188e3da0fe7cb1c4a37c17fac78550d03ac52d
--- /dev/null
+++ b/missions/2021_Lynx_Portugal_October/CAinterpCoeffs.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6a2a555a45cc83ed6e92551c9b1be33d07bf6fb13af81524f603ddd101d2cb60
+size 458
diff --git a/missions/2021_Lynx_Portugal_October/SRM_M2000RPortugal.mat b/missions/2021_Lynx_Portugal_October/SRM_M2000RPortugal.mat
new file mode 100644
index 0000000000000000000000000000000000000000..5b0ce0aac542de53f7c06a4c40e699cce31ec661
--- /dev/null
+++ b/missions/2021_Lynx_Portugal_October/SRM_M2000RPortugal.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:87123942f5655d390571b1fee3a0d4c658b87070cab647501d2127533ba5455f
+size 24181749
diff --git a/missions/2021_Lynx_Portugal_October/empty.mat b/missions/2021_Lynx_Portugal_October/empty.mat
new file mode 100644
index 0000000000000000000000000000000000000000..e5435be88ab6381b37545e78184065fc24978ceb
--- /dev/null
+++ b/missions/2021_Lynx_Portugal_October/empty.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:51a9ceb243bf287093c9f4a62846b734ed90d9a42535c1908d2367b0ad4e9a32
+size 6000781
diff --git a/missions/2021_Lynx_Portugal_October/emptyOld.mat b/missions/2021_Lynx_Portugal_October/emptyOld.mat
new file mode 100644
index 0000000000000000000000000000000000000000..b3ea379266fbf18e117fd6c877eeae6aae331827
--- /dev/null
+++ b/missions/2021_Lynx_Portugal_October/emptyOld.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:70180664d8d4c1cd96b95c2f62bc32a34a1b0fad430302c727c281d4359973ba
+size 6025068
diff --git a/missions/2021_Lynx_Portugal_October/full.mat b/missions/2021_Lynx_Portugal_October/full.mat
new file mode 100644
index 0000000000000000000000000000000000000000..0be86973e4bdf99e2d435fce584db86f1bd2202b
--- /dev/null
+++ b/missions/2021_Lynx_Portugal_October/full.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d88cd8510e7aa53be4865aff9aacdefcaf64fcf5eecd72599123e9cde09a69a5
+size 1990203
diff --git a/missions/2021_Lynx_Portugal_October/fullOld.mat b/missions/2021_Lynx_Portugal_October/fullOld.mat
new file mode 100644
index 0000000000000000000000000000000000000000..023448f7bc9e0aba558e71e25495f6609b773d63
--- /dev/null
+++ b/missions/2021_Lynx_Portugal_October/fullOld.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ffebc04ffcefd81db780eebea15caafe698d250041c009db259dc38cb97b8904
+size 1997639
diff --git a/missions/2021_Lynx_Portugal_October/simulationsData.m b/missions/2021_Lynx_Portugal_October/simulationsData.m
new file mode 100644
index 0000000000000000000000000000000000000000..0dca0961c95fcc2e5a2d9a39ebc15b7908c0415f
--- /dev/null
+++ b/missions/2021_Lynx_Portugal_October/simulationsData.m
@@ -0,0 +1,254 @@
+%{
+
+ROCKET DATA FILE
+
+Project name:                         Lynx
+Year of launch:                       2021
+Location of Launch:      Pont de Sor (POR)
+Date of launch:                 13/10/2021
+
+Last update:                    25/10/2021
+
+Copyright © 2021, Skyward Experimental Rocketry, AFD department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+% launchpad - Pont de Sor
+settings.z0 = 160;                                      % [m] Launchpad Altitude
+settings.lpin = 1.150;                                  % [m] Distance from base of second pin
+settings.lrampa = 12 - settings.lpin;                   % [m] LaunchPad route (total available route)
+settings.lat0 = 39.388727;                              % Launchpad latitude
+settings.lon0 = -8.287842;                              % Launchpad longitude
+settings.g0 = gravitywgs84(settings.z0, settings.lat0); % Gravity costant at launch latitude and altitude
+
+settings.Local = [settings.z0, 299.15, 99634.8516];     % vector [altLocal, TLocal, PLocal, rhoLocal], up to 4 components
+
+%% ENGINE DETAILS
+% load motors data 
+filename = '../Motors.mat';
+load(filename);
+
+engineName = 'M2000RPortugal';
+crossCase = true;   % true if you want a Cesaroni case with an Aerotech reload, and viceversa
+
+iMotor = [Motors.MotorName] == engineName;
+settings.motor.expTime = Motors(iMotor).t;
+settings.motor.expThrust = Motors(iMotor).T;
+settings.motor.expM = Motors(iMotor).m;
+settings.mp = Motors(iMotor).mp;                        % [kg]   Propellant Mass                                                
+settings.tb = Motors(iMotor).t(end) ;                   % [s]    Burning time
+if not(crossCase)
+    mc = Motors(iMotor).mc;                             % [kg]   Case mass
+else
+    mc = Motors(iMotor).mcCross;
+end
+settings.mNoMot = 17.94;
+settings.ms = settings.mNoMot + mc;                     % [kg]   Structural Mass
+settings.m0 = settings.ms + settings.mp;                % [kg]   Total Mass
+settings.mnc = 0;                                       % [kg]   Nosecone Mass
+
+%% GEOMETRY/DATCOM DETAILS
+% This parameters should be the same parameters set up in MISSILE DATCOM
+% simulation.
+
+%%% Boat-tail
+% The boat-tail is simulated as a trunc of cone with the major diameter
+% equal to the rocket one (settings.C)
+settings.boatL = 0;                   % [m] Boat-tail length
+settings.boatD = 0;                   % [m] Boat-tail base diameter
+
+%%% Rocket ogive and central body
+settings.C = 0.15;                    % [m] Caliber (Fuselage Diameter)
+settings.S = pi*settings.C^2/4;       % [m^2] Cross-sectional Surface
+settings.xcg = [1.24, 1.084];         % [m] CG postion [full, empty] from nosecone base
+settings.Lnose = 0.26;                % [m] Nosecone Length
+settings.OgType = 'KARMAN';
+settings.NosePower = 0;               % [-] Nosecone power type parameter
+settings.rocketLength = 2.496;        % [m] Rocket Length
+settings.Lcenter = settings.rocketLength -...
+    settings.Lnose - settings.boatL;  % [m] fuselage length
+
+%%% finset
+settings.Chord1 = 0.345;              % [m] attached chord length
+settings.Chord2 = 0.12;               % [m] free chord length
+settings.Height = 0.117;              % [m] fin heigth
+settings.deltaXLE = 0.225;            % [m] start of Chord 2 measured from start of Chord 1
+settings.Npanel = 3;                  % [m] number of fins
+settings.Ler = [0 0];                 % [deg] Leading edge radius at each span station
+settings.d = 0.008;                   % [m] rocket tip-fin distance
+settings.zupRaw = 0.0015;             % [m] fin semi-thickness
+settings.LmaxuRaw = 0.0015;           % [m] Fraction of chord from leading edge to max thickness
+settings.Tg = 80;                     % [°C] Carbon fiber glass transition
+settings.Tmax = 100;                  % [°C] carbon fiber max temperature
+
+%%% protub data
+settings.xprot = 1.44;                % [m] axial position from nosecone base
+settings.nloc = 3;                    % [-] number of brakes
+settings.lprot = 0.005;               % [m] brakes thickness
+settings.wprot = 0.088;               % [m] brakes width
+
+%% AIRBRAKES CONTROL SETTINGS
+settings.tControl = settings.tb;                           % [s] time after which the airbrakes can be used
+settings.MachControl = 0.8;                                % [-] Maximum Mach at which airbrakes can be used
+settings.vServo = 150*pi/180;                              % [rad/s] Servo-motor angular velocity
+% First entry, must be always 0 !
+settings.hprot = linspace(0, 0.0387, 3);                   % [m] Block airbrakes opening coordinate
+
+%% MASS GEOMERTY DETAILS
+% x-axis: along the fuselage
+% y-axis: right wing
+% z-axis: downward
+
+% inertias for full configuration (with all the propellant embarqued) obtained with CAD's
+settings.Ixxf = 0.0786;                      % [kg*m^2] Inertia to x-axis
+settings.Iyyf = 12.972;                      % [kg*m^2] Inertia to y-axis
+settings.Izzf = 12.972;                      % [kg*m^2] Inertia to z-axis
+
+% inertias for empty configuration (all the propellant consumed) obtained with CAD's
+settings.Ixxe = 0.0722;                      % [kg*m^2] Inertia to x-axis
+settings.Iyye = 10.041;                      % [kg*m^2] Inertia to y-axis
+settings.Izze = 10.042;                      % [kg*m^2] Inertia to z-axis
+
+%% AERODYNAMICS DETAILS
+% These coefficients are obtained using MISSILE DATCOM
+% (after parsing if using fortran datcom)
+% The files are stored in the ../data/rocketname folder with the following rule:
+% full.mat | empty.mat
+% Relative Path of the data files (default: ../data/). Remember the trailing slash!!
+
+% Coeffs in full.mat is a 5D matrix (6D in empty.mat, extra dimension is aerobrakes extension) 
+% given by Datcom that contains the aerodynamics
+% coefficient computed for the input parameters (Alphas,Betas,Altitudes,Machs)
+% Note: All the parameters (Alphas,Betas,Altitudes,Machs) must be the same for
+% empty and full configuration
+
+% DATA_PATH = '../data/';
+
+filename_coeffs = strcat('SRM_', engineName, '.mat');
+
+if not(isfield(settings,'autoMatProt_flag')) || settings.autoMatProt_flag == false
+    % Coefficients
+    Coeffs = load(filename_coeffs);
+    settings.Coeffs = Coeffs.CoeffsTot;
+
+    % Geometry
+    settings.Geometry = Coeffs.Geometry;
+    settings.Geometry.xcg = Coeffs.Geometry.xcg';
+
+    % State
+    settings.State.Alphas = Coeffs.State.Alphas';
+    settings.State.Betas = Coeffs.State.Betas';
+    settings.State.Altitudes = Coeffs.State.Altitudes';
+    settings.State.Machs = Coeffs.State.Machs';
+    settings.State.Controls = Coeffs.State.hprot';
+    settings.State.xcgTime = Coeffs.State.xcgTime';
+    clear('Coeffs');
+end
+
+
+%% PARACHUTES DETAILS
+settings.NdescentStages = 1; 
+settings.Npara = 2; 
+
+% parachute 1
+settings.para(1, 1).name = 'Rocket DROGUE'; 
+settings.para(1, 1).S = 0.6;       % [m^2]   Surface
+settings.para(1, 1).mass = 0.2;    % [kg]   Parachute Mass
+settings.para(1, 1).CD = 0.75;     % [/] Parachute Drag Coefficient
+settings.para(1, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 1).delay = 1;     % [s] drogue opening delay
+settings.para(1, 1).z_cut = 350;   % [m] Final altitude of the parachute
+
+% parachute 2
+settings.para(2, 1).name = 'Rocket MAIN chute'; 
+settings.para(2, 1).S = 10.54;      % [m^2]   Surface
+settings.para(2, 1).mass = 1.5;    % [kg]   Parachute Mass
+settings.para(2, 1).CD = 0.7;      % [/] Parachute Drag Coefficient
+settings.para(2, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(2, 1).z_cut = 0;     % [m] Final altitude of the parachute
+
+% parachute 1
+settings.para(1, 1).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 1).L = 1.5;         % [m] Shock Chord Length
+settings.para(1, 1).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 1).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 1).Vexit = 5;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 1).CX = 1.2;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(2, 1).L = 6;         % [m] Shock Chord Length
+settings.para(2, 1).K = 3000;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(2, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(2, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(2, 1).nf = 8.7;      % [/] Adimensional Opening Time
+
+%% INTEGRATION OPTIONS
+settings.ode.finalTime =  2000;    % [s] Final integration time
+
+% create an option structure for the integrations:
+
+% - AbsTol is the threshold below which the value of the solution becomes unimportant
+% - RelTol is the tolerance betweeen two consecutive values
+% - Events is the event function that defines when the integration must be
+% - stopped (it has to be created)
+% - InitialStep is the highest value tried by the solver
+
+settings.ode.optionsasc1 = odeset('Events', @eventApogee, 'InitialStep', 1);          %ODE options for ascend
+settings.ode.optionsasc2 = odeset('InitialStep', 1);                                  %ODE options for due to the opening delay of the parachute
+settings.ode.optionsascMultipleAB = odeset('Events', @eventAB,  'InitialStep', 1);    %ODE options for opening of the airbrakes 
+settings.ode.optionspara = odeset('Events', @eventParaCut);                           %ODE options for the parachutes
+settings.ode.optionsdesc = odeset('Events', @eventLanding);                           %ODE options for ballistic descent
+settings.ode.optionspad = odeset('Events', @eventPad);                                %ODE options for the launchpad phase
+
+% Settings for descent 6dof simulation
+settings.ode.optionsDrogue6DOF = odeset('Events', @eventParaCut, 'AbsTol', 1e-6,'RelTol', 1e-6);         %ODE options for due to cutting of the drogue chute
+settings.ode.optionsMainExt6DOF = odeset('Events', @eventMainExit, 'AbsTol', 1e-6,'RelTol', 1e-6);       %ODE options for due to the extraction of the main chute
+settings.ode.optionsMain6DOF = odeset('Events', @eventLanding, 'AbsTol', 1e-6,'RelTol', 1e-6);           %ODE options to terminate descent phase
+
+
+%% STOCHASTIC DETAILS
+%%% launch probability details
+settings.stoch.prob.xLim = 2e3;                    % Max ovest displacement [m]
+settings.stoch.prob.VLim = 50;                     % Max drogue velocity [Pa]
+settings.stoch.prob.XCPLim = 1.5;                  % Min XCP
+
+%%% Safe Ellipse - roccaraso
+settings.prob.SafeEllipse(1).name = 'Prova';
+settings.prob.SafeEllipse(1).a = 1100;
+settings.prob.SafeEllipse(1).b = 2800;
+settings.prob.SafeEllipse(1).x0  = 0;
+settings.prob.SafeEllipse(1).y0 = -300;
+settings.prob.SafeEllipse(1).alpha = 10;
+
+
+%% LANDING POINTS
+% satellite maps of the landing zone 
+% delta limit on the coordinates (for the landing map)
+settings.limLat = 0.04; 
+settings.limLon = 0.025;
+
+%% HRE - compatibility 
+settings.HREmot = false;
+
+%%% efflux data
+settings.motor.Ae = 0; 
+settings.motor.Pe = zeros(size(settings.motor.expTime)); 
+
+%%% intertias 
+settings.Ixx = (settings.motor.expTime./settings.tb).*(settings.Ixxe - settings.Ixxf) + settings.Ixxf; 
+settings.Iyy = (settings.motor.expTime./settings.tb).*(settings.Iyye - settings.Iyyf) + settings.Iyyf;
+settings.Izz = (settings.motor.expTime./settings.tb).*(settings.Izze - settings.Izzf) + settings.Izzf;
+settings.I = [settings.Ixx; settings.Iyy; settings.Izz]; 
+settings.Idot = ones(size(settings.I)).*(settings.I(:, end) - settings.I(:, 1))./settings.tb; 
+
+%%% mass & center of gravity
+settings.mTotalTime = settings.ms + settings.motor.expM; 
+settings.xcg = (settings.motor.expTime./settings.tb).*(settings.xcg(2) - settings.xcg(1)) + settings.xcg(1);
+
+%%% times 
+settings.timeEngineCut = inf; 
diff --git a/missions/2021_Lynx_Roccaraso_September/CAinterpCoeffs.mat b/missions/2021_Lynx_Roccaraso_September/CAinterpCoeffs.mat
new file mode 100644
index 0000000000000000000000000000000000000000..eb37825e739b5dc283c7ce9e591ed61751a6a982
--- /dev/null
+++ b/missions/2021_Lynx_Roccaraso_September/CAinterpCoeffs.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:adcfae194b0f7d75bedfb91b70c7ea4525042c6e377c23518fd48e668f24da0f
+size 456
diff --git a/missions/2021_Lynx_Roccaraso_September/SRM_L1350-CS.mat b/missions/2021_Lynx_Roccaraso_September/SRM_L1350-CS.mat
new file mode 100644
index 0000000000000000000000000000000000000000..43296b54f287f70c5bb6b89f311a6b6b83e007ec
--- /dev/null
+++ b/missions/2021_Lynx_Roccaraso_September/SRM_L1350-CS.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5cbeac2535fe0db37b9d24f88ed614bcfd905dc7d5b1e5d985611c7f6aaee6bd
+size 24177568
diff --git a/missions/2021_Lynx_Roccaraso_September/empty.mat b/missions/2021_Lynx_Roccaraso_September/empty.mat
new file mode 100644
index 0000000000000000000000000000000000000000..d71f3c53d5310b2bdec3ee3ba00f1b1ea6435114
--- /dev/null
+++ b/missions/2021_Lynx_Roccaraso_September/empty.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fac92d733506f69356f3292ff56317cb772ec152fb96c2ca8fec7392859a7777
+size 5963847
diff --git a/missions/2021_Lynx_Roccaraso_September/emptyOld.mat b/missions/2021_Lynx_Roccaraso_September/emptyOld.mat
new file mode 100644
index 0000000000000000000000000000000000000000..8eda3550033802ab3ae5df7aa0bb73fdde7b8b66
--- /dev/null
+++ b/missions/2021_Lynx_Roccaraso_September/emptyOld.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:bea10bde5762dd127da5ba1e0272d4401ad0a8aec23fb7248e957174984c528f
+size 5963847
diff --git a/missions/2021_Lynx_Roccaraso_September/full.mat b/missions/2021_Lynx_Roccaraso_September/full.mat
new file mode 100644
index 0000000000000000000000000000000000000000..5df8e40133a7d89eb1ad16fa1279241d21434fec
--- /dev/null
+++ b/missions/2021_Lynx_Roccaraso_September/full.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fd7153d747e2e034070d9c60b49f7ebdb24dedc3d705c2f12cb8090b3ecd3ad9
+size 1983766
diff --git a/missions/2021_Lynx_Roccaraso_September/fullOld.mat b/missions/2021_Lynx_Roccaraso_September/fullOld.mat
new file mode 100644
index 0000000000000000000000000000000000000000..2d1b0e39f737074279ed6a1508c34fc731c13c72
--- /dev/null
+++ b/missions/2021_Lynx_Roccaraso_September/fullOld.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7a0599dfba3a292302fbe0819fb4f5ee87d7a202897de4499ef27ea95d0d9c5f
+size 1983766
diff --git a/missions/2021_Lynx_Roccaraso_September/simulationsData.m b/missions/2021_Lynx_Roccaraso_September/simulationsData.m
new file mode 100644
index 0000000000000000000000000000000000000000..cc3fafe9859ba7e9c93ee03739d3ef84e0932a5c
--- /dev/null
+++ b/missions/2021_Lynx_Roccaraso_September/simulationsData.m
@@ -0,0 +1,268 @@
+%{
+
+ROCKET DATA FILE
+
+Project name:                         Lynx
+Year of launch:                       2021
+Location of Launch:        Roccaraso (ITA)
+Date of Launch:                 18/09/2021
+
+Last update:                    18/09/2021 
+
+Copyright © 2021, Skyward Experimental Rocketry, AFD department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+% launchpad - roccaraso
+settings.z0 = 1414;                                     % [m] Launchpad Altitude
+settings.lpin = 1.150;                                  % [m] Distance from base of second pin
+settings.lrampa = 5.9 - settings.lpin;                  % [m] LaunchPad route (total available route)
+settings.lat0 = 41.8086605;                             % Launchpad latitude
+settings.lon0 = 14.0543387;                             % Launchpad longitude
+settings.g0 = gravitywgs84(settings.z0, settings.lat0); % Gravity costant at launch latitude and altitude
+
+settings.Local = [settings.z0, 293.15, 85787.047];      % vector [altLocal, TLocal, PLocal, rhoLocal], up to 4 components
+
+%% ENGINE DETAILS
+% load motors data 
+filename = '../Motors.mat';
+load(filename);
+
+engineName = 'L1350-CS';
+crossCase = false;   % true if you want a Cesaroni case with an Aerotech reload, and viceversa
+
+iMotor = [Motors.MotorName] == engineName;
+settings.motor.expTime = Motors(iMotor).t;
+settings.motor.expThrust = Motors(iMotor).T;
+settings.motor.expM = Motors(iMotor).m;
+settings.mp = 1.950;                                    % [kg]   Propellant Mass                                                
+settings.tb = Motors(iMotor).t(end) ;                   % [s]    Burning time
+if not(crossCase)
+    mc = Motors(iMotor).mc;                             % [kg]   Case mass
+else
+    mc = Motors(iMotor).mcCross;
+end
+settings.mNoMot = 17.85 - mc;
+settings.ms = settings.mNoMot + mc;                     % [kg]   Structural Mass
+settings.m0 = settings.ms + settings.mp;                % [kg]   Total Mass
+settings.mnc = 0;                                       % [kg]   Nosecone Mass
+
+%% GEOMETRY/DATCOM DETAILS
+% This parameters should be the same parameters set up in MISSILE DATCOM
+% simulation.
+
+%%% Boat-tail
+% The boat-tail is simulated as a trunc of cone with the major diameter
+% equal to the rocket one (settings.C)
+settings.boatL = 0;                         % [m] Boat-tail length
+settings.boatD = 0;                         % [m] Boat-tail base diameter
+
+%%% Rocket ogive and central body
+settings.C = 0.15;                          % [m] Caliber (Fuselage Diameter)
+settings.S = pi*settings.C^2/4;             % [m^2] Cross-sectional Surface
+settings.xcg = [1.05, 0.98];                % [m] CG postion [full, empty] from nosecone base
+settings.Lnose = 0.26;                      % [m] Nosecone Length
+settings.OgType = 'KARMAN';
+settings.NosePower = 0;                     % [-] Nosecone power type parameter
+settings.rocketLength = 2.270;              % [m] Rocket Length
+settings.Lcenter = settings.rocketLength -...
+    settings.Lnose - settings.boatL;        % [m] fuselage length
+
+%%% finset
+settings.Chord1 = 0.345;                    % [m] attached chord length
+settings.Chord2 = 0.12;                     % [m] free chord length
+settings.Height = 0.117;                    % [m] fin heigth
+settings.deltaXLE = 0.225;                  % [m] start of Chord 2 measured from start of Chord 1
+settings.Npanel = 3;                        % [m] number of fins
+settings.Ler = [0 0];                       % [deg] Leading edge radius at each span station
+settings.d = 0.008;                         % [m] rocket tip-fin distance
+settings.zupRaw = 0.0012;                   % [m] fin semi-thickness
+settings.LmaxuRaw = 0.0012;                 % [m] Fraction of chord from leading edge to max thickness
+
+%%% protub data
+settings.xprot = 1.44;                      % [m] axial position from nosecone base
+settings.nloc = 3;                          % [-] number of brakes
+settings.lprot = 0.008;                     % [m] brakes thickness
+settings.wprot = 0.088;                     % [m] brakes width
+
+%% AIRBRAKES CONTROL SETTINGS
+settings.tControl = 3.7;                                   % [s] time after which the airbrakes can be used
+settings.MachControl = 0.7;                                % [-] Maximum Mach at which airbrakes can be used
+settings.vServo = 150*pi/180;                              % [rad/s] Servo-motor angular velocity
+% First entry, must be always 0 !
+settings.hprot = [0, 0.0264, 0.0387];                      % [m] Block airbrakes opening coordinate
+
+%% MASS GEOMERTY DETAILS
+% x-axis: along the fuselage
+% y-axis: right wing
+% z-axis: downward
+
+% inertias for full configuration (with all the propellant embarqued) obtained with CAD's
+settings.Ixxf = 0.06;                       % [kg*m^2] Inertia to x-axis
+settings.Iyyf = 7.511;                      % [kg*m^2] Inertia to y-axis
+settings.Izzf = 7.512;                      % [kg*m^2] Inertia to z-axis
+
+% inertias for empty configuration (all the propellant consumed) obtained with CAD's
+settings.Ixxe = 0.06;                       % [kg*m^2] Inertia to x-axis
+settings.Iyye = 6.436;                      % [kg*m^2] Inertia to y-axis
+settings.Izze = 6.437;                      % [kg*m^2] Inertia to z-axis
+
+%% AERODYNAMICS DETAILS
+% These coefficients are obtained using MISSILE DATCOM
+% (after parsing if using fortran datcom)
+% The files are stored in the ../data/rocketname folder with the following rule:
+% full.mat | empty.mat
+% Relative Path of the data files (default: ../data/). Remember the trailing slash!!
+
+% Coeffs in full.mat is a 5D matrix (6D in empty.mat, extra dimension is aerobrakes extension) 
+% given by Datcom that contains the aerodynamics
+% coefficient computed for the input parameters (Alphas,Betas,Altitudes,Machs)
+% Note: All the parameters (Alphas,Betas,Altitudes,Machs) must be the same for
+% empty and full configuration
+
+% DATA_PATH = '../data/';
+
+filename_coeffs = strcat('SRM_', engineName, '.mat');
+
+if not(isfield(settings,'autoMatProt_flag')) || settings.autoMatProt_flag == false
+    % Coefficients
+    Coeffs = load(filename_coeffs);
+    settings.Coeffs = Coeffs.CoeffsTot;
+
+    % Geometry
+    settings.Geometry = Coeffs.Geometry;
+    settings.Geometry.xcg = Coeffs.Geometry.xcg';
+
+    % State
+    settings.State.Alphas = Coeffs.State.Alphas';
+    settings.State.Betas = Coeffs.State.Betas';
+    settings.State.Altitudes = Coeffs.State.Altitudes';
+    settings.State.Machs = Coeffs.State.Machs';
+    settings.State.Controls = Coeffs.State.hprot';
+    settings.State.xcgTime = Coeffs.State.xcgTime';
+    clear('Coeffs');
+end
+%% PARACHUTES DETAILS
+
+settings.NdescentStages = 1; 
+settings.Npara = 2; 
+
+% parachute 1
+settings.para(1, 1).name = 'Rocket DROGUE'; 
+settings.para(1, 1).S = 0.6;       % [m^2]   Surface
+settings.para(1, 1).mass = 0.2;    % [kg]   Parachute Mass
+settings.para(1, 1).CD = 0.75;     % [/] Parachute Drag Coefficient
+settings.para(1, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 1).delay = 1;     % [s] drogue opening delay
+settings.para(1, 1).z_cut = 450;   % [m] Final altitude of the parachute
+
+% parachute 2
+settings.para(2, 1).name = 'Rocket MAIN chute'; 
+settings.para(2, 1).S = 10.54;      % [m^2]   Surface
+settings.para(2, 1).mass = 1.5;    % [kg]   Parachute Mass
+settings.para(2, 1).CD = 0.7;      % [/] Parachute Drag Coefficient
+settings.para(2, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(2, 1).z_cut = 0;     % [m] Final altitude of the parachute
+
+% parachute 1
+settings.para(1, 1).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 1).L = 1.5;         % [m] Shock Chord Length
+settings.para(1, 1).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 1).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 1).Vexit = 5;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 1).CX = 1.2;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(2, 1).L = 6;         % [m] Shock Chord Length
+settings.para(2, 1).K = 3000;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(2, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(2, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(2, 1).nf = 8.7;      % [/] Adimensional Opening Time
+
+%% INTEGRATION OPTIONS
+settings.ode.finalTime =  2000;    % [s] Final integration time
+
+% create an option structure for the integrations:
+
+% - AbsTol is the threshold below which the value of the solution becomes unimportant
+% - RelTol is the tolerance betweeen two consecutive values
+% - Events is the event function that defines when the integration must be
+% - stopped (it has to be created)
+% - InitialStep is the highest value tried by the solver
+
+settings.ode.optionsasc1 = odeset('Events', @eventApogee, 'InitialStep', 1);          %ODE options for ascend
+settings.ode.optionsasc2 = odeset('InitialStep', 1);                                  %ODE options for due to the opening delay of the parachute
+settings.ode.optionsascMultipleAB = odeset('Events', @eventAB,  'InitialStep', 1);    %ODE options for opening of the airbrakes 
+settings.ode.optionspara = odeset('Events', @eventParaCut);                           %ODE options for the parachutes
+settings.ode.optionsdesc = odeset('Events', @eventLanding);                           %ODE options for ballistic descent
+settings.ode.optionspad = odeset('Events', @eventPad);                                %ODE options for the launchpad phase
+
+% Settings for descent 6dof simulation
+settings.ode.optionsDrogue6DOF = odeset('Events', @eventParaCut, 'AbsTol', 1e-6,'RelTol', 1e-6);         %ODE options for due to cutting of the drogue chute
+settings.ode.optionsMainExt6DOF = odeset('Events', @eventMainExit, 'AbsTol', 1e-6,'RelTol', 1e-6);       %ODE options for due to the extraction of the main chute
+settings.ode.optionsMain6DOF = odeset('Events', @eventLanding, 'AbsTol', 1e-6,'RelTol', 1e-6);           %ODE options to terminate descent phase
+
+
+%% STOCHASTIC DETAILS
+%%% launch probability details
+settings.stoch.prob.xLim = 2e3;                    % Max ovest displacement [m]
+settings.stoch.prob.VLim = 80;                     % Max drogue velocity [Pa]
+settings.stoch.prob.XCPLim = 1.5;                  % Min XCP
+
+%%% Ellipse - roccaraso
+settings.prob.SafeEllipse(1).name = 'Ellipse';
+settings.prob.SafeEllipse(1).a = 1100;
+settings.prob.SafeEllipse(1).b = 1600;
+settings.prob.SafeEllipse(1).x0  = 100;
+settings.prob.SafeEllipse(1).y0 = -600;
+settings.prob.SafeEllipse(1).alpha = 260;
+
+%%% Safe Ellipse - roccaraso
+settings.prob.SafeEllipse(2).name = 'Safe Ellipse';
+settings.prob.SafeEllipse(2).a = 900;
+settings.prob.SafeEllipse(2).b = 1500;
+settings.prob.SafeEllipse(2).x0  = 100;
+settings.prob.SafeEllipse(2).y0 = -600;
+settings.prob.SafeEllipse(2).alpha = 260;
+
+%%% Very Safe Ellipse - roccaraso
+settings.prob.SafeEllipse(3).name = 'Very Safe Ellipse';
+settings.prob.SafeEllipse(3).a = 700;
+settings.prob.SafeEllipse(3).b = 1000;
+settings.prob.SafeEllipse(3).x0  = 100;
+settings.prob.SafeEllipse(3).y0 = -600;
+settings.prob.SafeEllipse(3).alpha = 260;
+
+
+%% LANDING POINTS
+% satellite maps of the landing zone 
+% delta limit on the coordinates (for the landing map)
+settings.limLat = 0.04; 
+settings.limLon = 0.025;
+
+
+%% HRE - compatibility 
+settings.HREmot = false;
+
+%%% efflux data
+settings.motor.Ae = 0; 
+settings.motor.Pe = zeros(size(settings.motor.expTime)); 
+
+%%% intertias 
+settings.Ixx = (settings.motor.expTime./settings.tb).*(settings.Ixxe - settings.Ixxf) + settings.Ixxf; 
+settings.Iyy = (settings.motor.expTime./settings.tb).*(settings.Iyye - settings.Iyyf) + settings.Iyyf;
+settings.Izz = (settings.motor.expTime./settings.tb).*(settings.Izze - settings.Izzf) + settings.Izzf;
+settings.I = [settings.Ixx; settings.Iyy; settings.Izz]; 
+settings.Idot = ones(size(settings.I)).*(settings.I(:, end) - settings.I(:, 1))./settings.tb; 
+
+%%% mass & center of gravity
+settings.mTotalTime = settings.ms + settings.motor.expM; 
+settings.xcg = (settings.motor.expTime./settings.tb).*(settings.xcg(2) - settings.xcg(1)) + settings.xcg(1);
+
+%%% times 
+settings.timeEngineCut = inf; 
\ No newline at end of file
diff --git a/missions/2022_Pyxis_Portugal_October/CAinterpCoeffs.mat b/missions/2022_Pyxis_Portugal_October/CAinterpCoeffs.mat
new file mode 100644
index 0000000000000000000000000000000000000000..35c1f32148d642acff28292a20fb0ee3954422ea
--- /dev/null
+++ b/missions/2022_Pyxis_Portugal_October/CAinterpCoeffs.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:54e3989ecd832cd0c69affea0c2077c3e58fa5fd87ad148547951cbb4e0dcc53
+size 455
diff --git a/missions/2022_Pyxis_Portugal_October/SRM_M1520-BS.mat b/missions/2022_Pyxis_Portugal_October/SRM_M1520-BS.mat
new file mode 100644
index 0000000000000000000000000000000000000000..ba3ce71a53765bf1a7297ad8581afa4775b1f5ce
--- /dev/null
+++ b/missions/2022_Pyxis_Portugal_October/SRM_M1520-BS.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:e85ad67766f4e9a958d7df544c3edad059e596bc4519a56981f935d83520dc3a
+size 25907517
diff --git a/missions/2022_Pyxis_Portugal_October/empty.mat b/missions/2022_Pyxis_Portugal_October/empty.mat
new file mode 100644
index 0000000000000000000000000000000000000000..4b63319d7b1c7a0bce65f7da733b8fe7172f3730
--- /dev/null
+++ b/missions/2022_Pyxis_Portugal_October/empty.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ff12ab4a04bb9ad22bb2b62eaff083dc1546a3494a93087aff01093be17f8fd5
+size 5977953
diff --git a/missions/2022_Pyxis_Portugal_October/full.mat b/missions/2022_Pyxis_Portugal_October/full.mat
new file mode 100644
index 0000000000000000000000000000000000000000..16d6f2702dcff0fafbfc5655f18c2c4bff9a75e5
--- /dev/null
+++ b/missions/2022_Pyxis_Portugal_October/full.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:afeba525b60fe2c99159ba1efd32b67a6de64aa5a589900b69c558be4f2e57dd
+size 1984810
diff --git a/missions/2022_Pyxis_Portugal_October/simulationsData.m b/missions/2022_Pyxis_Portugal_October/simulationsData.m
new file mode 100644
index 0000000000000000000000000000000000000000..75e03f84c194a2c333b7479bfe7498f8e01473fa
--- /dev/null
+++ b/missions/2022_Pyxis_Portugal_October/simulationsData.m
@@ -0,0 +1,290 @@
+%{
+
+ROCKET DATA FILE
+
+Project name:                        Pyxis
+Year of launch:                       2022
+Location of Launch:      Pont de Sor (POR)
+Date of launch:                        TBD
+
+Last update:                    09/12/2021
+
+Copyright © 2021, Skyward Experimental Rocketry, AFD department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+% launchpad - Pont de Sor
+settings.z0 = 160;                                      % [m] Launchpad Altitude
+settings.lpin = 1.176;                                  % [m] Distance from base of second pin
+settings.lrampa = 12 - settings.lpin;                   % [m] LaunchPad route (total available route)
+settings.lat0 = 39.389733;                              % Launchpad latitude
+settings.lon0 = -8.288992;                              % Launchpad longitude
+settings.g0 = gravitywgs84(settings.z0, settings.lat0); % Gravity costant at launch latitude and altitude
+
+settings.Local = [];                                    % vector [altLocal, TLocal, PLocal, rhoLocal], up to 4 components, the first one is mandatory
+
+%% ENGINE DETAILS
+% load motors data 
+filename = '../Motors.mat';
+load(filename);
+
+engineName = 'M1520-BS';
+crossCase = false;   % true if you want a Cesaroni case with an Aerotech reload, and viceversa
+
+iMotor = [Motors.MotorName] == engineName;
+settings.motor.expTime = Motors(iMotor).t;
+settings.motor.expThrust = Motors(iMotor).T;
+settings.motor.expM = Motors(iMotor).m;
+settings.mp = Motors(iMotor).mp;                  % [kg]   Propellant Mass
+settings.tb = Motors(iMotor).t(end) ;             % [s]    Burning time
+if not(crossCase)
+    mc = Motors(iMotor).mc;                       % [kg]   Case mass
+else
+    mc = Motors(iMotor).mcCross;
+end
+settings.mNoMot = 14;
+settings.ms = settings.mNoMot + mc;               % [kg]   Structural Mass
+settings.m0 = settings.ms + settings.mp;          % [kg]   Total Mass
+settings.mnc = 0;                                 % [kg]   Nosecone Mass
+
+%% GEOMETRY/DATCOM DETAILS
+% This parameters should be the same parameters set up in MISSILE DATCOM
+% simulation.
+
+%%% Boat-tail
+% The boat-tail is simulated as a trunc of cone with the major diameter
+% equal to the rocket one (settings.C)
+settings.boatL = 0.07;                           % [m] Boat-tail length
+settings.boatD = 0.122;                          % [m] Boat-tail base diameter
+
+%%% Rocket ogive and central body
+settings.C = 0.15;                               % [m] Caliber (Fuselage Diameter)
+settings.S = pi*settings.C^2/4;                  % [m^2] Cross-sectional Surface
+settings.xcg = [1.237, 1.099];                   % [m] CG postion [full, empty] from nosecone base
+settings.Lnose = 0.31;                           % [m] Nosecone Length
+settings.OgType = 'MHAACK';                      % [-] Nosecone shape
+settings.NosePower = 0;                          % [-] Nosecone power type parameter
+settings.pMod = 1.293318;                        % [-] P coefficient for modified nosecone shapes
+settings.cMod = 0.1864832;                       % [-] C coefficient for modified nosecone shapes
+settings.Lcenter = 2.166 - settings.boatL;       % [m] fuselage length
+settings.rocketLength = settings.Lcenter + ...
+settings.boatL + settings.Lnose;                 % [m] Rocket Length
+
+%%% Pitot tube
+settings.PITOTCLI = 0;                           % [m] Pitot initial conic section length
+settings.PITOTCLB = 0;                           % [m] Pitot final conic section length
+settings.PITOTCDI = 0;                           % [m] Pitot initial conic section diameter
+settings.PITOTCDB = 0;                           % [m] Pitot final conic section diameter
+settings.PITOTL = 0;                             % [m] Pitot tube length
+settings.PITOTD = 0;                             % [m] Pitot tube diameter
+
+%%% finset 
+settings.Chord1 = 0.34;                          % [m] attached chord length
+settings.Chord2 = 0.10;                          % [m] free chord length
+settings.Height = 0.15;                          % [m] fin heigth
+settings.deltaXLE = 0.21;                        % [m] start of Chord 2 measured from start of Chord 1
+settings.Npanel = 3;                             % [m] number of fins
+settings.Ler = [0 0];                            % [deg] Leading edge radius at each span station
+settings.d = 0.012;                              % [m] rocket tip-fin distance
+settings.zupRaw = 0.00175;                       % [m] fin semi-thickness
+settings.LmaxuRaw = 0.00175;                     % [m] Fraction of chord from leading edge to max thickness
+
+%%% protub data
+settings.xprot = 1.5135;                         % [m] axial position from nosecone base
+settings.nloc = 3;                               % [-] number of brakes
+settings.lprot = 0.008;                          % [m] brakes thickness
+settings.wprot = 0.100;                          % [m] brakes width
+
+%% AIRBRAKES CONTROL SETTINGS
+settings.tControl = settings.tb;                           % [s] time after which the airbrakes can be used
+settings.MachControl = 0.8;                                % [-] Maximum Mach at which airbrakes can be used
+settings.vServo = 150*pi/180;                              % [rad/s] Servo-motor angular velocity
+% First entry, must be always 0 !
+settings.hprot = linspace(0, 0.0373, 3);                   % [m] Block airbrakes opening coordinate
+
+%% DESCENT STAGES
+settings.NdescentStages = 2;                     % [-] Number of descent stages
+% Followimg are for stage 2, 3, ..
+settings.msStages = 4.034;                       % [kg] Mass of stage 2, 3, .. without chutes
+
+%% MASS GEOMERTY DETAILS
+% x-axis: along the fuselage
+% y-axis: right wing
+% z-axis: downward
+
+% inertias for full configuration (with all the propellant embarqued) obtained with CAD's
+settings.Ixxf = 0.06303;                       % [kg*m^2] Inertia to x-axis
+settings.Iyyf = 9.62497;                       % [kg*m^2] Inertia to y-axis
+settings.Izzf = 9.62524;                       % [kg*m^2] Inertia to z-axis
+
+% inertias for empty configuration (all the propellant consumed) obtained with CAD's
+settings.Ixxe = 0.05854;                       % [kg*m^2] Inertia to x-axis
+settings.Iyye = 7.74243;                       % [kg*m^2] Inertia to y-axis
+settings.Izze = 7.74271;                       % [kg*m^2] Inertia to z-axis
+
+%% AERODYNAMICS DETAILS
+% These coefficients are obtained using MISSILE DATCOM
+% (after parsing if using fortran datcom)
+% The files are stored in the ../data/rocketname folder with the following rule:
+% full.mat | empty.mat
+% Relative Path of the data files (default: ../data/). Remember the trailing slash!!
+
+% Coeffs in full.mat is a 5D matrix (6D in empty.mat, extra dimension is aerobrakes extension) 
+% given by Datcom that contains the aerodynamics
+% coefficient computed for the input parameters (Alphas,Betas,Altitudes,Machs)
+% Note: All the parameters (Alphas,Betas,Altitudes,Machs) must be the same for
+% empty and full configuration
+
+% DATA_PATH = '../data/';
+
+filename_coeffs = strcat('SRM_', engineName, '.mat');
+
+if not(isfield(settings,'autoMatProt_flag')) || settings.autoMatProt_flag == false
+    % Coefficients
+    Coeffs = load(filename_coeffs);
+    settings.Coeffs = Coeffs.CoeffsTot;
+
+    % Geometry
+    settings.Geometry = Coeffs.Geometry;
+    settings.Geometry.xcg = Coeffs.Geometry.xcg';
+
+    % State
+    settings.State.Alphas = Coeffs.State.Alphas';
+    settings.State.Betas = Coeffs.State.Betas';
+    settings.State.Altitudes = Coeffs.State.Altitudes';
+    settings.State.Machs = Coeffs.State.Machs';
+    settings.State.Controls = Coeffs.State.hprot';
+    settings.State.xcgTime = Coeffs.State.xcgTime';
+    clear('Coeffs');
+end
+
+%% PARACHUTES DETAILS
+settings.Npara = [2 1];
+
+% parachute 1
+settings.para(1, 1).name = "Rocket DROGUE";
+settings.para(1, 1).S = 0.6;       % [m^2]   Surface
+settings.para(1, 1).mass = 0.2;    % [kg]   Parachute Mass
+settings.para(1, 1).CD = 0.75;     % [/] Parachute Drag Coefficient
+settings.para(1, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 1).delay = 1;     % [s] drogue opening delay
+settings.para(1, 1).z_cut = 450;   % [m] Final altitude of the parachute
+
+% parachute 2
+settings.para(2, 1).name = "Rocket MAIN chute";
+settings.para(2, 1).S = 10.54;     % [m^2]   Surface
+settings.para(2, 1).mass = 1.5;    % [kg]   Parachute Mass
+settings.para(2, 1).CD = 0.7;      % [/] Parachute Drag Coefficient
+settings.para(2, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(2, 1).z_cut = 0;     % [m] Final altitude of the parachute
+
+% parachute 1
+settings.para(1, 1).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 1).L = 1.5;       % [m] Shock Chord Length
+settings.para(1, 1).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 1).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 1).Vexit = 5;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 1).CX = 1.2;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(2, 1).L = 6;         % [m] Shock Chord Length
+settings.para(2, 1).K = 3000;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(2, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(2, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(2, 1).nf = 8.7;      % [/] Adimensional Opening Time
+
+
+%%% PAYLOAD CHUTES
+% parachute 1
+settings.para(1, 2).name = "Payload DROGUE";
+settings.para(1, 2).S = 0.11;      % [m^2]   Surface
+settings.para(1, 2).mass = 0.05;   % [kg]   Parachute Mass
+settings.para(1, 2).CD = 1.2;      % [/] Parachute Drag Coefficient
+settings.para(1, 2).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 2).delay = 1;     % [s] drogue opening delay
+settings.para(1, 2).z_cut = 450;   % [m] Final altitude of the parachute
+settings.para(1, 2).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 2).L = 1.5;       % [m] Shock Chord Length
+settings.para(1, 2).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 2).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 2).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 2).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 2).Vexit = 10;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 2).name = "Payload AIRFOIL";
+settings.para(2, 2).mass = 0.45;    % [kg]   Parachute Mass
+settings.para(2, 2).z_cut = 0;      % [m] Final altitude of the parachute
+
+
+%% INTEGRATION OPTIONS
+settings.ode.finalTime =  2000;    % [s] Final integration time
+
+% create an option structure for the integrations:
+
+% - AbsTol is the threshold below which the value of the solution becomes unimportant
+% - RelTol is the tolerance betweeen two consecutive values
+% - Events is the event function that defines when the integration must be
+% - stopped (it has to be created)
+% - InitialStep is the highest value tried by the solver
+
+settings.ode.optionsasc1 = odeset('Events', @eventApogee, 'InitialStep', 1);          %ODE options for ascend
+settings.ode.optionsasc2 = odeset('InitialStep', 1);                                  %ODE options for due to the opening delay of the parachute
+settings.ode.optionsascMultipleAB = odeset('Events', @eventAB,  'InitialStep', 1);    %ODE options for opening of the airbrakes 
+settings.ode.optionspara = odeset('Events', @eventParaCut);                           %ODE options for the parachutes
+settings.ode.optionsdesc = odeset('Events', @eventLanding);                           %ODE options for ballistic descent
+settings.ode.optionspad = odeset('Events', @eventPad);                                %ODE options for the launchpad phase
+
+% Settings for descent 6dof simulation
+settings.ode.optionsDrogue6DOF = odeset('Events', @eventParaCut, 'AbsTol', 1e-6,'RelTol', 1e-6);         %ODE options for due to cutting of the drogue chute
+settings.ode.optionsMainExt6DOF = odeset('Events', @eventMainExit, 'AbsTol', 1e-6,'RelTol', 1e-6);       %ODE options for due to the extraction of the main chute
+settings.ode.optionsMain6DOF = odeset('Events', @eventLanding, 'AbsTol', 1e-6,'RelTol', 1e-6);           %ODE options to terminate descent phase
+
+
+%% STOCHASTIC DETAILS
+%%% launch probability details
+settings.stoch.prob.xLim = 2e3;                    % Max ovest displacement [m]
+settings.stoch.prob.VLim = 50;                     % Max drogue velocity [Pa]
+settings.stoch.prob.XCPLim = 1.5;                  % Min XCP
+
+%%% Safe Ellipse - roccaraso
+settings.prob.SafeEllipse(1).name = 'prova';
+settings.prob.SafeEllipse(1).a = 1100;
+settings.prob.SafeEllipse(1).b = 2800;
+settings.prob.SafeEllipse(1).x0  = 0;
+settings.prob.SafeEllipse(1).y0 = -300;
+settings.prob.SafeEllipse(1).alpha = 10;
+
+
+%% LANDING POINTS
+% satellite maps of the landing zone 
+% delta limit on the coordinates (for the landing map)
+settings.limLat = 0.04; 
+settings.limLon = 0.025;
+
+
+%% HRE - compatibility 
+settings.HREmot = false;
+
+%%% efflux data
+settings.motor.Ae = 0; 
+settings.motor.Pe = zeros(size(settings.motor.expTime)); 
+
+%%% intertias 
+settings.Ixx = (settings.motor.expTime./settings.tb).*(settings.Ixxe - settings.Ixxf) + settings.Ixxf; 
+settings.Iyy = (settings.motor.expTime./settings.tb).*(settings.Iyye - settings.Iyyf) + settings.Iyyf;
+settings.Izz = (settings.motor.expTime./settings.tb).*(settings.Izze - settings.Izzf) + settings.Izzf;
+settings.I = [settings.Ixx; settings.Iyy; settings.Izz]; 
+settings.Idot = ones(size(settings.I)).*(settings.I(:, end) - settings.I(:, 1))./settings.tb; 
+
+%%% mass & center of gravity
+settings.mTotalTime = settings.ms + settings.motor.expM; 
+settings.xcg = (settings.motor.expTime./settings.tb).*(settings.xcg(2) - settings.xcg(1)) + settings.xcg(1);
+
+%%% times 
+settings.timeEngineCut = inf; 
diff --git a/missions/2022_Pyxis_Roccaraso_September/CAinterpCoeffs.mat b/missions/2022_Pyxis_Roccaraso_September/CAinterpCoeffs.mat
new file mode 100644
index 0000000000000000000000000000000000000000..a4ec1b74175a7921991773cf4ae7dd5737742d61
--- /dev/null
+++ b/missions/2022_Pyxis_Roccaraso_September/CAinterpCoeffs.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:976afed6632e9ac123ce5f4d0532ec90b85127c3d5a473e33e6ba694d629ab53
+size 458
diff --git a/missions/2022_Pyxis_Roccaraso_September/SRM_L820-SK.mat b/missions/2022_Pyxis_Roccaraso_September/SRM_L820-SK.mat
new file mode 100644
index 0000000000000000000000000000000000000000..e858da672b23ee0909ff9f45cb21fc5e358902eb
--- /dev/null
+++ b/missions/2022_Pyxis_Roccaraso_September/SRM_L820-SK.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:421b4a99bb14d477ef2d43d0ddcca24a498669a727e30cb4edf872fad9578c20
+size 24161424
diff --git a/missions/2022_Pyxis_Roccaraso_September/empty.mat b/missions/2022_Pyxis_Roccaraso_September/empty.mat
new file mode 100644
index 0000000000000000000000000000000000000000..f00b420e19b5c064c38e37f60206069dbb6fe617
--- /dev/null
+++ b/missions/2022_Pyxis_Roccaraso_September/empty.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:106a2e15cb9d632a41e2e13121a674f7bd4f8ff0515a115ccf763e416c06beb3
+size 5999853
diff --git a/missions/2022_Pyxis_Roccaraso_September/full.mat b/missions/2022_Pyxis_Roccaraso_September/full.mat
new file mode 100644
index 0000000000000000000000000000000000000000..76291fc51930a7339db6c162dafd5de12b27ec10
--- /dev/null
+++ b/missions/2022_Pyxis_Roccaraso_September/full.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7b4ae7e4f377f148416fe2faf8d2c0f8bd8a99f136821451d9e99c35146ad945
+size 1989864
diff --git a/missions/2022_Pyxis_Roccaraso_September/simulationsData.m b/missions/2022_Pyxis_Roccaraso_September/simulationsData.m
new file mode 100644
index 0000000000000000000000000000000000000000..1eed47747e1ca53c230e68b9bcf10d1e53321b65
--- /dev/null
+++ b/missions/2022_Pyxis_Roccaraso_September/simulationsData.m
@@ -0,0 +1,285 @@
+%{
+
+ROCKET DATA FILE
+
+Project name:                        Pyxis
+Year of launch:                       2022
+Location of Launch:        Roccaraso (ITA)
+Date of launch:                        TBD
+
+Last update:                    20/11/2021
+
+Copyright © 2021, Skyward Experimental Rocketry, AFD department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+% launchpad - Roccaraso
+settings.z0 = 1414;                                     % [m] Launchpad Altitude
+settings.lpin = 1.150;                                  % [m] Distance from base of second pin
+settings.lrampa = 5.9 - settings.lpin;                  % [m] LaunchPad route (total available route)
+settings.lat0 = 41.8084579;                             % Launchpad latitude
+settings.lon0 = 14.0546408;                             % Launchpad longitude
+settings.g0 = gravitywgs84(settings.z0, settings.lat0); % Gravity costant at launch latitude and altitude
+
+settings.Local = [];                                    % vector [altLocal, TLocal, PLocal, rhoLocal], up to 4 components, the first one is mandatory
+
+%% ENGINE DETAILS
+% load motors data 
+filename = '../Motors.mat';
+load(filename);
+
+engineName = 'L820-SK';
+crossCase = false;   % true if you want a Cesaroni case with an Aerotech reload, and viceversa
+
+iMotor = [Motors.MotorName] == engineName;
+settings.motor.expTime = Motors(iMotor).t;
+settings.motor.expThrust = Motors(iMotor).T;
+settings.motor.expM = Motors(iMotor).m;
+settings.mp = Motors(iMotor).mp;                  % [kg]   Propellant Mass
+settings.tb = Motors(iMotor).t(end) ;             % [s]    Burning time
+if not(crossCase)
+    mc = Motors(iMotor).mc;                       % [kg]   Case mass
+else
+    mc = Motors(iMotor).mcCross;
+end
+settings.mNoMot = 13.9;
+settings.ms = settings.mNoMot + mc;               % [kg]   Structural Mass
+settings.m0 = settings.ms + settings.mp;          % [kg]   Total Mass
+settings.mnc = 0;                                 % [kg]   Nosecone Mass
+
+%% GEOMETRY/DATCOM DETAILS
+% This parameters should be the same parameters set up in MISSILE DATCOM
+% simulation.
+
+%%% Boat-tail
+% The boat-tail is simulated as a trunc of cone with the major diameter
+% equal to the rocket one (settings.C)
+settings.boatL = 0.07;                           % [m] Boat-tail length
+settings.boatD = 0.122;                          % [m] Boat-tail base diameter
+
+%%% Rocket ogive and central body
+settings.C = 0.15;                               % [m] Caliber (Fuselage Diameter)
+settings.S = pi*settings.C^2/4;                  % [m^2] Cross-sectional Surface
+settings.xcg = [1.106, 1.013];                   % [m] CG postion [full, empty] from nosecone base
+settings.Lnose = 0.31;                           % [m] Nosecone Length
+settings.OgType = 'MHAACK';                      % [-] Nosecone shape
+settings.NosePower = 0;                          % [-] Nosecone power type parameter
+settings.pMod = 1.293318;                        % [-] P coefficient for modified nosecone shapes
+settings.cMod = 0.1864832;                       % [-] C coefficient for modified nosecone shapes
+settings.Lcenter = 2.166 - settings.boatL;       % [m] fuselage length
+settings.rocketLength = settings.Lcenter + ...
+settings.boatL + settings.Lnose;                 % [m] Rocket Length
+
+%%% Pitot tube
+settings.PITOTCLI = 0;                           % [m] Pitot initial conic section length
+settings.PITOTCLB = 0;                           % [m] Pitot final conic section length
+settings.PITOTCDI = 0;                           % [m] Pitot initial conic section diameter
+settings.PITOTCDB = 0;                           % [m] Pitot final conic section diameter
+settings.PITOTL = 0;                             % [m] Pitot tube length
+settings.PITOTD = 0;                             % [m] Pitot tube diameter
+
+%%% finset 
+settings.Chord1 = 0.34;                          % [m] attached chord length
+settings.Chord2 = 0.10;                          % [m] free chord length
+settings.Height = 0.15;                          % [m] fin heigth
+settings.deltaXLE = 0.21;                        % [m] start of Chord 2 measured from start of Chord 1
+settings.Npanel = 3;                             % [m] number of fins
+settings.Ler = [0 0];                            % [deg] Leading edge radius at each span station
+settings.d = 0.012;                              % [m] rocket tip-fin distance
+settings.zupRaw = 0.00175;                       % [m] fin semi-thickness
+settings.LmaxuRaw = 0.00175;                     % [m] Fraction of chord from leading edge to max thickness
+
+%%% protub data
+settings.xprot = 1.5135;                         % [m] axial position from nosecone base
+settings.nloc = 3;                               % [-] number of brakes
+settings.lprot = 0.008;                          % [m] brakes thickness
+settings.wprot = 0.100;                          % [m] brakes width
+
+%% DESCENT STAGES
+settings.NdescentStages = 2;                     % [-] Number of descent stages
+% Followimg are for stage 2, 3, ..
+settings.msStages = 4.034;                       % [kg] Mass of stage 2, 3, .. without chutes
+
+
+%% AIRBRAKES CONTROL SETTINGS
+settings.tControl = settings.tb;                       % [s] time after which the airbrakes can be used
+settings.MachControl = 0.8;                                % [-] Maximum Mach at which airbrakes can be used
+settings.vServo = 150*pi/180;                              % [rad/s] Servo-motor angular velocity
+% First entry, must be always 0 !
+settings.hprot = linspace(0, 0.0373, 3);                   % [m] Block airbrakes opening coordinate
+
+%% MASS GEOMERTY DETAILS
+% x-axis: along the fuselage
+% y-axis: right wing
+% z-axis: downward
+
+% inertias for full configuration (with all the propellant embarqued) obtained with CAD's
+settings.Ixxf = 0.05669;                       % [kg*m^2] Inertia to x-axis
+settings.Iyyf = 7.61807;                       % [kg*m^2] Inertia to y-axis
+settings.Izzf = 7.61838;                       % [kg*m^2] Inertia to z-axis
+
+% inertias for empty configuration (all the propellant consumed) obtained with CAD's
+settings.Ixxe = 0.05529;                       % [kg*m^2] Inertia to x-axis
+settings.Iyye = 6.42474;                       % [kg*m^2] Inertia to y-axis
+settings.Izze = 6.42506;                       % [kg*m^2] Inertia to z-axis
+
+%% AERODYNAMICS DETAILS
+% These coefficients are obtained using MISSILE DATCOM
+% (after parsing if using fortran datcom)
+% The files are stored in the ../data/rocketname folder with the following rule:
+% full.mat | empty.mat
+% Relative Path of the data files (default: ../data/). Remember the trailing slash!!
+
+% Coeffs in full.mat is a 5D matrix (6D in empty.mat, extra dimension is aerobrakes extension) 
+% given by Datcom that contains the aerodynamics
+% coefficient computed for the input parameters (Alphas,Betas,Altitudes,Machs)
+% Note: All the parameters (Alphas,Betas,Altitudes,Machs) must be the same for
+% empty and full configuration
+
+% DATA_PATH = '../data/';
+
+filename_coeffs = strcat('SRM_', engineName, '.mat');
+
+if not(isfield(settings,'autoMatProt_flag')) || settings.autoMatProt_flag == false
+    % Coefficients
+    Coeffs = load(filename_coeffs);
+    settings.Coeffs = Coeffs.CoeffsTot;
+
+    % Geometry
+    settings.Geometry = Coeffs.Geometry;
+    settings.Geometry.xcg = Coeffs.Geometry.xcg';
+
+    % State
+    settings.State.Alphas = Coeffs.State.Alphas';
+    settings.State.Betas = Coeffs.State.Betas';
+    settings.State.Altitudes = Coeffs.State.Altitudes';
+    settings.State.Machs = Coeffs.State.Machs';
+    settings.State.Controls = Coeffs.State.hprot';
+    settings.State.xcgTime = Coeffs.State.xcgTime';
+    clear('Coeffs');
+end
+%% PARACHUTES DETAILS
+settings.Npara = [2 1];
+
+%%% MAIN ROCKET CHUTES
+% parachute 1
+settings.para(1, 1).name = "Rocket DROGUE";
+settings.para(1, 1).S = 0.4;       % [m^2]   Surface
+settings.para(1, 1).mass = 0.05;    % [kg]   Parachute Mass
+settings.para(1, 1).CD = 0.78;     % [/] Parachute Drag Coefficient
+settings.para(1, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 1).delay = 1;     % [s] drogue opening delay
+settings.para(1, 1).z_cut = 350;   % [m] Final altitude of the parachute
+settings.para(1, 1).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 1).L = 1.5;       % [m] Shock Chord Length
+settings.para(1, 1).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 1).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 1).Vexit = 0;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 1).name = "Rocket MAIN chute";
+settings.para(2, 1).S = 5.5;       % [m^2]   Surface
+settings.para(2, 1).mass = 0.6;    % [kg]   Parachute Mass
+settings.para(2, 1).CD = 0.85;     % [/] Parachute Drag Coefficient
+settings.para(2, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(2, 1).z_cut = 0;     % [m] Final altitude of the parachute
+settings.para(2, 1).CX = 1.2;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(2, 1).L = 6;         % [m] Shock Chord Length
+settings.para(2, 1).K = 3000;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(2, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(2, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(2, 1).nf = 8.7;      % [/] Adimensional Opening Time
+
+%%% PAYLOAD CHUTES
+% parachute 1
+settings.para(1, 2).name = "Payload DROGUE";
+settings.para(1, 2).S = 0.11;      % [m^2]   Surface
+settings.para(1, 2).mass = 0.05;   % [kg]   Parachute Mass
+settings.para(1, 2).CD = 1.2;      % [/] Parachute Drag Coefficient
+settings.para(1, 2).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 2).delay = 1;     % [s] drogue opening delay
+settings.para(1, 2).z_cut = 450;   % [m] Final altitude of the parachute
+settings.para(1, 2).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 2).L = 1.5;       % [m] Shock Chord Length
+settings.para(1, 2).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 2).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 2).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 2).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 2).Vexit = 10;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 2).name = "Payload AIRFOIL";
+settings.para(2, 2).mass = 0.45;    % [kg]   Parachute Mass
+settings.para(2, 2).z_cut = 0;      % [m] Final altitude of the parachute
+
+%% INTEGRATION OPTIONS
+settings.ode.finalTime =  2000;    % [s] Final integration time
+
+% create an option structure for the integrations:
+
+% - AbsTol is the threshold below which the value of the solution becomes unimportant
+% - RelTol is the tolerance betweeen two consecutive values
+% - Events is the event function that defines when the integration must be
+% - stopped (it has to be created)
+% - InitialStep is the highest value tried by the solver
+
+settings.ode.optionsasc1 = odeset('Events', @eventApogee, 'InitialStep', 1);          %ODE options for ascend
+settings.ode.optionsasc2 = odeset('InitialStep', 1);                                  %ODE options for due to the opening delay of the parachute
+settings.ode.optionsascMultipleAB = odeset('Events', @eventAB,  'InitialStep', 1);    %ODE options for opening of the airbrakes 
+settings.ode.optionspara = odeset('Events', @eventParaCut);                           %ODE options for the parachutes
+settings.ode.optionsdesc = odeset('Events', @eventLanding);                           %ODE options for ballistic descent
+settings.ode.optionspad = odeset('Events', @eventPad);                                %ODE options for the launchpad phase
+
+% Settings for descent 6dof simulation
+settings.ode.optionsDrogue6DOF = odeset('Events', @eventParaCut, 'AbsTol', 1e-6,'RelTol', 1e-6);         %ODE options for due to cutting of the drogue chute
+settings.ode.optionsMainExt6DOF = odeset('Events', @eventMainExit, 'AbsTol', 1e-6,'RelTol', 1e-6);       %ODE options for due to the extraction of the main chute
+settings.ode.optionsMain6DOF = odeset('Events', @eventLanding, 'AbsTol', 1e-6,'RelTol', 1e-6);           %ODE options to terminate descent phase
+
+
+%% STOCHASTIC DETAILS
+%%% launch probability details
+settings.stoch.prob.xLim = 2e3;                    % Max ovest displacement [m]
+settings.stoch.prob.VLim = 50;                     % Max drogue velocity [Pa]
+settings.stoch.prob.XCPLim = 1.5;                  % Min XCP
+
+%%% Safe Ellipse - roccaraso
+settings.prob.SafeEllipse(1).name = 'prova';
+settings.prob.SafeEllipse(1).a = 1100;
+settings.prob.SafeEllipse(1).b = 2800;
+settings.prob.SafeEllipse(1).x0  = 0;
+settings.prob.SafeEllipse(1).y0 = -300;
+settings.prob.SafeEllipse(1).alpha = 10;
+
+
+%% LANDING POINTS
+% satellite maps of the landing zone 
+% delta limit on the coordinates (for the landing map)
+settings.limLat = 0.04; 
+settings.limLon = 0.025;
+
+
+%% HRE - compatibility 
+settings.HREmot = false;
+
+%%% efflux data
+settings.motor.Ae = 0; 
+settings.motor.Pe = zeros(size(settings.motor.expTime)); 
+
+%%% intertias 
+settings.Ixx = (settings.motor.expTime./settings.tb).*(settings.Ixxe - settings.Ixxf) + settings.Ixxf; 
+settings.Iyy = (settings.motor.expTime./settings.tb).*(settings.Iyye - settings.Iyyf) + settings.Iyyf;
+settings.Izz = (settings.motor.expTime./settings.tb).*(settings.Izze - settings.Izzf) + settings.Izzf;
+settings.I = [settings.Ixx; settings.Iyy; settings.Izz]; 
+settings.Idot = ones(size(settings.I)).*(settings.I(:, end) - settings.I(:, 1))./settings.tb; 
+
+%%% mass & center of gravity
+settings.mTotalTime = settings.ms + settings.motor.expM; 
+settings.xcg = (settings.motor.expTime./settings.tb).*(settings.xcg(2) - settings.xcg(1)) + settings.xcg(1);
+
+%%% times 
+settings.timeEngineCut = inf; 
diff --git a/missions/2023_Gemini_Portugal_October/CAinterpCoeffs.mat b/missions/2023_Gemini_Portugal_October/CAinterpCoeffs.mat
new file mode 100644
index 0000000000000000000000000000000000000000..32ffc03943f52cd70a072f67f915ebe02d4543db
--- /dev/null
+++ b/missions/2023_Gemini_Portugal_October/CAinterpCoeffs.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5a45b7cfe18472d3149d527626b602e6a31313ba9864d0832f734f9856c6b8c2
+size 458
diff --git a/missions/2023_Gemini_Portugal_October/HRE_FURIA-EUROC-T04T03.mat b/missions/2023_Gemini_Portugal_October/HRE_FURIA-EUROC-T04T03.mat
new file mode 100644
index 0000000000000000000000000000000000000000..27b88aa8a47a71a0c85270eb03e9762adbb32e33
--- /dev/null
+++ b/missions/2023_Gemini_Portugal_October/HRE_FURIA-EUROC-T04T03.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:39019a254a01975fd2e9ef33bc2f99f60f27c3e59fae92ab686f01a46ba0b008
+size 60378278
diff --git a/missions/2023_Gemini_Portugal_October/HRE_FURIA-SFT6.mat b/missions/2023_Gemini_Portugal_October/HRE_FURIA-SFT6.mat
new file mode 100644
index 0000000000000000000000000000000000000000..02c6f1eb4a01d10a03b91b82bd09aff0f83fc0b7
--- /dev/null
+++ b/missions/2023_Gemini_Portugal_October/HRE_FURIA-SFT6.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:bfe5ec74ea1cd080bb24101812563ab3e85328497130f7393d80907f1ef100fa
+size 102505382
diff --git a/missions/2023_Gemini_Portugal_October/HighAOA.mat b/missions/2023_Gemini_Portugal_October/HighAOA.mat
new file mode 100644
index 0000000000000000000000000000000000000000..970ec438cf653e3f8dc07dbd8f1bf031bc5e323f
--- /dev/null
+++ b/missions/2023_Gemini_Portugal_October/HighAOA.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d0d238c8493ecb3d1da1a7cbc3952661d8b34572fb60e879433d4b5df4b94437
+size 3758932
diff --git a/missions/2023_Gemini_Portugal_October/correction.mat b/missions/2023_Gemini_Portugal_October/correction.mat
new file mode 100644
index 0000000000000000000000000000000000000000..d5580214cf4bfaa3110160bedc3fe4a99be0ad2e
--- /dev/null
+++ b/missions/2023_Gemini_Portugal_October/correction.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2f17d0880f3c5805b93db253906cacf9de060d528b057eb9d2f143e8e7715d8e
+size 104840
diff --git a/missions/2023_Gemini_Portugal_October/rocketStress.mat b/missions/2023_Gemini_Portugal_October/rocketStress.mat
new file mode 100644
index 0000000000000000000000000000000000000000..be32b2f2c9eefea083ddd726496755e6a87cae95
--- /dev/null
+++ b/missions/2023_Gemini_Portugal_October/rocketStress.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:702471060e612caee67f825581a61bb46ee51a24027b02cfe3eb143c5c3d4388
+size 1139447
diff --git a/missions/2023_Gemini_Portugal_October/simulationsData.m b/missions/2023_Gemini_Portugal_October/simulationsData.m
new file mode 100644
index 0000000000000000000000000000000000000000..5f62802dcbccea003ff1862ad56da4693c42d8e2
--- /dev/null
+++ b/missions/2023_Gemini_Portugal_October/simulationsData.m
@@ -0,0 +1,369 @@
+%{
+
+ROCKET DATA FILE
+
+Project name:                       Gemini
+Year of launch:                       2023
+Location of Launch:      Pont de Sor (POR)
+Date of launch:                        TBD
+
+Last update:                    09/12/2023
+
+Copyright © 2021, Skyward Experimental Rocketry, AFD department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+% launchpad - Pont de Sor
+settings.z0 = 160;                                            % [m] Launchpad Altitude
+settings.lpin1 = 0.5603;                                      % [m] Distance from the upper pin to the upper tank cap
+settings.lpin2 = 0.2055;                                      % [m] Distance from the lower pin to the lower tank cap  
+settings.lrampa = 12;                                         % [m] Total launchpad length 
+settings.lat0 = 39.388727;                                    % Launchpad latitude
+settings.lon0 = -8.287842;                                    % Launchpad longitude
+settings.g0 = gravitywgs84(settings.z0, settings.lat0);       % Gravity costant at launch latitude and altitude
+
+settings.Local = [];                                          % vector [altLocal, TLocal, PLocal, rhoLocal], up to 4 components
+
+%% ENGINE DETAILS
+% load motors data
+settings.HREmot = true;                                       % non runna con false perche mancano le inerzie totali nel simulationsData
+
+filename = '../HREmotors.mat';
+load(filename);
+
+%engineName = 'HRE_FURIA-EUROC-T04T03';
+engineName = 'HRE_FURIA-SFT6';
+
+iMotor = strcmp({HREmotors.MotorName}, engineName);
+settings.motor.expTime = HREmotors(iMotor).t;                 % [s] Engine time vector
+settings.motor.expThrust = HREmotors(iMotor).T;               % [N] Engine thrust vector
+settings.motor.expM = HREmotors(iMotor).m;                    % [Kg] Propellant Mass (in time)
+settings.motor.L = HREmotors(iMotor).L;                       % [m] Engine length
+settings.motor.Ltank = HREmotors(iMotor).Ltank;               % [m] Tank length
+settings.motor.xcg = HREmotors(iMotor).xcg;                   % [m] Engine xcg from tank tip
+settings.motor.mc = HREmotors(iMotor).mc;                     % [kg] Engine Structural Mass
+settings.motor.Pe = HREmotors(iMotor).Pe;                     % [Pa] Eflux pressure
+settings.motor.Ae = HREmotors(iMotor).Ae;                     % [Pa] Eflux Area
+settings.motor.mFuel = HREmotors(iMotor).mFu;                 % [kg] Fuel (grain) initial mass
+settings.motor.mOx = HREmotors(iMotor).mOx;                   % [kg] Oxidizer initial mass
+settings.mp = HREmotors(iMotor).mp;                           % [kg] Initial propellant Mass
+settings.tb = HREmotors(iMotor).t(end);                       % [s]  Burning time
+
+% Engine fuselage data
+settings.motor.mFuselage = HREmotors(iMotor).mFus;            % [kg] Fuselage of the engine only
+settings.motor.xcgFuselage = (settings.motor.L - ...
+settings.motor.Ltank)/2 + settings.motor.Ltank;               % [m]  xcg of the engine fuselage only from tank tip
+
+% Engine cut
+settings.timeEngineCut = inf;                                 % [s] Moment in time in wich the engine will be cut off
+settings.tIGN = 0.4;                                          % [s] ignition transient
+settings.tCO = 0.3;                                           % [s] cut-off transient
+
+%% Launchpad synthesis
+
+% settings.lpin = settings.lpin1 + settings.motor.Ltank ... 
+%                    + settings.lpin2;                         % [m] Distance between the two launchpad pins
+
+settings.lpin = 2.4 + 0.2;                                    % [m] distance of the upper pin from the rail base (upper pin-boat + boat-rail base)
+settings.lrampa = settings.lrampa - settings.lpin;            % [m] LaunchPad route (total available route)
+
+%% ROCKET DETAILS
+%%% Boat-tail
+% The boat-tail is simulated as a trunc of cone with the major diameter
+% equal to the rocket one (settings.C)
+settings.boatL = 0.07;                                        % [m] Boat-tail length
+settings.boatD = 0.122;                                       % [m] Boat-tail base diameter
+
+%%% Rocket ogive and central body
+settings.C = 0.150;                                           % [m] Caliber (Fuselage Diameter)
+settings.S = pi*settings.C^2/4;                               % [m^2] Cross-sectional Surface
+settings.Lnose = 0.32;                                        % [m] Nosecone Length
+settings.OgType = 'MHAACK';                                   % [-] Nosecone shape
+settings.NosePower = 3/4;                                     % [-] Nosecone power type parameter
+settings.pMod = 1.250152e+00;                                 % [-] P coefficient for modified nosecone shapes
+settings.cMod = 1.799127e-01;                                 % [-] C coefficient for modified nosecone shapes
+
+settings.LcenterRocket = 1.517;                               % [m] Central body length
+settings.mNoMot = 16.423;                                       % [kg] Mass of everything above engine
+settings.xcgNoMot = 1.254;                                    % [m] XCG of everything above engine (from nosecone base)
+
+% SYNTHESIS WITH ENGINE FUSELAGE MASS AND XCG SEPARATELY GIVEN
+settings.addMfusMot = false;
+
+if settings.addMfusMot
+    settings.xcgNoMot = (settings.mNoMot*settings.xcgNoMot + settings.motor.mFuselage*(settings.motor.xcgFuselage + settings.LcenterRocket))...
+                        /(settings.mNoMot + settings.motor.mFuselage);
+    settings.mNoMot = settings.mNoMot + settings.motor.mFuselage;
+end
+
+%% SYNTHESIS
+settings.Lcenter = settings.LcenterRocket + settings.motor.L; % [m] Total centerbody length
+settings.mnc = 0;                                             % [kg] Nosecone Mass
+settings.ms = settings.mNoMot + settings.motor.mc + ...
+settings.mnc;                                                 % [kg] Total structural Mass
+settings.m0 = settings.ms + settings.motor.expM(1);           % [kg] Total initial Mass
+settings.mTotalTime = settings.mNoMot + settings.motor.mc + ...
+settings.mnc + settings.motor.expM;                           % Total mass (in time)
+
+% Total XCG of the rocket [m] !! FROM NOSECONE BASE !!
+settings.xcg = ((settings.motor.expM + settings.motor.mc) .* (settings.LcenterRocket + settings.motor.xcg) + ...
+    settings.mNoMot * settings.xcgNoMot) ./ (settings.motor.expM + settings.motor.mc + settings.mNoMot);
+
+settings.Lcenter = settings.Lcenter - settings.boatL;         % correcting the rocket length in order to account the boat tail
+
+%% AERODYNAMICS
+%%% Pitot tube
+settings.PITOTCLI = 0;                                        % [m] Pitot initial conic section length
+settings.PITOTCLB = 0;                                        % [m] Pitot final conic section length
+settings.PITOTCDI = 0;                                        % [m] Pitot initial conic section diameter
+settings.PITOTCDB = 0;                                        % [m] Pitot final conic section diameter
+settings.PITOTL = 0;                                          % [m] Pitot tube length
+settings.PITOTD = 0;                                          % [m] Pitot tube diameter
+
+%%% finset
+settings.Chord1 = 0.30;                                       % [m] attached chord length
+settings.Chord2 = 0.14;                                       % [m] free chord length
+settings.Height = 0.11;                                       % [m] fin heigth
+settings.deltaXLE = 0.13;                                     % [m] start of Chord 2 measured from start of Chord 1
+settings.Npanel = 3;                                          % [m] number of fins
+settings.Ler = [0 0];                                         % [deg] Leading edge radius at each span station
+settings.d = 0.012;                                           % [m] distance between end of root chord and end of center body
+settings.zupRaw = 0.00175;                                    % [m] fin semi-thickness
+settings.LmaxuRaw = 0.00175;                                  % [m] Fraction of chord from leading edge to max thickness
+
+%%% protub data
+settings.xprot = 1.517;                                       % [m] axial position from nosecone base
+settings.nloc = 3;                                            % [-] number of brakes
+settings.lprot = 0.008;                                       % [m] brakes thickness
+settings.wprot = 0.1002754821;                                % [m] brakes width (normal)
+
+%% AIRBRAKES CONTROL SETTINGS
+settings.tControl = settings.tb;                              % [s] time after which the airbrakes can be used
+settings.MachControl = 0.8;                                   % [-] Maximum Mach at which airbrakes can be used
+settings.vServo = 150*pi/180;                                 % [rad/s] Servo-motor angular velocity
+% First entry, must be always 0 !
+settings.hprot = linspace(0, 0.0363, 3);                      % [m] Block airbrakes opening coordinate
+
+%% MASS GEOMERTY DETAILS
+% x-axis: along the fuselage
+% y-axis: right wing
+% z-axis: downward
+
+%%% Inertias for fuselage only (no engine)
+settings.rocketIxx = 0.06535397;                              % [kg*m^2] Inertia to x-axis
+settings.rocketIyy = 17.21019828;                             % [kg*m^2] Inertia to y-axis
+settings.rocketIzz = 17.21056483;                             % [kg*m^2] Inertia to z-axis
+
+%%% Inertias for engine only
+engineIxx = HREmotors(iMotor).Ixx;                            % [kg*m^2] Inertia to x-axis
+engineIyy = HREmotors(iMotor).Iyy;                            % [kg*m^2] Inertia to y-axis
+engineIzz = HREmotors(iMotor).Izz;                            % [kg*m^2] Inertia to z-axis
+
+%%% Total inertias
+settings.Ixx = settings.rocketIxx + engineIxx;
+
+settings.Iyy = settings.rocketIyy + ((settings.xcgNoMot - settings.xcg).^2).*settings.mNoMot + ...
+    (engineIyy + ( (settings.motor.xcg + settings.LcenterRocket - settings.xcg).^2 ) .*(settings.motor.expM + settings.motor.mc));
+
+settings.Izz = settings.rocketIzz + ((settings.xcgNoMot - settings.xcg).^2).*settings.mNoMot + ...
+    (engineIzz + ( (settings.motor.xcg + settings.LcenterRocket - settings.xcg).^2).*(settings.motor.expM + settings.motor.mc));
+
+settings.I = [settings.Ixx; settings.Iyy; settings.Izz];
+
+%%% Inertias derivatives (avoiding calculations in ascent.m)
+settings.Idot = diff(settings.I')'./diff(settings.motor.expTime);
+settings.Idot(:, end+1) = settings.Idot(:, end);
+
+%% APPLY CORRECTIONS
+
+if exist("correction.mat","file")
+    corr = load("correction.mat");
+    settings.mTotalTime = settings.mTotalTime + corr.deltaM;
+    settings.xcg = settings.xcg + reshapeVector(settings.motor.expTime, corr.time, corr.deltaXcg);
+    settings.ms = settings.ms + corr.deltaM;
+end
+
+%% AERODYNAMICS DETAILS
+% These coefficients are obtained using MISSILE DATCOM, check simulator/README.md for further details.
+
+% DATA_PATH = '../data/';
+
+filename_coeffs = strcat(engineName, '.mat');
+
+if not(isfield(settings,'autoMatProt_flag')) || settings.autoMatProt_flag == false
+    % Coefficients
+    Coeffs = load(filename_coeffs);
+    settings.Coeffs = Coeffs.CoeffsTot;
+
+    % Geometry
+    settings.Geometry = Coeffs.Geometry;
+    settings.Geometry.xcg = Coeffs.Geometry.xcg';
+
+    % State
+    settings.State.Alphas = Coeffs.State.Alphas';
+    settings.State.Betas = Coeffs.State.Betas';
+    settings.State.Altitudes = Coeffs.State.Altitudes';
+    settings.State.Machs = Coeffs.State.Machs';
+    settings.State.Controls = Coeffs.State.hprot';
+    settings.State.xcgTime = Coeffs.State.xcgTime';
+    clear('Coeffs');
+
+    CoeffHighAOA = load('HighAOA.mat'); 
+    settings.highAOA.Coeffs = CoeffHighAOA.CoeffsTot; 
+    settings.highAOA.State.Alphas = CoeffHighAOA.State.Alphas';
+    settings.highAOA.State.Betas = CoeffHighAOA.State.Betas'; 
+    settings.highAOA.State.Altitudes = CoeffHighAOA.State.Altitudes'; 
+    settings.highAOA.State.Machs = CoeffHighAOA.State.Machs';
+    settings.highAOA.State.Controls = CoeffHighAOA.State.hprot';
+    settings.highAOA.State.xcgTime = CoeffHighAOA.State.xcgTime';
+    clear('CoeffHighAOA');
+
+    settings.highAOA.tb = settings.tb; 
+    settings.highAOA.timeEngineCut = settings.timeEngineCut; 
+
+end
+
+%% DESCENT STAGES
+settings.NdescentStages = 2;                                  % [-] Number of descent stages
+% Followimg are for stage 2, 3, ..
+settings.msStages = 4.412;                                    % [kg] Mass of stage 2, 3, .. without chutes
+
+%% PARACHUTES DETAILS
+
+settings.Npara = [2 1];
+
+% parachute 1
+settings.para(1, 1).name = 'DROGUE chute';
+settings.para(1, 1).S = 1.2219;    % [m^2]   Surface
+settings.para(1, 1).mass = 0.15;   % [kg]   Parachute Mass
+settings.para(1, 1).CD = 0.96;     % [/] Parachute Drag Coefficient
+settings.para(1, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 1).delay = 1;     % [s] drogue opening delay
+settings.para(1, 1).z_cut = 350;   % [m] Final altitude of the parachute
+
+% parachute 2
+settings.para(2, 1).name = 'MAIN chute';
+settings.para(2, 1).S = 7.34;      % [m^2]   Surface
+settings.para(2, 1).mass = 1.05;    % [kg]   Parachute Mass
+settings.para(2, 1).CD = 2.11;     % [/] Parachute Drag Coefficient
+settings.para(2, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(2, 1).z_cut = 0;     % [m] Final altitude of the parachute
+
+% parachute 1
+settings.para(1, 1).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 1).L = 1.5;       % [m] Shock Chord Length
+settings.para(1, 1).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 1).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 1).Vexit = 5;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 1).CX = 1.2;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(2, 1).L = 6;         % [m] Shock Chord Length
+settings.para(2, 1).K = 3000;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(2, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(2, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(2, 1).nf = 8.7;      % [/] Adimensional Opening Time
+
+%%% PAYLOAD CHUTES
+% parachute 1
+settings.para(1, 2).name = "Payload DROGUE";
+settings.para(1, 2).S = 0.11;      % [m^2]   Surface
+settings.para(1, 2).mass = 0.15;   % [kg]   Parachute Mass
+settings.para(1, 2).CD = 1.2;      % [/] Parachute Drag Coefficient
+settings.para(1, 2).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 2).delay = 1;     % [s] drogue opening delay
+settings.para(1, 2).z_cut = 450;   % [m] Final altitude of the parachute
+settings.para(1, 2).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 2).L = 1.5;       % [m] Shock Chord Length
+settings.para(1, 2).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 2).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 2).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 2).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 2).Vexit = 10;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 2).name = "Payload AIRFOIL";
+settings.para(2, 2).mass = 0.50;    % [kg]   Parachute Mass
+settings.para(2, 2).z_cut = 0;      % [m] Final altitude of the parachute
+
+
+%% PARAFOIL DETAILS
+% Geometry
+settings.payload.paraPrev = [2, 1];         % Reference to the previous parachute in the multistage descent 
+                                            % (in this case, the payload drogue: 2nd stage, 1st parachute)
+
+settings.payload.mass = 4.2;                 % [kg]  mass 
+settings.payload.b    = 2.55/2;              % [m]   semiwingspan  - vela nuova: 2.55/2; - vela vecchia: 2.06/2;
+settings.payload.c    = 0.8;                 % [m]   mean aero chord
+settings.payload.S    = 2.04;                % [m^2] payload surface - vela nuova 2.04; - vela vecchia: 1.64;  
+settings.payload.inertia = [0.42, 0,   0.03;
+                            0,    0.4,    0; 
+                            0.03, 0, 0.053]; % [kg m^2] [3x3] inertia matrix payload+payload 
+
+settings.payload.inverseInertia = inv(settings.payload.inertia);
+
+% Aerodynamic constants
+settings.payload.CD0       =  0.25; 
+settings.payload.CDAlpha2  =  0.12;
+settings.payload.CL0       =  0.091;
+settings.payload.CLAlpha   =  0.9;
+settings.payload.Cm0       =  0.35; 
+settings.payload.CmAlpha   = -0.72;
+settings.payload.Cmq       = -1.49;
+settings.payload.CLDeltaA  = -0.0035;
+settings.payload.Cnr       = -0.27;
+settings.payload.CnDeltaA  =  0.0115;
+settings.payload.deltaSMax =  0.1;
+settings.payload.CDDeltaA  = 0.01;
+settings.payload.Clp       = -0.84;
+settings.payload.ClPhi     = -0.1;
+settings.payload.ClDeltaA  = 0.01;
+
+%% INTEGRATION OPTIONS
+settings.ode.finalTime =  2000;    % [s] Final integration time
+
+% create an option structure for the integrations:
+
+% - AbsTol is the threshold below which the value of the solution becomes unimportant
+% - RelTol is the tolerance betweeen two consecutive values
+% - Events is the event function that defines when the integration must be
+% - stopped (it has to be created)
+% - InitialStep is the highest value tried by the solver
+
+settings.ode.optionsasc1 = odeset('Events', @eventApogee, 'InitialStep', 1);          %ODE options for ascend
+settings.ode.optionsasc2 = odeset('InitialStep', 1);                                  %ODE options for due to the opening delay of the parachute
+settings.ode.optionsascMultipleAB = odeset('Events', @eventAB,  'InitialStep', 1);    %ODE options for opening of the airbrakes
+settings.ode.optionspara = odeset('Events', @eventParaCut);                           %ODE options for the parachutes
+settings.ode.optionsdesc = odeset('Events', @eventLanding,...
+                                  'RelTol', 1e-3, 'AbsTol', 1e-3);                    %ODE options for ballistic descent
+settings.ode.optionspad = odeset('Events', @eventPad);                                %ODE options for the launchpad phase
+
+% Settings for descent 6dof simulation
+settings.ode.optionsDrogue6DOF = odeset('Events', @eventParaCut, 'AbsTol', 1e-6,'RelTol', 1e-6);         %ODE options for due to cutting of the drogue chute
+settings.ode.optionsMainExt6DOF = odeset('Events', @eventMainExit, 'AbsTol', 1e-6,'RelTol', 1e-6);       %ODE options for due to the extraction of the main chute
+settings.ode.optionsMain6DOF = odeset('Events', @eventLanding, 'AbsTol', 1e-6,'RelTol', 1e-6);           %ODE options to terminate descent phase
+
+%% STOCHASTIC DETAILS
+%%% launch probability details
+settings.stoch.prob.xLim = 2e3;                    % Max ovest displacement [m]
+settings.stoch.prob.VLim = 50;                     % Max drogue velocity [Pa]
+settings.stoch.prob.XCPLim = 1.5;                  % Min XCP
+
+%%% Safe Ellipse - roccaraso
+settings.prob.SafeEllipse(1).name = 'prova';
+settings.prob.SafeEllipse(1).a = 1100;
+settings.prob.SafeEllipse(1).b = 2800;
+settings.prob.SafeEllipse(1).x0  = 0;
+settings.prob.SafeEllipse(1).y0 = -300;
+settings.prob.SafeEllipse(1).alpha = 10;
+
+%% LANDING POINTS
+% satellite maps of the landing zone
+% delta limit on the coordinates (for the landing map)
+settings.limLat = 0.04;
+settings.limLon = 0.025;
diff --git a/missions/2023_Gemini_Roccaraso_September/CAinterpCoeffs.mat b/missions/2023_Gemini_Roccaraso_September/CAinterpCoeffs.mat
new file mode 100644
index 0000000000000000000000000000000000000000..1a3167d1fe0e92f3104651e7ee85982890eb4def
--- /dev/null
+++ b/missions/2023_Gemini_Roccaraso_September/CAinterpCoeffs.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:808d0a54031039114c060054f3c7f7b78a6c605bd1e6de58426ed97bb7625c75
+size 457
diff --git a/missions/2023_Gemini_Roccaraso_September/HRE_FURIA-R-SFT6.mat b/missions/2023_Gemini_Roccaraso_September/HRE_FURIA-R-SFT6.mat
new file mode 100644
index 0000000000000000000000000000000000000000..b4fe7f40100a274ac7ef733e753e24917d61d0b8
--- /dev/null
+++ b/missions/2023_Gemini_Roccaraso_September/HRE_FURIA-R-SFT6.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6aef0549361528025e4525ad22a9938063135c5467b5a93e6e8e1661789868ba
+size 61487003
diff --git a/missions/2023_Gemini_Roccaraso_September/HRE_FURIA-Rv1-T04T03.mat b/missions/2023_Gemini_Roccaraso_September/HRE_FURIA-Rv1-T04T03.mat
new file mode 100644
index 0000000000000000000000000000000000000000..be040a693c7ee2211c684aa16b77184cff959ce6
--- /dev/null
+++ b/missions/2023_Gemini_Roccaraso_September/HRE_FURIA-Rv1-T04T03.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:051848cbf7eccd159de5adf816c7de26b39bcc3c8a7a4a19dfe9bfdc48a185da
+size 60384667
diff --git a/missions/2023_Gemini_Roccaraso_September/HRE_FURIA-Rv2-T04T03.mat b/missions/2023_Gemini_Roccaraso_September/HRE_FURIA-Rv2-T04T03.mat
new file mode 100644
index 0000000000000000000000000000000000000000..1b0539df112d175470fd64a10a1dc5966190fbff
--- /dev/null
+++ b/missions/2023_Gemini_Roccaraso_September/HRE_FURIA-Rv2-T04T03.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9c59bce16eeac7f1575d5d77786dd4ebdd7204b79421ebaec9798ddc13b478f7
+size 60388201
diff --git a/missions/2023_Gemini_Roccaraso_September/HighAOA.mat b/missions/2023_Gemini_Roccaraso_September/HighAOA.mat
new file mode 100644
index 0000000000000000000000000000000000000000..3eb85734ff2af661f8ca372a7001fc92c1e83355
--- /dev/null
+++ b/missions/2023_Gemini_Roccaraso_September/HighAOA.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:34c11f637bb1e61172650442e662d5c781ca39a8ef50c90a44e9558672a50847
+size 2036925
diff --git a/missions/2023_Gemini_Roccaraso_September/correction.mat b/missions/2023_Gemini_Roccaraso_September/correction.mat
new file mode 100644
index 0000000000000000000000000000000000000000..7309b6986781dfd4e983cb133d1c59095d8c34ff
--- /dev/null
+++ b/missions/2023_Gemini_Roccaraso_September/correction.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:579c4f01f5575c48f2a94f7111e23086f0f5dd2b84f5fc6a82dee17d8ed6101d
+size 164586
diff --git a/missions/2023_Gemini_Roccaraso_September/rocketStress.mat b/missions/2023_Gemini_Roccaraso_September/rocketStress.mat
new file mode 100644
index 0000000000000000000000000000000000000000..6d05191cd837b186ff576e1563fc6cd9530bc143
--- /dev/null
+++ b/missions/2023_Gemini_Roccaraso_September/rocketStress.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:1ee95ae166010ba5926a705841239d54958e20b99ee1a73b82b70c79d4c327a4
+size 773137
diff --git a/missions/2023_Gemini_Roccaraso_September/simulationsData.m b/missions/2023_Gemini_Roccaraso_September/simulationsData.m
new file mode 100644
index 0000000000000000000000000000000000000000..9ad12523b74895fe9301ade2c6f86ae07c6eb9f0
--- /dev/null
+++ b/missions/2023_Gemini_Roccaraso_September/simulationsData.m
@@ -0,0 +1,335 @@
+%{
+
+ROCKET DATA FILE
+
+Project name:                       Gemini
+Year of launch:                       2023
+Location of Launch:       Roccasarso (ITA)
+Date of launch:                        TBD
+
+Last update:                    09/12/2023
+
+Copyright © 2021, Skyward Experimental Rocketry, AFD department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+% launchpad - Roccaraso
+settings.z0 = 1414;                                           % [m] Launchpad Altitude
+settings.lpin1 = 0.5603;                                      % [m] Distance from the upper pin to the upper tank cap
+settings.lpin2 = 0.2055;                                      % [m] Distance from the lower pin to the lower tank cap  
+settings.lrampa = 10;                                         % [m] Total launchpad length 
+settings.lat0 = 41.8084579;                                   % Launchpad latitude
+settings.lon0 = 14.0546408;                                   % Launchpad longitude
+settings.g0 = gravitywgs84(settings.z0, settings.lat0);       % Gravity costant at launch latitude and altitude
+
+settings.Local = [];            %% ENGINE DETAILS
+
+%% load motors data
+settings.HREmot = true;                                       % non runna con false perche mancano le inerzie totali nel simulationsData
+
+filename = '../HREmotors.mat';
+load(filename);
+
+% engineName = 'HRE_FURIA-Rv2-T04T03'; 
+engineName = 'HRE_FURIA-R-SFT6';
+
+iMotor = strcmp({HREmotors.MotorName}, engineName);
+settings.motor.expTime = HREmotors(iMotor).t;                 % [s] Engine time vector
+settings.motor.expThrust = HREmotors(iMotor).T;               % [N] Engine thrust vector
+settings.motor.expM = HREmotors(iMotor).m;                    % [Kg] Propellant Mass (in time)
+settings.motor.L = HREmotors(iMotor).L;                       % [m] Engine length
+settings.motor.Ltank = HREmotors(iMotor).Ltank;               % [m] Tank length
+settings.motor.xcg = HREmotors(iMotor).xcg;                   % [m] Engine xcg from tank tip
+settings.motor.mc = HREmotors(iMotor).mc;                     % [kg] Engine Structural Mass
+settings.motor.Pe = HREmotors(iMotor).Pe;                     % [Pa] Eflux pressure
+settings.motor.Ae = HREmotors(iMotor).Ae;                     % [Pa] Eflux Area
+settings.motor.mFuel = HREmotors(iMotor).mFu;                 % [kg] Fuel (grain) initial mass
+settings.motor.mOx = HREmotors(iMotor).mOx;                   % [kg] Oxidizer initial mass
+settings.mp = HREmotors(iMotor).mp;                           % [kg] Initial propellant Mass
+settings.tb = HREmotors(iMotor).t(end) ;                      % [s]  Burning time
+
+% Engine fuselage data
+settings.motor.mFuselage = HREmotors(iMotor).mFus;            % [kg] Fuselage of the engine only
+settings.motor.xcgFuselage = (settings.motor.L - ...
+settings.motor.Ltank)/2 + settings.motor.Ltank;               % [m]  xcg of the engine fuselage only from tank tip
+
+% Engine cut
+settings.timeEngineCut = inf;                                 % [s] Moment in time in wich the engine will be cut off
+settings.tIGN = 0.4;                                          % [s] ignition transient
+settings.tCO = 0.3;                                           % [s] cut-off transient
+
+%% Launchpad synthesis
+
+%settings.lpin = settings.lpin1 + settings.motor.Ltank ... 
+%                    + settings.lpin2;                        % [m] Distance between the two launchpad pins
+
+settings.lpin = 1.926 + 0.2;                                  % [m] distance of the upper pin from the rail base (upper pin-boat + boat-rail base)
+settings.lrampa = settings.lrampa - settings.lpin;            % [m] LaunchPad route (total available route)
+
+%% ROCKET DETAILS
+%%% Boat-tail
+% The boat-tail is simulated as a trunc of cone with the major diameter
+% equal to the rocket one (settings.C)
+settings.boatL = 0.07;                                        % [m] Boat-tail length
+settings.boatD = 0.122;                                       % [m] Boat-tail base diameter
+
+%%% Rocket ogive and central body
+settings.C = 0.150;                                           % [m] Caliber (Fuselage Diameter)
+settings.S = pi*settings.C^2/4;                               % [m^2] Cross-sectional Surface
+settings.Lnose = 0.32;                                        % [m] Nosecone Length
+settings.OgType = 'MHAACK';                                   % [-] Nosecone shape
+settings.NosePower = 3/4;                                     % [-] Nosecone power type parameter
+settings.pMod = 1.250152e+00;                                 % [-] P coefficient for modified nosecone shapes
+settings.cMod = 1.799127e-01;                                 % [-] C coefficient for modified nosecone shapes
+
+settings.LcenterRocket = 1.517;                               % [m] Central body length
+settings.mNoMot = 16.423;                                       % [kg] Mass of everything above engine
+settings.xcgNoMot = 1.149;                                    % [m] XCG of everything above engine (from nosecone base)
+
+%SYNTHESIS WITH ENGINE FUSELAGE MASS AND XCG SEPARATELY GIVEN
+settings.addMfusMot = false;
+
+if settings.addMfusMot
+    settings.xcgNoMot = (settings.mNoMot*settings.xcgNoMot + settings.motor.mFuselage*(settings.motor.xcgFuselage + settings.LcenterRocket))...
+                        /(settings.mNoMot + settings.motor.mFuselage);
+    settings.mNoMot = settings.mNoMot + settings.motor.mFuselage;
+end
+
+%% SYNTHESIS
+settings.Lcenter = settings.LcenterRocket + settings.motor.L; % [m] Total centerbody length
+settings.mnc = 0;                                             % [kg] Nosecone Mass
+settings.ms = settings.mNoMot + settings.motor.mc + ...
+settings.mnc;                                                 % [kg] Total structural Mass
+settings.m0 = settings.ms + settings.motor.expM(1);           % [kg] Total initial Mass
+settings.mTotalTime = settings.mNoMot + settings.motor.mc + ...
+settings.mnc + settings.motor.expM;                           % Total mass (in time)
+
+% Total XCG of the rocket [m] !! FROM NOSECONE BASE !!
+settings.xcg = ( (settings.motor.expM + settings.motor.mc).*(settings.LcenterRocket + settings.motor.xcg) + ...
+    settings.mNoMot*settings.xcgNoMot)./(settings.motor.expM + settings.motor.mc + settings.mNoMot);
+
+settings.Lcenter = settings.Lcenter - settings.boatL;         % correcting the rocket length in order to account the boat tail
+
+%% AERODYNAMICS
+%%% Pitot tube
+settings.PITOTCLI = 0;                                        % [m] Pitot initial conic section length
+settings.PITOTCLB = 0;                                        % [m] Pitot final conic section length
+settings.PITOTCDI = 0;                                        % [m] Pitot initial conic section diameter
+settings.PITOTCDB = 0;                                        % [m] Pitot final conic section diameter
+settings.PITOTL = 0;                                          % [m] Pitot tube length
+settings.PITOTD = 0;                                          % [m] Pitot tube diameter
+
+%%% finset 
+settings.Chord1 = 0.30;                                       % [m] attached chord length
+settings.Chord2 = 0.14;                                       % [m] free chord length
+settings.Height = 0.11;                                       % [m] fin heigth
+settings.deltaXLE = 0.13;                                     % [m] start of Chord 2 measured from start of Chord 1
+settings.Npanel = 3;                                          % [m] number of fins
+settings.Ler = [0 0];                                         % [deg] Leading edge radius at each span station
+settings.d = 0.012;                                           % [m] distance between end of root chord and end of center body
+settings.zupRaw = 0.00175;                                    % [m] fin semi-thickness
+settings.LmaxuRaw = 0.00175;                                  % [m] Fraction of chord from leading edge to max thickness
+
+%%% protub data
+settings.xprot = 1.517;                                       % [m] axial position from nosecone base
+settings.nloc = 3;                                            % [-] number of brakes
+settings.lprot = 0.008;                                       % [m] brakes thickness
+settings.wprot = 0.1002754821;                                % [m] brakes width (normal)
+
+%% AIRBRAKES CONTROL SETTINGS
+settings.tControl = settings.tb;                              % [s] time after which the airbrakes can be used
+settings.MachControl = 0.8;                                   % [-] Maximum Mach at which airbrakes can be used
+settings.vServo = 150*pi/180;                                 % [rad/s] Servo-motor angular velocity
+% First entry, must be always 0 !
+settings.hprot = linspace(0, 0.0363, 3);                      % [m] Block airbrakes opening coordinate
+
+%% MASS GEOMERTY DETAILS
+% x-axis: along the fuselage
+% y-axis: right wing
+% z-axis: downward
+
+%%% Inertias for fuselage only (no engine)
+settings.rocketIxx = 0.06535397;                              % [kg*m^2] Inertia to x-axis
+settings.rocketIyy = 12.07664659;                             % [kg*m^2] Inertia to y-axis
+settings.rocketIzz = 12.07701314;                             % [kg*m^2] Inertia to z-axis
+
+%%% Inertias for engine only
+engineIxx = HREmotors(iMotor).Ixx;                            % [kg*m^2] Inertia to x-axis
+engineIyy = HREmotors(iMotor).Iyy;                            % [kg*m^2] Inertia to y-axis
+engineIzz = HREmotors(iMotor).Izz;                            % [kg*m^2] Inertia to z-axis
+
+%%% Total inertias
+settings.Ixx = settings.rocketIxx + engineIxx;
+
+settings.Iyy = settings.rocketIyy + ((settings.xcgNoMot - settings.xcg).^2) .* settings.mNoMot + ...
+    (engineIyy + ((settings.motor.xcg + settings.LcenterRocket - settings.xcg).^2 ) .* (settings.motor.expM + settings.motor.mc));
+
+settings.Izz = settings.rocketIzz + ((settings.xcgNoMot - settings.xcg).^2) .* settings.mNoMot + ...
+    (engineIzz + ((settings.motor.xcg + settings.LcenterRocket - settings.xcg).^2) .* (settings.motor.expM + settings.motor.mc));
+
+settings.I = [settings.Ixx; settings.Iyy; settings.Izz];
+
+%%% Inertias derivatives (avoiding calculations in ascent.m)
+settings.Idot = diff(settings.I')'./diff(settings.motor.expTime);
+settings.Idot(:, end+1) = settings.Idot(:, end);
+
+%% APPLY CORRECTIONS
+
+if exist("correction.mat","file")
+    corr = load("correction.mat");
+    settings.mTotalTime = settings.mTotalTime + corr.deltaM;
+    settings.xcg = settings.xcg + reshapeVector(settings.motor.expTime, corr.time, corr.deltaXcg);
+    settings.ms = settings.ms + corr.deltaM;
+end
+
+%% AERODYNAMICS DETAILS
+% These coefficients are obtained using MISSILE DATCOM, check simulator/README.md for further details.
+
+% DATA_PATH = '../data/';
+
+filename_coeffs = strcat(engineName, '.mat');
+
+if not(isfield(settings,'autoMatProt_flag')) || settings.autoMatProt_flag == false
+    % Coefficients
+    Coeffs = load(filename_coeffs);
+    settings.Coeffs = Coeffs.CoeffsTot;
+
+    % Geometry
+    settings.Geometry = Coeffs.Geometry;
+    settings.Geometry.xcg = Coeffs.Geometry.xcg';
+
+    % State
+    settings.State.Alphas = Coeffs.State.Alphas';
+    settings.State.Betas = Coeffs.State.Betas';
+    settings.State.Altitudes = Coeffs.State.Altitudes';
+    settings.State.Machs = Coeffs.State.Machs';
+    settings.State.Controls = Coeffs.State.hprot';
+    settings.State.xcgTime = Coeffs.State.xcgTime';
+    clear('Coeffs');
+
+    CoeffHighAOA = load('HighAOA.mat'); 
+    settings.highAOA.Coeffs = CoeffHighAOA.CoeffsTot; 
+    settings.highAOA.State.Alphas = CoeffHighAOA.State.Alphas';
+    settings.highAOA.State.Betas = CoeffHighAOA.State.Betas'; 
+    settings.highAOA.State.Altitudes = CoeffHighAOA.State.Altitudes'; 
+    settings.highAOA.State.Machs = CoeffHighAOA.State.Machs';
+    settings.highAOA.State.Controls = CoeffHighAOA.State.hprot';
+    settings.highAOA.State.xcgTime = CoeffHighAOA.State.xcgTime';
+    clear('CoeffHighAOA');
+
+    settings.highAOA.tb = settings.tb; 
+    settings.highAOA.timeEngineCut = settings.timeEngineCut; 
+end
+
+%% DESCENT STAGES
+settings.NdescentStages = 2;                                  % [-] Number of descent stages
+% Followimg are for stage 2, 3, ..
+settings.msStages = 4.412;                                    % [kg] Mass of stage 2, 3, .. without chutes
+
+%% PARACHUTES DETAILS
+
+settings.Npara = [2 1];
+
+% parachute 1
+settings.para(1, 1).name = 'DROGUE chute';
+settings.para(1, 1).S = 1.2219;    % [m^2]   Surface
+settings.para(1, 1).mass = 0.15;   % [kg]   Parachute Mass
+settings.para(1, 1).CD = 0.96;     % [/] Parachute Drag Coefficient
+settings.para(1, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 1).delay = 1;     % [s] drogue opening delay
+settings.para(1, 1).z_cut = 350;   % [m] Final altitude of the parachute
+
+% parachute 2
+settings.para(2, 1).name = 'MAIN chute';
+settings.para(2, 1).S = 7.34;      % [m^2]   Surface
+settings.para(2, 1).mass = 1.05;    % [kg]   Parachute Mass
+settings.para(2, 1).CD = 1.75;     % [/] Parachute Drag Coefficient
+settings.para(2, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(2, 1).z_cut = 0;     % [m] Final altitude of the parachute
+
+% parachute 1
+settings.para(1, 1).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 1).L = 1.5;       % [m] Shock Chord Length
+settings.para(1, 1).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 1).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 1).Vexit = 5;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 1).CX = 1.2;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(2, 1).L = 6;         % [m] Shock Chord Length
+settings.para(2, 1).K = 3000;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(2, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(2, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(2, 1).nf = 8.7;      % [/] Adimensional Opening Time
+
+%%% PAYLOAD CHUTES
+% parachute 1
+settings.para(1, 2).name = "Payload DROGUE";
+settings.para(1, 2).S = 0.11;      % [m^2]   Surface
+settings.para(1, 2).mass = 0.15;   % [kg]   Parachute Mass
+settings.para(1, 2).CD = 1.2;      % [/] Parachute Drag Coefficient
+settings.para(1, 2).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 2).delay = 1;     % [s] drogue opening delay
+settings.para(1, 2).z_cut = 300;   % [m] Final altitude of the parachute
+settings.para(1, 2).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 2).L = 1.5;       % [m] Shock Chord Length
+settings.para(1, 2).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 2).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 2).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 2).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 2).Vexit = 10;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 2).name = "Payload AIRFOIL";
+settings.para(2, 2).mass = 0.50;    % [kg]   Parachute Mass
+settings.para(2, 2).z_cut = 0;      % [m] Final altitude of the parachute
+
+
+%% INTEGRATION OPTIONS
+settings.ode.finalTime =  2000;    % [s] Final integration time
+
+% create an option structure for the integrations:
+
+% - AbsTol is the threshold below which the value of the solution becomes unimportant
+% - RelTol is the tolerance betweeen two consecutive values
+% - Events is the event function that defines when the integration must be
+% - stopped (it has to be created)
+% - InitialStep is the highest value tried by the solver
+
+settings.ode.optionsasc1 = odeset('Events', @eventApogee, 'InitialStep', 1);          %ODE options for ascend
+settings.ode.optionsasc2 = odeset('InitialStep', 1);                                  %ODE options for due to the opening delay of the parachute
+settings.ode.optionsascMultipleAB = odeset('Events', @eventAB,  'InitialStep', 1);    %ODE options for opening of the airbrakes
+settings.ode.optionspara = odeset('Events', @eventParaCut);                           %ODE options for the parachutes
+settings.ode.optionsdesc = odeset('Events', @eventLanding);                           %ODE options for ballistic descent
+settings.ode.optionspad = odeset('Events', @eventPad);                                %ODE options for the launchpad phase
+
+% Settings for descent 6dof simulation
+settings.ode.optionsDrogue6DOF = odeset('Events', @eventParaCut, 'AbsTol', 1e-6,'RelTol', 1e-6);         %ODE options for due to cutting of the drogue chute
+settings.ode.optionsMainExt6DOF = odeset('Events', @eventMainExit, 'AbsTol', 1e-6,'RelTol', 1e-6);       %ODE options for due to the extraction of the main chute
+settings.ode.optionsMain6DOF = odeset('Events', @eventLanding, 'AbsTol', 1e-6,'RelTol', 1e-6);           %ODE options to terminate descent phase
+
+%% STOCHASTIC DETAILS
+%%% launch probability details
+settings.stoch.prob.xLim = 2e3;                    % Max ovest displacement [m]
+settings.stoch.prob.VLim = 50;                     % Max drogue velocity [Pa]
+settings.stoch.prob.XCPLim = 1.5;                  % Min XCP
+
+%%% Safe Ellipse - roccaraso
+settings.prob.SafeEllipse(1).name = 'prova';
+settings.prob.SafeEllipse(1).a = 1100;
+settings.prob.SafeEllipse(1).b = 2800;
+settings.prob.SafeEllipse(1).x0  = 0;
+settings.prob.SafeEllipse(1).y0 = -300;
+settings.prob.SafeEllipse(1).alpha = 10;
+
+
+%% LANDING POINTS
+% satellite maps of the landing zone
+% delta limit on the coordinates (for the landing map)
+settings.limLat = 0.04;
+settings.limLon = 0.025;
diff --git a/missions/2024_Lyra_Portugal_October/HRE_ARM_7_S_HE-T03T03.mat b/missions/2024_Lyra_Portugal_October/HRE_ARM_7_S_HE-T03T03.mat
new file mode 100644
index 0000000000000000000000000000000000000000..ee1ef38bca87efad2801bbaa461516ef62616bb9
--- /dev/null
+++ b/missions/2024_Lyra_Portugal_October/HRE_ARM_7_S_HE-T03T03.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:dcb450574d33e83bb89bce75cf5dca2490f9338fb3a4dcc5ededa0113464a2ac
+size 60380731
diff --git a/missions/2024_Lyra_Portugal_October/HighAOA_HRE_ARM_7_S_HE-T03T03.mat b/missions/2024_Lyra_Portugal_October/HighAOA_HRE_ARM_7_S_HE-T03T03.mat
new file mode 100644
index 0000000000000000000000000000000000000000..812573bedd1d966c1fa3ab7aff43897e52e63af9
--- /dev/null
+++ b/missions/2024_Lyra_Portugal_October/HighAOA_HRE_ARM_7_S_HE-T03T03.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7ab760661aa3085525d8f4c652f9a833a29977b1405c361f0a05c1891d42cfb5
+size 5777728
diff --git a/missions/2024_Lyra_Portugal_October/rocketStress.mat b/missions/2024_Lyra_Portugal_October/rocketStress.mat
new file mode 100644
index 0000000000000000000000000000000000000000..7f20d2d5812dd5cb1a7182b281f7ebf0e977e824
--- /dev/null
+++ b/missions/2024_Lyra_Portugal_October/rocketStress.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a2af604b52fb12df3ffc96825ccb04e4ec3fa9bce8a0f70c3aa2cb2ae4100177
+size 1138617
diff --git a/missions/2024_Lyra_Portugal_October/simulationsData.m b/missions/2024_Lyra_Portugal_October/simulationsData.m
new file mode 100644
index 0000000000000000000000000000000000000000..ca6c19ccd77a7ed5df159785635517d18510e370
--- /dev/null
+++ b/missions/2024_Lyra_Portugal_October/simulationsData.m
@@ -0,0 +1,371 @@
+%{
+
+ROCKET DATA FILE
+
+Project name:                         Lyra
+Year of launch:                       2024
+Location of Launch:      Pont de Sor (POR)
+Date of launch:                        TBD
+
+Last update:                    01/10/2023
+
+Copyright © 2021, Skyward Experimental Rocketry, AFD department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+% launchpad - Pont de Sor
+settings.z0 = 160;                                            % [m] Launchpad Altitude
+settings.lpin1 = 0.5603;                                      % [m] Distance from the upper pin to the upper tank cap
+settings.lpin2 = 0.2055;                                      % [m] Distance from the lower pin to the lower tank cap  
+settings.lrampa = 12;                                         % [m] Total launchpad length 
+settings.lat0 = 39.388727;                                    % Launchpad latitude
+settings.lon0 = -8.287842;                                    % Launchpad longitude
+settings.g0 = gravitywgs84(settings.z0, settings.lat0);       % Gravity costant at launch latitude and altitude
+
+settings.Local = [];                                          % vector [altLocal, TLocal, PLocal, rhoLocal], up to 4 components
+
+%% ENGINE DETAILS
+% load motors data
+settings.HREmot = true;                                       % non runna con false perche mancano le inerzie totali nel simulationsData
+
+filename = '../HREmotors.mat';
+load(filename);
+
+
+ engineName = 'HRE_ARM_7_S_HE-T03T03';
+
+
+iMotor = strcmp({HREmotors.MotorName}, engineName);
+settings.motor.expTime = HREmotors(iMotor).t;                 % [s] Engine time vector
+settings.motor.expThrust = HREmotors(iMotor).T;               % [N] Engine thrust vector
+settings.motor.expM = HREmotors(iMotor).m;                    % [Kg] Propellant Mass (in time)
+settings.motor.L = HREmotors(iMotor).L;                       % [m] Engine length
+settings.motor.Ltank = HREmotors(iMotor).Ltank;               % [m] Tank length
+settings.motor.xcg = HREmotors(iMotor).xcg;                   % [m] Engine xcg from tank tip
+settings.motor.mc = HREmotors(iMotor).mc;                     % [kg] Engine Structural Mass
+settings.motor.Pe = HREmotors(iMotor).Pe;                     % [Pa] Eflux pressure
+settings.motor.Ae = HREmotors(iMotor).Ae;                     % [Pa] Eflux Area
+settings.motor.mFuel = HREmotors(iMotor).mFu;                 % [kg] Fuel (grain) initial mass
+settings.motor.mOx = HREmotors(iMotor).mOx;                   % [kg] Oxidizer initial mass
+settings.mp = HREmotors(iMotor).mp;                           % [kg] Initial propellant Mass
+settings.tb = HREmotors(iMotor).t(end);                       % [s]  Burning time
+
+% Engine fuselage data
+settings.motor.mFuselage = HREmotors(iMotor).mFus;            % [kg] Fuselage of the engine only
+settings.motor.xcgFuselage = (settings.motor.L - ...
+settings.motor.Ltank)/2 + settings.motor.Ltank;               % [m]  xcg of the engine fuselage only from tank tip
+
+% Engine cut
+settings.timeEngineCut = inf;                                 % [s] Moment in time in wich the engine will be cut off
+settings.tIGN = 0.3;                                          % [s] ignition transient
+settings.tCO = 0.3;                                           % [s] cut-off transient
+
+%% Launchpad synthesis
+
+% settings.lpin = settings.lpin1 + settings.motor.Ltank ... 
+%                    + settings.lpin2;                         % [m] Distance between the two launchpad pins
+
+settings.lpin = 2.4 + 0.2;                                    % [m] distance of the upper pin from the rail base (upper pin-boat + boat-rail base)
+settings.lrampa = settings.lrampa - settings.lpin;            % [m] LaunchPad route (total available route)
+
+%% ROCKET DETAILS
+%%% Boat-tail
+% The boat-tail is simulated as a trunc of cone with the major diameter
+% equal to the rocket one (settings.C)
+settings.boatL = 0.114;                                        % [m] Boat-tail length
+settings.boatD = 0.125;                                       % [m] Boat-tail base diameter
+settings.boatType = 1;                                        % [-] Boat-tail shape. 0: Cone, 1: Tangent Ogive
+
+%%% Rocket ogive and central body
+settings.C = 0.150;                                           % [m] Caliber (Fuselage Diameter)
+settings.S = pi*settings.C^2/4;                               % [m^2] Cross-sectional Surface
+settings.Lnose = 0.28;                                        % [m] Nosecone Length
+settings.OgType = 'MHAACK';                                   % [-] Nosecone shape
+settings.NosePower = 3/4;                                     % [-] Nosecone power type parameter
+settings.pMod = 1.180236e+00;                                 % [-] P coefficient for modified nosecone shapes
+settings.cMod = 1.143441e-01;                                 % [-] C coefficient for modified nosecone shapes
+
+settings.LcenterRocket = 1.517;                               % [m] Central body length
+settings.mNoMot = 16.952;                                     % [kg] Mass of everything above engine
+settings.xcgNoMot = 1.1700;                                   % [m] XCG of everything above engine (from nosecone base)
+
+% SYNTHESIS WITH ENGINE FUSELAGE MASS AND XCG SEPARATELY GIVEN
+settings.addMfusMot = false;
+
+if settings.addMfusMot
+    settings.xcgNoMot = (settings.mNoMot*settings.xcgNoMot + settings.motor.mFuselage*(settings.motor.xcgFuselage + settings.LcenterRocket))...
+                        /(settings.mNoMot + settings.motor.mFuselage);
+    settings.mNoMot = settings.mNoMot + settings.motor.mFuselage;
+end
+
+%% SYNTHESIS
+settings.Lcenter = settings.LcenterRocket + settings.motor.L; % [m] Total centerbody length
+settings.mnc = 0;                                             % [kg] Nosecone Mass
+settings.ms = settings.mNoMot + settings.motor.mc + ...
+settings.mnc;                                                 % [kg] Total structural Mass
+settings.m0 = settings.ms + settings.motor.expM(1);           % [kg] Total initial Mass
+settings.mTotalTime = settings.mNoMot + settings.motor.mc + ...
+settings.mnc + settings.motor.expM;                           % Total mass (in time)
+
+% Total XCG of the rocket [m] !! FROM NOSECONE BASE !!
+settings.xcg = ((settings.motor.expM + settings.motor.mc) .* (settings.LcenterRocket + settings.motor.xcg) + ...
+    settings.mNoMot * settings.xcgNoMot) ./ (settings.motor.expM + settings.motor.mc + settings.mNoMot);
+
+settings.Lcenter = settings.Lcenter - settings.boatL;         % correcting the rocket length in order to account the boat tail
+
+%% AERODYNAMICS
+%%% Pitot tube
+settings.PITOTCLI = 0;                                        % [m] Pitot initial conic section length
+settings.PITOTCLB = 0;                                        % [m] Pitot final conic section length
+settings.PITOTCDI = 0;                                        % [m] Pitot initial conic section diameter
+settings.PITOTCDB = 0;                                        % [m] Pitot final conic section diameter
+settings.PITOTL = 0;                                          % [m] Pitot tube length
+settings.PITOTD = 0;                                          % [m] Pitot tube diameter
+
+%%% finset
+settings.Chord1 = 0.27;                                       % [m] attached chord length
+settings.Chord2 = 0.11;                                       % [m] free chord length
+settings.Height = 0.13;                                       % [m] fin heigth
+settings.deltaXLE = 0.14;                                     % [m] start of Chord 2 measured from start of Chord 1
+settings.Npanel = 3;                                          % [m] number of fins
+settings.Ler = [0 0];                                         % [deg] Leading edge radius at each span station
+settings.d = -0.044;                                           % [m] distance between end of root chord and end of center body
+settings.zupRaw = 0.00175;                                    % [m] fin semi-thickness
+settings.LmaxuRaw = 0.00175;                                  % [m] Fraction of chord from leading edge to max thickness
+
+%%% protub data
+settings.xprot = 1.517;                                       % [m] axial position from nosecone base
+settings.nloc = 3;                                            % [-] number of brakes
+settings.lprot = 0.008;                                       % [m] brakes thickness
+settings.wprot = 0.1002754821;                                % [m] brakes width (normal)
+
+%% AIRBRAKES CONTROL SETTINGS
+settings.tControl = settings.tb;                              % [s] time after which the airbrakes can be used
+settings.MachControl = 0.8;                                   % [-] Maximum Mach at which airbrakes can be used
+settings.vServo = 150*pi/180;                                 % [rad/s] Servo-motor angular velocity
+% First entry, must be always 0 !
+settings.hprot = linspace(0, 0.0363, 3);                      % [m] Block airbrakes opening coordinate
+
+%% MASS GEOMERTY DETAILS
+% x-axis: along the fuselage
+% y-axis: right wing
+% z-axis: downward
+
+%%% Inertias for fuselage only (no engine)
+settings.rocketIxx = 0.06535397;                              % [kg*m^2] Inertia to x-axis
+settings.rocketIyy = 17.21019828;                             % [kg*m^2] Inertia to y-axis
+settings.rocketIzz = 17.21056483;                             % [kg*m^2] Inertia to z-axis
+
+%%% Inertias for engine only
+engineIxx = HREmotors(iMotor).Ixx;                            % [kg*m^2] Inertia to x-axis
+engineIyy = HREmotors(iMotor).Iyy;                            % [kg*m^2] Inertia to y-axis
+engineIzz = HREmotors(iMotor).Izz;                            % [kg*m^2] Inertia to z-axis
+
+%%% Total inertias
+settings.Ixx = settings.rocketIxx + engineIxx;
+
+settings.Iyy = settings.rocketIyy + ((settings.xcgNoMot - settings.xcg).^2).*settings.mNoMot + ...
+    (engineIyy + ( (settings.motor.xcg + settings.LcenterRocket - settings.xcg).^2 ) .*(settings.motor.expM + settings.motor.mc));
+
+settings.Izz = settings.rocketIzz + ((settings.xcgNoMot - settings.xcg).^2).*settings.mNoMot + ...
+    (engineIzz + ( (settings.motor.xcg + settings.LcenterRocket - settings.xcg).^2).*(settings.motor.expM + settings.motor.mc));
+
+settings.I = [settings.Ixx; settings.Iyy; settings.Izz];
+
+%%% Inertias derivatives (avoiding calculations in ascent.m)
+settings.Idot = diff(settings.I')'./diff(settings.motor.expTime);
+settings.Idot(:, end+1) = settings.Idot(:, end);
+
+%% APPLY CORRECTIONS
+
+if exist("correction.mat","file")
+    corr = load("correction.mat");
+    settings.mTotalTime = settings.mTotalTime + corr.deltaM;
+    settings.xcg = settings.xcg + reshapeVector(settings.motor.expTime, corr.time, corr.deltaXcg);
+    settings.ms = settings.ms + corr.deltaM;
+end
+
+%% AERODYNAMICS DETAILS
+% These coefficients are obtained using MISSILE DATCOM, check simulator/README.md for further details.
+
+% DATA_PATH = '../data/';
+
+filename_coeffs = strcat(engineName, '.mat');
+
+if not(isfield(settings,'autoMatProt_flag')) || settings.autoMatProt_flag == false
+    % Coefficients
+    Coeffs = load(filename_coeffs);
+    settings.Coeffs = Coeffs.CoeffsTot;
+
+    % Geometry
+    settings.Geometry = Coeffs.Geometry;
+    settings.Geometry.xcg = Coeffs.Geometry.xcg';
+
+    % State
+    settings.State.Alphas = Coeffs.State.Alphas';
+    settings.State.Betas = Coeffs.State.Betas';
+    settings.State.Altitudes = Coeffs.State.Altitudes';
+    settings.State.Machs = Coeffs.State.Machs';
+    settings.State.Controls = Coeffs.State.hprot';
+    settings.State.xcgTime = Coeffs.State.xcgTime';
+    clear('Coeffs');
+    
+    CoeffHighAOA = load(strcat('HighAOA_', filename_coeffs)); 
+    settings.highAOA.Coeffs = CoeffHighAOA.CoeffsTot; 
+    settings.highAOA.State.Alphas = CoeffHighAOA.State.Alphas';
+    settings.highAOA.State.Betas = CoeffHighAOA.State.Betas'; 
+    settings.highAOA.State.Altitudes = CoeffHighAOA.State.Altitudes'; 
+    settings.highAOA.State.Machs = CoeffHighAOA.State.Machs';
+    settings.highAOA.State.Controls = CoeffHighAOA.State.hprot';
+    settings.highAOA.State.xcgTime = CoeffHighAOA.State.xcgTime';
+    clear('CoeffHighAOA');
+
+    settings.highAOA.tb = settings.tb; 
+    settings.highAOA.timeEngineCut = settings.timeEngineCut; 
+
+end
+
+%% DESCENT STAGES
+settings.NdescentStages = 2;                                  % [-] Number of descent stages
+% Followimg are for stage 2, 3, ..
+settings.msStages = 4.412;                                    % [kg] Mass of stage 2, 3, .. without chutes
+
+%% PARACHUTES DETAILS
+
+settings.Npara = [2 1];
+
+% parachute 1
+settings.para(1, 1).name = 'DROGUE chute';
+settings.para(1, 1).S = 0.6;    % [m^2]   Surface
+settings.para(1, 1).mass = 0.15;   % [kg]   Parachute Mass
+settings.para(1, 1).CD = 0.75;     % [/] Parachute Drag Coefficient
+settings.para(1, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 1).delay = 1;     % [s] drogue opening delay
+settings.para(1, 1).z_cut = 350;   % [m] Final altitude of the parachute
+
+% parachute 2
+settings.para(2, 1).name = 'MAIN chute';
+settings.para(2, 1).S = 7.34;      % [m^2]   Surface
+settings.para(2, 1).mass = 1.05;    % [kg]   Parachute Mass
+settings.para(2, 1).CD = 1.75;     % [/] Parachute Drag Coefficient
+settings.para(2, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(2, 1).z_cut = 0;     % [m] Final altitude of the parachute
+
+% parachute 1
+settings.para(1, 1).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 1).L = 1.5;       % [m] Shock Chord Length
+settings.para(1, 1).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 1).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 1).Vexit = 5;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 1).CX = 1.2;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(2, 1).L = 6;         % [m] Shock Chord Length
+settings.para(2, 1).K = 3000;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(2, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(2, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(2, 1).nf = 8.7;      % [/] Adimensional Opening Time
+
+%%% PAYLOAD CHUTES
+% parachute 1
+settings.para(1, 2).name = "Payload DROGUE";
+settings.para(1, 2).S = 0.11;      % [m^2]   Surface
+settings.para(1, 2).mass = 0.15;   % [kg]   Parachute Mass
+settings.para(1, 2).CD = 1.2;      % [/] Parachute Drag Coefficient
+settings.para(1, 2).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 2).delay = 1;     % [s] drogue opening delay
+settings.para(1, 2).z_cut = 450;   % [m] Final altitude of the parachute
+settings.para(1, 2).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 2).L = 1.5;       % [m] Shock Chord Length
+settings.para(1, 2).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 2).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 2).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 2).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 2).Vexit = 10;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 2).name = "Payload AIRFOIL";
+settings.para(2, 2).mass = 0.50;    % [kg]   Parachute Mass
+settings.para(2, 2).z_cut = 0;      % [m] Final altitude of the parachute
+
+
+%% PARAFOIL DETAILS
+% Geometry
+settings.payload.paraPrev = [2, 1];         % Reference to the previous parachute in the multistage descent 
+                                            % (in this case, the payload drogue: 2nd stage, 1st parachute)
+
+settings.payload.mass = 4.2;                 % [kg]  mass 
+settings.payload.b    = 2.55/2;              % [m]   semiwingspan  - vela nuova: 2.55/2; - vela vecchia: 2.06/2;
+settings.payload.c    = 0.8;                 % [m]   mean aero chord
+settings.payload.S    = 2.04;                % [m^2] payload surface - vela nuova 2.04; - vela vecchia: 1.64;  
+settings.payload.inertia = [0.42, 0,   0.03;
+                            0,    0.4,    0; 
+                            0.03, 0, 0.053]; % [kg m^2] [3x3] inertia matrix payload+payload 
+
+settings.payload.inverseInertia = inv(settings.payload.inertia);
+
+% Aerodynamic constants
+settings.payload.CD0       =  0.25; 
+settings.payload.CDAlpha2  =  0.12;
+settings.payload.CL0       =  0.091;
+settings.payload.CLAlpha   =  0.9;
+settings.payload.Cm0       =  0.35; 
+settings.payload.CmAlpha   = -0.72;
+settings.payload.Cmq       = -1.49;
+settings.payload.CLDeltaA  = -0.0035;
+settings.payload.Cnr       = -0.27;
+settings.payload.CnDeltaA  =  0.0115;
+settings.payload.deltaSMax =  0.1;
+settings.payload.CDDeltaA  = 0.01;
+settings.payload.Clp       = -0.84;
+settings.payload.ClPhi     = -0.1;
+settings.payload.ClDeltaA  = 0.01;
+
+%% INTEGRATION OPTIONS
+settings.ode.finalTime =  2000;    % [s] Final integration time
+
+% create an option structure for the integrations:
+
+% - AbsTol is the threshold below which the value of the solution becomes unimportant
+% - RelTol is the tolerance betweeen two consecutive values
+% - Events is the event function that defines when the integration must be
+% - stopped (it has to be created)
+% - InitialStep is the highest value tried by the solver
+
+settings.ode.optionsasc1 = odeset('Events', @eventApogee, 'InitialStep', 1);          %ODE options for ascend
+settings.ode.optionsasc2 = odeset('InitialStep', 1);                                  %ODE options for due to the opening delay of the parachute
+settings.ode.optionsascMultipleAB = odeset('Events', @eventAB,  'InitialStep', 1);    %ODE options for opening of the airbrakes
+settings.ode.optionspara = odeset('Events', @eventParaCut);                           %ODE options for the parachutes
+settings.ode.optionsdesc = odeset('Events', @eventLanding,...
+                                  'RelTol', 1e-3, 'AbsTol', 1e-3);                    %ODE options for ballistic descent
+settings.ode.optionspad = odeset('Events', @eventPad);                                %ODE options for the launchpad phase
+
+% Settings for descent 6dof simulation
+settings.ode.optionsDrogue6DOF = odeset('Events', @eventParaCut, 'AbsTol', 1e-6,'RelTol', 1e-6);         %ODE options for due to cutting of the drogue chute
+settings.ode.optionsMainExt6DOF = odeset('Events', @eventMainExit, 'AbsTol', 1e-6,'RelTol', 1e-6);       %ODE options for due to the extraction of the main chute
+settings.ode.optionsMain6DOF = odeset('Events', @eventLanding, 'AbsTol', 1e-6,'RelTol', 1e-6);           %ODE options to terminate descent phase
+
+%% STOCHASTIC DETAILS
+%%% launch probability details
+settings.stoch.prob.xLim = 2e3;                    % Max ovest displacement [m]
+settings.stoch.prob.VLim = 50;                     % Max drogue velocity [Pa]
+settings.stoch.prob.XCPLim = 1.5;                  % Min XCP
+
+%%% Safe Ellipse - roccaraso
+settings.prob.SafeEllipse(1).name = 'prova';
+settings.prob.SafeEllipse(1).a = 1100;
+settings.prob.SafeEllipse(1).b = 2800;
+settings.prob.SafeEllipse(1).x0  = 0;
+settings.prob.SafeEllipse(1).y0 = -300;
+settings.prob.SafeEllipse(1).alpha = 10;
+
+%% LANDING POINTS
+% satellite maps of the landing zone
+% delta limit on the coordinates (for the landing map)
+settings.limLat = 0.04;
+settings.limLon = 0.025;
diff --git a/missions/2024_Lyra_Roccaraso_September/simulationsData.m b/missions/2024_Lyra_Roccaraso_September/simulationsData.m
new file mode 100644
index 0000000000000000000000000000000000000000..f3ee141cbbbcfaf8f0d285c4f1d3141bcb7bca8e
--- /dev/null
+++ b/missions/2024_Lyra_Roccaraso_September/simulationsData.m
@@ -0,0 +1,335 @@
+%{
+
+ROCKET DATA FILE
+
+Project name:                       Crater
+Year of launch:                       2024
+Location of Launch:       Roccasarso (ITA)
+Date of launch:                        TBD
+
+Last update:                    01/10/2024
+
+Copyright © 2021, Skyward Experimental Rocketry, AFD department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+% launchpad - Roccaraso
+settings.z0 = 1414;                                           % [m] Launchpad Altitude
+settings.lpin1 = 000;                                      % [m] Distance from the upper pin to the upper tank cap
+settings.lpin2 = 000;                                      % [m] Distance from the lower pin to the lower tank cap  
+settings.lrampa = 000;                                         % [m] Total launchpad length 
+settings.lat0 = 41.8084579;                                   % Launchpad latitude
+settings.lon0 = 14.0546408;                                   % Launchpad longitude
+settings.g0 = gravitywgs84(settings.z0, settings.lat0);       % Gravity costant at launch latitude and altitude
+
+settings.Local = [];            %% ENGINE DETAILS
+
+%% load motors data
+settings.HREmot = true;                                       % non runna con false perche mancano le inerzie totali nel simulationsData
+
+filename = '../HREmotors.mat';
+load(filename);
+
+% engineName = 'HRE_FURIA-Rv2-T04T03'; 
+engineName = 'TBD';
+
+iMotor = strcmp({HREmotors.MotorName}, engineName);
+settings.motor.expTime = HREmotors(iMotor).t;                 % [s] Engine time vector
+settings.motor.expThrust = HREmotors(iMotor).T;               % [N] Engine thrust vector
+settings.motor.expM = HREmotors(iMotor).m;                    % [Kg] Propellant Mass (in time)
+settings.motor.L = HREmotors(iMotor).L;                       % [m] Engine length
+settings.motor.Ltank = HREmotors(iMotor).Ltank;               % [m] Tank length
+settings.motor.xcg = HREmotors(iMotor).xcg;                   % [m] Engine xcg from tank tip
+settings.motor.mc = HREmotors(iMotor).mc;                     % [kg] Engine Structural Mass
+settings.motor.Pe = HREmotors(iMotor).Pe;                     % [Pa] Eflux pressure
+settings.motor.Ae = HREmotors(iMotor).Ae;                     % [Pa] Eflux Area
+settings.motor.mFuel = HREmotors(iMotor).mFu;                 % [kg] Fuel (grain) initial mass
+settings.motor.mOx = HREmotors(iMotor).mOx;                   % [kg] Oxidizer initial mass
+settings.mp = HREmotors(iMotor).mp;                           % [kg] Initial propellant Mass
+settings.tb = HREmotors(iMotor).t(end) ;                      % [s]  Burning time
+
+% Engine fuselage data
+settings.motor.mFuselage = HREmotors(iMotor).mFus;            % [kg] Fuselage of the engine only
+settings.motor.xcgFuselage = (settings.motor.L - ...
+settings.motor.Ltank)/2 + settings.motor.Ltank;               % [m]  xcg of the engine fuselage only from tank tip
+
+% Engine cut
+settings.timeEngineCut = inf;                                 % [s] Moment in time in wich the engine will be cut off
+settings.tIGN = 000; %0.4;                                          % [s] ignition transient
+settings.tCO = 000; %0.3;                                           % [s] cut-off transient
+
+%% Launchpad synthesis
+
+%settings.lpin = settings.lpin1 + settings.motor.Ltank ... 
+%                    + settings.lpin2;                        % [m] Distance between the two launchpad pins
+
+settings.lpin = 000;%1.926 + 0.2;                             % [m] distance of the upper pin from the rail base (upper pin-boat + boat-rail base)
+settings.lrampa = settings.lrampa - settings.lpin;            % [m] LaunchPad route (total available route)
+
+%% ROCKET DETAILS
+%%% Boat-tail
+% The boat-tail is simulated as a trunc of cone with the major diameter
+% equal to the rocket one (settings.C)
+settings.boatL = 000;                                        % [m] Boat-tail length
+settings.boatD = 000;                                       % [m] Boat-tail base diameter
+
+%%% Rocket ogive and central body
+settings.C = 000;                                           % [m] Caliber (Fuselage Diameter)
+settings.S = pi*settings.C^2/4;                               % [m^2] Cross-sectional Surface
+settings.Lnose = 000;                                        % [m] Nosecone Length
+settings.OgType = 'MHAACK';                                   % [-] Nosecone shape
+settings.NosePower = 000; %3/4;                                     % [-] Nosecone power type parameter
+settings.pMod = 000; %1.250152e+00;                                 % [-] P coefficient for modified nosecone shapes
+settings.cMod = 000; %1.799127e-01;                                 % [-] C coefficient for modified nosecone shapes
+
+settings.LcenterRocket = 000; %1.517;                               % [m] Central body length
+settings.mNoMot = 000; %16.423;                                       % [kg] Mass of everything above engine
+settings.xcgNoMot = 000; %1.149;                                    % [m] XCG of everything above engine (from nosecone base)
+
+%SYNTHESIS WITH ENGINE FUSELAGE MASS AND XCG SEPARATELY GIVEN
+settings.addMfusMot = false;
+
+if settings.addMfusMot
+    settings.xcgNoMot = (settings.mNoMot*settings.xcgNoMot + settings.motor.mFuselage*(settings.motor.xcgFuselage + settings.LcenterRocket))...
+                        /(settings.mNoMot + settings.motor.mFuselage);
+    settings.mNoMot = settings.mNoMot + settings.motor.mFuselage;
+end
+
+%% SYNTHESIS
+settings.Lcenter = settings.LcenterRocket + settings.motor.L; % [m] Total centerbody length
+settings.mnc = 0;                                             % [kg] Nosecone Mass
+settings.ms = settings.mNoMot + settings.motor.mc + ...
+settings.mnc;                                                 % [kg] Total structural Mass
+settings.m0 = settings.ms + settings.motor.expM(1);           % [kg] Total initial Mass
+settings.mTotalTime = settings.mNoMot + settings.motor.mc + ...
+settings.mnc + settings.motor.expM;                           % Total mass (in time)
+
+% Total XCG of the rocket [m] !! FROM NOSECONE BASE !!
+settings.xcg = ( (settings.motor.expM + settings.motor.mc).*(settings.LcenterRocket + settings.motor.xcg) + ...
+    settings.mNoMot*settings.xcgNoMot)./(settings.motor.expM + settings.motor.mc + settings.mNoMot);
+
+settings.Lcenter = settings.Lcenter - settings.boatL;         % correcting the rocket length in order to account the boat tail
+
+%% AERODYNAMICS
+%%% Pitot tube
+settings.PITOTCLI = 0;                                        % [m] Pitot initial conic section length
+settings.PITOTCLB = 0;                                        % [m] Pitot final conic section length
+settings.PITOTCDI = 0;                                        % [m] Pitot initial conic section diameter
+settings.PITOTCDB = 0;                                        % [m] Pitot final conic section diameter
+settings.PITOTL = 0;                                          % [m] Pitot tube length
+settings.PITOTD = 0;                                          % [m] Pitot tube diameter
+
+%%% finset 
+settings.Chord1 = 0.30;                                       % [m] attached chord length
+settings.Chord2 = 0.14;                                       % [m] free chord length
+settings.Height = 0.11;                                       % [m] fin heigth
+settings.deltaXLE = 0.13;                                     % [m] start of Chord 2 measured from start of Chord 1
+settings.Npanel = 3;                                          % [m] number of fins
+settings.Ler = [0 0];                                         % [deg] Leading edge radius at each span station
+settings.d = 0.012;                                           % [m] distance between end of root chord and end of center body
+settings.zupRaw = 0.00175;                                    % [m] fin semi-thickness
+settings.LmaxuRaw = 0.00175;                                  % [m] Fraction of chord from leading edge to max thickness
+
+%%% protub data
+settings.xprot = 1.517;                                       % [m] axial position from nosecone base
+settings.nloc = 3;                                            % [-] number of brakes
+settings.lprot = 0.008;                                       % [m] brakes thickness
+settings.wprot = 0.1002754821;                                % [m] brakes width (normal)
+
+%% AIRBRAKES CONTROL SETTINGS
+settings.tControl = settings.tb;                              % [s] time after which the airbrakes can be used
+settings.MachControl = 0.8;                                   % [-] Maximum Mach at which airbrakes can be used
+settings.vServo = 150*pi/180;                                 % [rad/s] Servo-motor angular velocity
+% First entry, must be always 0 !
+settings.hprot = linspace(0, 0.0363, 3);                      % [m] Block airbrakes opening coordinate
+
+%% MASS GEOMERTY DETAILS
+% x-axis: along the fuselage
+% y-axis: right wing
+% z-axis: downward
+
+%%% Inertias for fuselage only (no engine)
+settings.rocketIxx = 0.06535397;                              % [kg*m^2] Inertia to x-axis
+settings.rocketIyy = 12.07664659;                             % [kg*m^2] Inertia to y-axis
+settings.rocketIzz = 12.07701314;                             % [kg*m^2] Inertia to z-axis
+
+%%% Inertias for engine only
+engineIxx = HREmotors(iMotor).Ixx;                            % [kg*m^2] Inertia to x-axis
+engineIyy = HREmotors(iMotor).Iyy;                            % [kg*m^2] Inertia to y-axis
+engineIzz = HREmotors(iMotor).Izz;                            % [kg*m^2] Inertia to z-axis
+
+%%% Total inertias
+settings.Ixx = settings.rocketIxx + engineIxx;
+
+settings.Iyy = settings.rocketIyy + ((settings.xcgNoMot - settings.xcg).^2) .* settings.mNoMot + ...
+    (engineIyy + ((settings.motor.xcg + settings.LcenterRocket - settings.xcg).^2 ) .* (settings.motor.expM + settings.motor.mc));
+
+settings.Izz = settings.rocketIzz + ((settings.xcgNoMot - settings.xcg).^2) .* settings.mNoMot + ...
+    (engineIzz + ((settings.motor.xcg + settings.LcenterRocket - settings.xcg).^2) .* (settings.motor.expM + settings.motor.mc));
+
+settings.I = [settings.Ixx; settings.Iyy; settings.Izz];
+
+%%% Inertias derivatives (avoiding calculations in ascent.m)
+settings.Idot = diff(settings.I')'./diff(settings.motor.expTime);
+settings.Idot(:, end+1) = settings.Idot(:, end);
+
+%% APPLY CORRECTIONS
+
+if exist("correction.mat","file")
+    corr = load("correction.mat");
+    settings.mTotalTime = settings.mTotalTime + corr.deltaM;
+    settings.xcg = settings.xcg + reshapeVector(settings.motor.expTime, corr.time, corr.deltaXcg);
+    settings.ms = settings.ms + corr.deltaM;
+end
+
+%% AERODYNAMICS DETAILS
+% These coefficients are obtained using MISSILE DATCOM, check simulator/README.md for further details.
+
+% DATA_PATH = '../data/';
+
+filename_coeffs = strcat(engineName, '.mat');
+
+if not(isfield(settings,'autoMatProt_flag')) || settings.autoMatProt_flag == false
+    % Coefficients
+    Coeffs = load(filename_coeffs);
+    settings.Coeffs = Coeffs.CoeffsTot;
+
+    % Geometry
+    settings.Geometry = Coeffs.Geometry;
+    settings.Geometry.xcg = Coeffs.Geometry.xcg';
+
+    % State
+    settings.State.Alphas = Coeffs.State.Alphas';
+    settings.State.Betas = Coeffs.State.Betas';
+    settings.State.Altitudes = Coeffs.State.Altitudes';
+    settings.State.Machs = Coeffs.State.Machs';
+    settings.State.Controls = Coeffs.State.hprot';
+    settings.State.xcgTime = Coeffs.State.xcgTime';
+    clear('Coeffs');
+
+    CoeffHighAOA = load('HighAOA.mat'); 
+    settings.highAOA.Coeffs = CoeffHighAOA.CoeffsTot; 
+    settings.highAOA.State.Alphas = CoeffHighAOA.State.Alphas';
+    settings.highAOA.State.Betas = CoeffHighAOA.State.Betas'; 
+    settings.highAOA.State.Altitudes = CoeffHighAOA.State.Altitudes'; 
+    settings.highAOA.State.Machs = CoeffHighAOA.State.Machs';
+    settings.highAOA.State.Controls = CoeffHighAOA.State.hprot';
+    settings.highAOA.State.xcgTime = CoeffHighAOA.State.xcgTime';
+    clear('CoeffHighAOA');
+
+    settings.highAOA.tb = settings.tb; 
+    settings.highAOA.timeEngineCut = settings.timeEngineCut; 
+end
+
+%% DESCENT STAGES
+settings.NdescentStages = 2;                                  % [-] Number of descent stages
+% Followimg are for stage 2, 3, ..
+settings.msStages = 4.412;                                    % [kg] Mass of stage 2, 3, .. without chutes
+
+%% PARACHUTES DETAILS
+
+settings.Npara = [2 1];
+
+% parachute 1
+settings.para(1, 1).name = 'DROGUE chute';
+settings.para(1, 1).S = 1.2219;    % [m^2]   Surface
+settings.para(1, 1).mass = 0.15;   % [kg]   Parachute Mass
+settings.para(1, 1).CD = 0.96;     % [/] Parachute Drag Coefficient
+settings.para(1, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 1).delay = 1;     % [s] drogue opening delay
+settings.para(1, 1).z_cut = 350;   % [m] Final altitude of the parachute
+
+% parachute 2
+settings.para(2, 1).name = 'MAIN chute';
+settings.para(2, 1).S = 7.34;      % [m^2]   Surface
+settings.para(2, 1).mass = 1.05;    % [kg]   Parachute Mass
+settings.para(2, 1).CD = 1.75;     % [/] Parachute Drag Coefficient
+settings.para(2, 1).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(2, 1).z_cut = 0;     % [m] Final altitude of the parachute
+
+% parachute 1
+settings.para(1, 1).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 1).L = 1.5;       % [m] Shock Chord Length
+settings.para(1, 1).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 1).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 1).Vexit = 5;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 1).CX = 1.2;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(2, 1).L = 6;         % [m] Shock Chord Length
+settings.para(2, 1).K = 3000;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(2, 1).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(2, 1).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(2, 1).nf = 8.7;      % [/] Adimensional Opening Time
+
+%%% PAYLOAD CHUTES
+% parachute 1
+settings.para(1, 2).name = "Payload DROGUE";
+settings.para(1, 2).S = 0.11;      % [m^2]   Surface
+settings.para(1, 2).mass = 0.15;   % [kg]   Parachute Mass
+settings.para(1, 2).CD = 1.2;      % [/] Parachute Drag Coefficient
+settings.para(1, 2).CL = 0;        % [/] Parachute Lift Coefficient
+settings.para(1, 2).delay = 1;     % [s] drogue opening delay
+settings.para(1, 2).z_cut = 300;   % [m] Final altitude of the parachute
+settings.para(1, 2).CX = 1.4;      % [/] Parachute Longitudinal Drag Coefficient
+settings.para(1, 2).L = 1.5;       % [m] Shock Chord Length
+settings.para(1, 2).K = 7200;      % [N/m^2] Shock Chord Elastic Constant
+settings.para(1, 2).C = 0;         % [Ns/m] Shock Chord Dynamic Coefficient
+settings.para(1, 2).m = 1;         % [m^2/s] Coefficient of the surface vs. time opening model
+settings.para(1, 2).nf = 12;       % [/] Adimensional Opening Time
+settings.para(1, 2).Vexit = 10;     % [m/s] Expulsion Speed
+
+% parachute 2
+settings.para(2, 2).name = "Payload AIRFOIL";
+settings.para(2, 2).mass = 0.50;    % [kg]   Parachute Mass
+settings.para(2, 2).z_cut = 0;      % [m] Final altitude of the parachute
+
+
+%% INTEGRATION OPTIONS
+settings.ode.finalTime =  2000;    % [s] Final integration time
+
+% create an option structure for the integrations:
+
+% - AbsTol is the threshold below which the value of the solution becomes unimportant
+% - RelTol is the tolerance betweeen two consecutive values
+% - Events is the event function that defines when the integration must be
+% - stopped (it has to be created)
+% - InitialStep is the highest value tried by the solver
+
+settings.ode.optionsasc1 = odeset('Events', @eventApogee, 'InitialStep', 1);          %ODE options for ascend
+settings.ode.optionsasc2 = odeset('InitialStep', 1);                                  %ODE options for due to the opening delay of the parachute
+settings.ode.optionsascMultipleAB = odeset('Events', @eventAB,  'InitialStep', 1);    %ODE options for opening of the airbrakes
+settings.ode.optionspara = odeset('Events', @eventParaCut);                           %ODE options for the parachutes
+settings.ode.optionsdesc = odeset('Events', @eventLanding);                           %ODE options for ballistic descent
+settings.ode.optionspad = odeset('Events', @eventPad);                                %ODE options for the launchpad phase
+
+% Settings for descent 6dof simulation
+settings.ode.optionsDrogue6DOF = odeset('Events', @eventParaCut, 'AbsTol', 1e-6,'RelTol', 1e-6);         %ODE options for due to cutting of the drogue chute
+settings.ode.optionsMainExt6DOF = odeset('Events', @eventMainExit, 'AbsTol', 1e-6,'RelTol', 1e-6);       %ODE options for due to the extraction of the main chute
+settings.ode.optionsMain6DOF = odeset('Events', @eventLanding, 'AbsTol', 1e-6,'RelTol', 1e-6);           %ODE options to terminate descent phase
+
+%% STOCHASTIC DETAILS
+%%% launch probability details
+settings.stoch.prob.xLim = 2e3;                    % Max ovest displacement [m]
+settings.stoch.prob.VLim = 50;                     % Max drogue velocity [Pa]
+settings.stoch.prob.XCPLim = 1.5;                  % Min XCP
+
+%%% Safe Ellipse - roccaraso
+settings.prob.SafeEllipse(1).name = 'prova';
+settings.prob.SafeEllipse(1).a = 1100;
+settings.prob.SafeEllipse(1).b = 2800;
+settings.prob.SafeEllipse(1).x0  = 0;
+settings.prob.SafeEllipse(1).y0 = -300;
+settings.prob.SafeEllipse(1).alpha = 10;
+
+
+%% LANDING POINTS
+% satellite maps of the landing zone
+% delta limit on the coordinates (for the landing map)
+settings.limLat = 0.04;
+settings.limLon = 0.025;
diff --git a/missions/HREmotors.mat b/missions/HREmotors.mat
new file mode 100644
index 0000000000000000000000000000000000000000..1ffa2dbb3d70ca3e4350f7c90acc37ec7549e8ee
--- /dev/null
+++ b/missions/HREmotors.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fbb0e6fb46b31cadbc084461987aa8f45d9128c7e6b8d59483b88809e1495f3f
+size 6730820
diff --git a/missions/Motors.mat b/missions/Motors.mat
new file mode 100644
index 0000000000000000000000000000000000000000..4a1e5652675f11a71eb9459430174f255ba6d74b
--- /dev/null
+++ b/missions/Motors.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5d321100dbc7d5a66a7a456ffce103f26a2b2055706f70533f5212b400ef1003
+size 242474
diff --git a/missions/README.md b/missions/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..722b350741149ea19998590a3d590a43b5915abb
--- /dev/null
+++ b/missions/README.md
@@ -0,0 +1,74 @@
+# data
+This folder contains all the engines and the list of the rockets which are available within the repository.
+
+
+## Gemini_Portugal_October_2023 ![](https://img.shields.io/badge/-work%20in%20progress-green)
+**Project name:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;  &nbsp;Gemini <br />
+**Year of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 2023 <br />
+**Location of Launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Pont de Sor (POR) <br />
+**Date of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;  TBD <br />
+<br /><br />
+
+
+## Gemini_Roccaraso_September_2023 ![](https://img.shields.io/badge/-work%20in%20progress-green)
+**Project name:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;Gemini <br />
+**Year of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 2023 <br />
+**Location of Launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&nbsp; Roccaraso (ITA) <br />
+**Date of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;  TBD <br />
+<br /><br />
+
+
+## Pyxis_Portugal_October_2022
+**Project name:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Pyxis <br />
+**Year of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 2022 <br />
+**Location of Launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Pont de Sor (POR) <br />
+**Date of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;13/10/2022 <br />
+<br /><br />
+
+
+## Pyxis_Roccaraso_September_2022
+**Project name:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Pyxis <br />
+**Year of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 2022 <br />
+**Location of Launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Roccaraso (ITA) <br />
+**Date of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;  17/09/2022 <br />
+<br /><br />
+
+
+## Lynx_Portugal_October_2021
+**Project name:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Lynx <br />
+**Year of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 2021 <br />
+**Location of Launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Pont de Sor (POR) <br />
+**Date of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&nbsp;&nbsp; 13/10/2021 <br />
+<br /><br />
+
+
+
+## Lynx_Roccaraso_September_2021
+**Project name:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Lynx <br />
+**Year of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 2021 <br />
+**Location of Launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Roccaraso (ITA) <br />
+**Date of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 18/09/2021 <br />
+<br /><br />
+
+
+## HermesV1_Roccaraso_November_2019
+**Project name:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Lynx <br />
+**Year of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 2019 <br />
+**Location of Launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Roccaraso (ITA) <br />
+**Date of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 15/11/2019 <br /> <br /> <br /> <br />
+
+
+
+## R2A_Sardinia_July_2017
+**Project name:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; R2A <br />
+**Year of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 2017 <br />
+**Location of Launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Sardinia (ITA) <br />
+**Date of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;07/2017 <br /> <br /> <br /> <br />
+
+
+
+## R1X_Roccaraso_July_2014
+**Project name:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; R1X <br />
+**Year of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 2014 <br />
+**Location of Launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Roccaraso (ITA) <br />
+**Date of launch:** &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 07/2014 <br /> <br /> <br /> <br />
diff --git a/roccarasoTerrainData/zdata.mat b/roccarasoTerrainData/zdata.mat
index 9879fb3f3bab6a8c03b79eb4870c30289572c773..2c148ef7fe9ed43f473f54fc7ef2c8e5cdad7e46 100644
Binary files a/roccarasoTerrainData/zdata.mat and b/roccarasoTerrainData/zdata.mat differ