diff --git a/README.md b/README.md
index bf26e07e302b42aed284da9b9d111bdf06ef211a..5756d44a3e69d84a2f6beaa8d1402d6cb9ccf94e 100644
--- a/README.md
+++ b/README.md
@@ -1,90 +1,20 @@
-# Gemini Flight Logs new
+# Gemini Flight Logs
-## Getting started
+### Launch conditions
+* Launchpad Coordinates: 41.8084579 N, 14.0546408 E;
+* Launchpad Elevation: 79°
+* Launchpad Azimuth: 178°
+* Launch datetime: 19-Sep-2023 13:38:45
-To make it easy for you to get started with GitLab, here's a list of recommended next steps.
+### Atmospheric conditions
+* Wind
+ | Altitude agl [m] | 0-300 | 300-400 | 400-600 | 600-900 | 900-1300 | 1300-2000|
+ |:----------------:|:-----:|:-------:|:-------:|:-------:|:--------:|:--------:|
+ | Magnitude [m/s] | 3 | 6 | 7 | 9 | 12 | 14 |
+ | Azimuth [deg] | 300 | 270 | 270 | 270 | 270 | 270 |
+* Temperature: 22° C
+* Pressure: 86303 Pa
-Already a pro? Just edit this README.md and make it your own. Want to make it easy? [Use the template at the bottom](#editing-this-readme)!
-## Add your files
-- [ ] [Create](https://docs.gitlab.com/ee/user/project/repository/web_editor.html#create-a-file) or [upload](https://docs.gitlab.com/ee/user/project/repository/web_editor.html#upload-a-file) files
-- [ ] [Add files using the command line](https://docs.gitlab.com/ee/gitlab-basics/add-file.html#add-a-file-using-the-command-line) or push an existing Git repository with the following command:
-```
-cd existing_repo
-git remote add origin https://git.skywarder.eu/skyward/flight-logs/gemini-flight-logs-new.git
-git branch -M main
-git push -uf origin main
-```
-
-## Integrate with your tools
-
-- [ ] [Set up project integrations](https://git.skywarder.eu/skyward/flight-logs/gemini-flight-logs-new/-/settings/integrations)
-
-## Collaborate with your team
-
-- [ ] [Invite team members and collaborators](https://docs.gitlab.com/ee/user/project/members/)
-- [ ] [Create a new merge request](https://docs.gitlab.com/ee/user/project/merge_requests/creating_merge_requests.html)
-- [ ] [Automatically close issues from merge requests](https://docs.gitlab.com/ee/user/project/issues/managing_issues.html#closing-issues-automatically)
-- [ ] [Enable merge request approvals](https://docs.gitlab.com/ee/user/project/merge_requests/approvals/)
-- [ ] [Automatically merge when pipeline succeeds](https://docs.gitlab.com/ee/user/project/merge_requests/merge_when_pipeline_succeeds.html)
-
-## Test and Deploy
-
-Use the built-in continuous integration in GitLab.
-
-- [ ] [Get started with GitLab CI/CD](https://docs.gitlab.com/ee/ci/quick_start/index.html)
-- [ ] [Analyze your code for known vulnerabilities with Static Application Security Testing(SAST)](https://docs.gitlab.com/ee/user/application_security/sast/)
-- [ ] [Deploy to Kubernetes, Amazon EC2, or Amazon ECS using Auto Deploy](https://docs.gitlab.com/ee/topics/autodevops/requirements.html)
-- [ ] [Use pull-based deployments for improved Kubernetes management](https://docs.gitlab.com/ee/user/clusters/agent/)
-- [ ] [Set up protected environments](https://docs.gitlab.com/ee/ci/environments/protected_environments.html)
-
-***
-
-# Editing this README
-
-When you're ready to make this README your own, just edit this file and use the handy template below (or feel free to structure it however you want - this is just a starting point!). Thank you to [makeareadme.com](https://www.makeareadme.com/) for this template.
-
-## Suggestions for a good README
-Every project is different, so consider which of these sections apply to yours. The sections used in the template are suggestions for most open source projects. Also keep in mind that while a README can be too long and detailed, too long is better than too short. If you think your README is too long, consider utilizing another form of documentation rather than cutting out information.
-
-## Name
-Choose a self-explaining name for your project.
-
-## Description
-Let people know what your project can do specifically. Provide context and add a link to any reference visitors might be unfamiliar with. A list of Features or a Background subsection can also be added here. If there are alternatives to your project, this is a good place to list differentiating factors.
-
-## Badges
-On some READMEs, you may see small images that convey metadata, such as whether or not all the tests are passing for the project. You can use Shields to add some to your README. Many services also have instructions for adding a badge.
-
-## Visuals
-Depending on what you are making, it can be a good idea to include screenshots or even a video (you'll frequently see GIFs rather than actual videos). Tools like ttygif can help, but check out Asciinema for a more sophisticated method.
-
-## Installation
-Within a particular ecosystem, there may be a common way of installing things, such as using Yarn, NuGet, or Homebrew. However, consider the possibility that whoever is reading your README is a novice and would like more guidance. Listing specific steps helps remove ambiguity and gets people to using your project as quickly as possible. If it only runs in a specific context like a particular programming language version or operating system or has dependencies that have to be installed manually, also add a Requirements subsection.
-
-## Usage
-Use examples liberally, and show the expected output if you can. It's helpful to have inline the smallest example of usage that you can demonstrate, while providing links to more sophisticated examples if they are too long to reasonably include in the README.
-
-## Support
-Tell people where they can go to for help. It can be any combination of an issue tracker, a chat room, an email address, etc.
-
-## Roadmap
-If you have ideas for releases in the future, it is a good idea to list them in the README.
-
-## Contributing
-State if you are open to contributions and what your requirements are for accepting them.
-
-For people who want to make changes to your project, it's helpful to have some documentation on how to get started. Perhaps there is a script that they should run or some environment variables that they need to set. Make these steps explicit. These instructions could also be useful to your future self.
-
-You can also document commands to lint the code or run tests. These steps help to ensure high code quality and reduce the likelihood that the changes inadvertently break something. Having instructions for running tests is especially helpful if it requires external setup, such as starting a Selenium server for testing in a browser.
-
-## Authors and acknowledgment
-Show your appreciation to those who have contributed to the project.
-
-## License
-For open source projects, say how it is licensed.
-
-## Project status
-If you have run out of energy or time for your project, put a note at the top of the README saying that development has slowed down or stopped completely. Someone may choose to fork your project or volunteer to step in as a maintainer or owner, allowing your project to keep going. You can also make an explicit request for maintainers.
diff --git a/RoccarasoFlight/data/ANEMOMOMETER/anemometerData.mat b/RoccarasoFlight/data/ANEMOMOMETER/anemometerData.mat
new file mode 100644
index 0000000000000000000000000000000000000000..9b5bbf7af0487271490b243fa79bdc94b8cff6d0
Binary files /dev/null and b/RoccarasoFlight/data/ANEMOMOMETER/anemometerData.mat differ
diff --git a/RoccarasoFlight/postprocessing/AFD/configSimPostp.m b/RoccarasoFlight/postprocessing/AFD/configSimPostp.m
new file mode 100644
index 0000000000000000000000000000000000000000..9fe20883058a0db8d6d726e301fb609225137b35
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/AFD/configSimPostp.m
@@ -0,0 +1,192 @@
+%{
+
+CONFIG - This script sets up all the parameters for the simulation
+All the parameters are stored in the "settings" structure.
+
+REVISIONS:
+- #0 16/04/2016, Release, Francesco Colombi
+
+Copyright © 2021, Skyward Experimental Rocketry, AFD department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+%% MISSION FILE
+% Choose the mision you want to simulate from rocketsData folder
+
+
+%% LOAD DATA
+run("..\..\..\..\msa-toolkit\data\Gemini_Roccaraso_September_2023\simulationsData.m");
+
+%% LAUNCH SETUP
+% launchpad directions
+% for a single run the maximum and the minimum value of the following angles must be the same.
+settings.OMEGAmin = 79*pi/180; % [rad] Minimum Elevation Angle, user input in degrees (ex. 80)
+settings.OMEGAmax = 79*pi/180; % [rad] Maximum Elevation Angle, user input in degrees (ex. 80)
+settings.PHImin = 178*pi/180; % [rad] Minimum Azimuth Angle from North Direction, user input in degrees (ex. 90)
+settings.PHImax = 178*pi/180; % [rad] Maximum Azimuth Angle from North Direction, user input in degrees (ex. 90)
+
+% !! ATTENTION !! The following 2 work just for the stochastichs simulations with constant wind model only
+settings.upwind = false; % If true, phi is selected opposite to the wind direction
+settings.PHIsigma = 0*pi/180; % [deg] If upwind is true, you can select a variance for the direction
+
+%% ROCCARASO TERRAIN
+% !! ATTENTION !! the following works only at Roccaraso
+settings.terrain = false;
+
+if settings.terrain
+ settings.lat0 = 41.8084579; % Launchpad latitude
+ settings.lon0 = 14.0546408; % Launchpad longitude
+ settings.funZ = funZGen('zData.mat', settings.lat0, settings.lon0);
+end
+
+%% ENGINE SETTINGS
+% Following set the instant of time at which the engine is cut off. Set Inf
+% if you want to have a complete burn.
+
+if settings.HREmot
+ settings.timeEngineCut = (inf) * settings.tb; % [s] Moment in time in wich the engine will be cut off
+ settings = settingsEngineCut(settings); % updatind settings parameters
+end
+
+%% AEROBRAKES SETTINGS
+% Multiple air-brakes and smooth opening simulation
+settings.multipleAB = true; % If true, multiple and smooth airbrakes opening will be simulated
+
+% If FALSE:
+% - settings.control: only the first value will be computed;
+% - settings.dtControl: is not read.
+% If TRUE:
+% - settings.control: define the sequence of air-brakes opening
+% configuration. Closed air-brakes are simulated
+% untill the conditions settings.tControl and
+% settings.machControl are both verified, check
+% simulationsData.m;
+% - settings.dtControl: define the usage time of the i-th
+% configuration. Its length must be -1 the
+% length of settings.control
+
+settings.control = [1 3]; % aerobrakes, 1-2-3 for 0%, 50% or 100% opened
+settings.dtControl = res.time(indexOpeningARB) - res.time(indexCO); % aerobrakes, configurations usage time
+
+%% PARAFOIL SIMULATION
+settings.parafoilSim = false; % [-] True if parafoil the parafoli open loop descentn needs to be performed
+
+%% DESCENT PHASE MODEL
+settings.descent6DOF = false; % set to true in order to start a 6DOF parachute descent phase
+
+%% WIND DETAILS
+% select which model you want to use:
+
+%%%%% Matlab Wind Model
+settings.wind.model = false;
+% matlab hswm model, wind model on altitude based on historical data
+
+% input Day and Hour as arrays to run stochastic simulations
+settings.wind.DayMin = 105; % [d] Minimum Day of the launch
+settings.wind.DayMax = 105; % [d] Maximum Day of the launch
+settings.wind.HourMin = 4; % [h] Minimum Hour of the day
+settings.wind.HourMax = 4; % [h] Maximum Hour of the day
+settings.wind.ww = 0; % [m/s] Vertical wind speed
+
+%%%%% Input wind
+
+%%% Roccaraso 2022 wind
+% settings.wind.inputMagnitude = [3 * ones(1, 18), 7* ones(1, 15), 15* ones(1, 10), 19 * ones(1, 20), 21 * ones(1, 37)];
+% settings.wind.inputAzimut = [300 * ones(1, 10) , 250 * ones(1, 90)]; % wind azimut angle at each altitude (toward wind incoming direction) [deg]
+%---------------------------------------------------
+
+settings.wind.input = true;
+% Wind is generated for every altitude interpolating with the coefficient defined below
+
+settings.wind.inputGround = 3; % wind magnitude at the ground [m/s]
+settings.wind.inputAlt = linspace(0, 2000, 100); % altitude vector [m]
+
+settings.wind.inputMagnitude = [3 * ones(1,15), 6 * ones(1, 5), 7 * ones(1, 10), 9 * ones(1, 15), 12 * ones(1, 20), 14 * ones(1, 35)];
+settings.wind.inputAzimut = [300 * ones(1, 15), 270 * ones(1, 85)]; % wind azimut angle at each altitude (toward wind incoming direction) [deg]
+
+settings.wind.inputMult = (settings.wind.inputMagnitude./settings.wind.inputGround - 1) * 100;
+
+%settings.wind.inputAzimut = [300 * ones(1, 10) , 250 * ones(1, 90)]; % wind azimut angle at each altitude (toward wind incoming direction) [deg]
+%settings.wind.inputMult = [0 66 133 183 200 250]; % percentage of increasing magnitude at each altitude
+%settings.wind.inputAzimut = randi([170 280], 1, 100); % wind azimut angle at each altitude (toward wind incoming direction) [deg]
+
+settings.wind.inputUncertainty = [0, 0];
+% settings.wind.inputUncertainty = [a,b]; wind uncertanties:
+% - a, wind magnitude percentage uncertanty: magn = magn *(1 +- a)
+% - b, wind direction band uncertanty: dir = dir 1 +- b
+
+
+%%%%% Variable wind model
+settings.wind.variable = false;
+if settings.wind.variable
+ % Generate a gaussian distribution at ground level and
+ % a uniform distribution in altitude, linear interpolation between
+ % The azimuth in the high layer is +-deltaAz wrt the Azimut at ground
+ settings.wind.var.magMeanG = 4.7;
+ settings.wind.var.magSigmaG = 1.2840;
+ settings.wind.var.azMinG = 0*pi/180;
+ settings.wind.var.azMaxG = 359*pi/180;
+ settings.wind.var.hG = 100;
+ settings.wind.var.magMinH = 4;
+ settings.wind.var.magMaxH = 30;
+ settings.wind.var.deltaAzH = 5*pi/180;
+ settings.wind.var.hH = 1000;
+ % NOTE: wind azimuth angle indications (wind directed towards):
+ % 0 deg (use 360 instead of 0) -> North
+ % 90 deg -> East
+ % 180 deg -> South
+ % 270 deg -> West
+end
+
+%%%%% Random wind model
+% if all the model above are false
+
+% Wind is generated randomly from the minimum to the maximum parameters which defines the wind.
+% Setting the same values for min and max will fix the parameters of the wind.
+
+settings.wind.MagMin = 5; % [m/s] Minimum Magnitude
+settings.wind.MagMax = 5; % [m/s] Maximum Magnitude
+settings.wind.ElMin = 0*pi/180; % [rad] Minimum Elevation, user input in degrees (ex. 0)
+settings.wind.ElMax = 0*pi/180; % [rad] Maximum Elevation, user input in degrees (ex. 0) (Max == 90 Deg)
+settings.wind.AzMin = wrapTo360(180)*pi/180; % [rad] Minimum Azimuth, user input in degrees (ex. 90)
+settings.wind.AzMax = wrapTo360(180)*pi/180; % [rad] Maximum Azimuth, user input in degrees (ex. 90)
+
+% NOTE: wind azimuth angle indications (wind directed towards):
+% 0 deg (use 360 instead of 0) -> North
+% 90 deg -> East
+% 180 deg -> South
+% 270 deg -> West
+
+%% BALLISTIC SIMULATION
+% Set to True to run a ballistic (without parachutes) simulation
+settings.ballistic = false;
+
+%% PLOT DETAILS
+settings.plots = true;
+
+%% LANDING POINTS
+% satellite maps of the landing zone
+settings.landingMap = true; % 2D map
+settings.satellite3D = false; % 3D map
+
+%% OPEN ROCKET
+settings.openRocketSM = false;
+
+if settings.openRocketSM
+ openRocketSkywardPath = '../../OpenRocketSkyward/';
+ addpath(genpath(openRocketSkywardPath))
+ settings.mHaackFlag = isequal(settings.OgType, 'MHAACK');
+
+ vars.Mach = linspace(0.1, 0.9, 6);
+ vars.Alpha = sort([linspace(-20, 20, 17), -1, 1]);
+ vars.Beta = linspace(-20, 20, 20);
+ vars.Alt = linspace(0, 4000, 5);
+ vars.xcg = settings.xcg;
+ vars.hprot = settings.hprot;
+
+ settings.input = createDissileInput(settings, vars);
+end
+settings.stoch.N = 1;
diff --git a/RoccarasoFlight/postprocessing/AFD/estimationCA.m b/RoccarasoFlight/postprocessing/AFD/estimationCA.m
index a7e55b8a10119d8282516f527331ee3b77cabb04..d7f08131d7d65be148bd90a8ed5d58513e4c6e31 100644
--- a/RoccarasoFlight/postprocessing/AFD/estimationCA.m
+++ b/RoccarasoFlight/postprocessing/AFD/estimationCA.m
@@ -1,117 +1,149 @@
+function results = estimationCA(main, alt0)
+
+ %% DECLARE CONSTANTS
+ NPOINTS = 1000;
+ MSTRUCT = 25.2885; % [kg] Gemini structural mass
+ AREF = 0.15^2 * pi * 0.25;
+ PROP_FINAL_MASS = 0.7799 + 0.4445; % [kg] initial propellant mass
+ tIGN = 0.97500; % Ignition time in the thrust curve
+ tCO = 5.32702; % Cutoff time in the thrsut curve
+ mDot = 0.9660; % [kg/s] Averaged propellant mass flow rate
+ DELTA_TIME = 0.4988;
+ load("CA_alpha0.mat");
+
+ % acceleration raw
+ indexApoIMU = find(main.IMU(:, 1)> main.tApogee, 1, 'first');
+ accXB_raw = main.IMU(1:indexApoIMU ,2);
+ tAcc = main.IMU(1:indexApoIMU, 1);
+ accXB = movmean(accXB_raw, 30);
+
+ % Motor data
+ load("../RPS/cleanData/ROC2/engineFlightData.mat");
+ timeThrust = flightData.timeThrust;
+ indexTimeThrust = and(timeThrust > tIGN, timeThrust < tCO);
+ timeThrust = timeThrust(indexTimeThrust);
+ timeThrust = timeThrust - timeThrust(1) - DELTA_TIME;
+
+ Thrust = flightData.Thrust(indexTimeThrust);
+
+ tNodal = [timeThrust(1), -0.2248, 3.33322, timeThrust(end)];
+ indexTransIGN = and(timeThrust>=tNodal(1), timeThrust<tNodal(2));
+ indexTransCO = and(timeThrust>tNodal(3), timeThrust<=tNodal(4));
+
+ thrustTransIGN = Thrust(indexTransIGN);
+ thrustTransCO = Thrust(indexTransCO);
+
+ offsetIGN = thrustTransIGN(1);
+ offsetCO = thrustTransCO(end);
+
+ thrustTransIGN = thrustTransIGN - offsetIGN;
+ thrustTransCO = thrustTransCO - offsetCO;
+
+ thrustTransIGN = thrustTransIGN * (1 + offsetIGN/thrustTransIGN(end));
+ thrustTransCO = thrustTransCO * (1 + offsetCO/thrustTransCO(1));
+
+ Thrust = [thrustTransIGN; Thrust(and(timeThrust>=tNodal(2), timeThrust<=tNodal(3))); thrustTransCO];
+ mMotor = linspace(PROP_FINAL_MASS + mDot*(tCO - tIGN), PROP_FINAL_MASS, sum(indexTimeThrust));
+ timeThrust = timeThrust - timeThrust(1);
+
+ % Velocity
+ tNas = main.NASData(1:main.NASApogeeIndex, 1);
+ altNas = -main.NASData(1:main.NASApogeeIndex, 4);
+ velNas = vecnorm(main.NASData(1:main.NASApogeeIndex, 5:7), 2, 2);
+
+ load sensorPitot.mat
+ indNasPit = sensorTot_WithPitot.nas.time < main.tApogee;
+
+ tNasPit = sensorTot_WithPitot.nas.time(indNasPit);
+ velNasPit = vecnorm(sensorTot_WithPitot.nas.states(indNasPit, 4:6), 2,2);
+ altNasPit = -sensorTot_WithPitot.nas.states(indNasPit, 3);
+
+
+ % Normalize data
+ timeRef = linspace(max([tNas(1), tAcc(1), tNasPit(1)]), min([tNas(end), tAcc(end), tNasPit(end)]) , NPOINTS);
+ indexEndThrust = find(timeRef > timeThrust(end), 1, 'first') - 1;
+
+ accXBRef = interp1(tAcc, accXB, timeRef);
+
+ ThrustRef = zeros(1, NPOINTS);
+ ThrustRef(1:indexEndThrust) = interp1(timeThrust, Thrust, timeRef(1:indexEndThrust));
+
+ mRocket = (mMotor(end) + MSTRUCT) * ones(1, NPOINTS);
+ mRocket(1:indexEndThrust) = interp1(timeThrust, mMotor, timeRef(1:indexEndThrust)) + MSTRUCT;
+
+ velRef = interp1(tNas, movmean(velNas, 10), timeRef);
+
+ velNasPitRef = interp1(tNasPit, velNasPit, timeRef);
+
+ altRef = interp1(tNasPit, altNasPit, timeRef);
+
+ % Air density and sonic velocity
+ [~, vSonRef, ~, rhoRef] = atmosisa(altRef + alt0);
+ machRef = velRef./vSonRef;
+ machPitRef = velNasPitRef./vSonRef;
+
+ % index opening
+ indexOpeningARB = find(timeRef > 5.12, 1, 'first');
+
+ % Compute CA
+ CA = (ThrustRef - mRocket.*accXBRef) ./ (0.5 .* rhoRef .* velRef.^2 * AREF);
+
+ CA_NasPitot = (ThrustRef - mRocket.*accXBRef) ./ (0.5 .* rhoRef .* velNasPitRef.^2 * AREF);
+
+ tCO = tCO -1.8519;
+ indexCO = find(ThrustRef == 0, 1, 'first');
+
+ mRocketNoPwr = mRocket(indexCO:end);
+ CANoPwr = CA(indexCO:end);
+ CA_100ARB = zeros(1, length(timeRef));
+ CA_0ARB = zeros(1, length(timeRef));
+
+ for i = 1:length(timeRef)
+ [~, indexMach] = nearestValSorted(MachRef, machPitRef(i));
+ [~, indexAlt] = nearestValSorted(AltRef, altRef(i));
+ CA_100ARB(i) = CA_MD(indexMach, indexAlt, 3);
+ CA_0ARB(i) = CA_MD(indexMach, indexAlt, 1);
+ end
+
+ ThrustComp = 0.5 .* rhoRef(1:indexCO) .* velRef(1:indexCO).^2 .* AREF .* CA_0ARB(1:indexCO) + mRocket(1:indexCO).*accXBRef(1:indexCO);
+
+ figure
+ plot(timeRef(indexCO:end),CANoPwr);
+ hold on
+ plot(timeRef(indexCO:end), CA_NasPitot(indexCO:end));
+ plot(timeRef(indexCO:end), [CA_0ARB(indexCO:indexOpeningARB) CA_100ARB(indexOpeningARB+1:end)]);
+ xlabel('Time [s]');
+ ylabel('CA [-]');
+ title('CA from telemetry');
+ grid on
+ legend('NAS', 'NAS with PITOT', 'MD', 'Location','southoutside', 'Orientation', 'horizontal');
+
+ %% SAVE RESULTS
+
+ results.time = timeRef;
+ results.Thrust = ThrustRef;
+ results.ThrustComputed = ThrustComp;
+
+ results.velocity.vel_pitot = velNasPitRef;
+ results.velocity.val_noPitot = velRef;
+
+ results.alt = altRef;
+
+ results.CA.CA100 = CA_100ARB;
+ results.CA.CA0 = CA_0ARB;
+ results.CA.CA_pitot = CA_NasPitot;
+ results.CA.CA_noPitot = CA;
+
+ results.mass.propMass = mMotor;
+ results.mass.rocketMass = mRocket;
+
+ results.accXBody = accXBRef;
+
+ results.events.shutdonw = indexCO;
+ results.events.ARB_opening = indexOpeningARB;
-%% CA estimation
-NPOINTS = 1000;
-MSTRUCT = 25.2885; % [kg] Gemini structural mass
-AREF = 0.15^2 * pi * 0.25;
-PROP_FINAL_MASS = 0.7799 + 0.4445; % [kg] initial propellant mass
-tIGN = 0.97500; % Ignition time in the thrust curve
-tCO = 5.32702; % Cutoff time in the thrsut curve
-mDot = 0.9660; % [kg/s] Averaged propellant mass flow rate
-DELTA_TIME = 0.4988;
-load("CA_alpha0.mat");
-
-% acceleration raw
-indexApoIMU = find(main.IMU(:, 1)> main.tApogee, 1, 'first');
-accXB_raw = main.IMU(1:indexApoIMU ,2);
-tAcc = main.IMU(1:indexApoIMU, 1);
-accXB = movmean(accXB_raw, 30);
-
-% Motor data
-load("../RPS/cleanData/ROC-02/engineFlightData.mat");
-timeThrust = flightData.timeThrust;
-indexTimeThrust = and(timeThrust > tIGN, timeThrust < tCO);
-timeThrust = timeThrust(indexTimeThrust);
-timeThrust = timeThrust - timeThrust(1) - DELTA_TIME;
-
-Thrust = flightData.Thrust(indexTimeThrust);
-
-tNodal = [timeThrust(1), -0.2248, 3.33322, timeThrust(end)];
-indexTransIGN = and(timeThrust>=tNodal(1), timeThrust<tNodal(2));
-indexTransCO = and(timeThrust>tNodal(3), timeThrust<=tNodal(4));
-
-thrustTransIGN = Thrust(indexTransIGN);
-thrustTransCO = Thrust(indexTransCO);
-
-offsetIGN = thrustTransIGN(1);
-offsetCO = thrustTransCO(end);
-
-thrustTransIGN = thrustTransIGN - offsetIGN;
-thrustTransCO = thrustTransCO - offsetCO;
-
-thrustTransIGN = thrustTransIGN * (1 + offsetIGN/thrustTransIGN(end));
-thrustTransCO = thrustTransCO * (1 + offsetCO/thrustTransCO(1));
-
-Thrust = [thrustTransIGN; Thrust(and(timeThrust>=tNodal(2), timeThrust<=tNodal(3))); thrustTransCO];
-mMotor = linspace(PROP_FINAL_MASS + mDot*(tCO - tIGN), PROP_FINAL_MASS, sum(indexTimeThrust));
-timeThrust = timeThrust - timeThrust(1);
-
-% Velocity
-tNas = main.NASData(1:main.NASApogeeIndex, 1);
-altNas = -main.NASData(1:main.NASApogeeIndex, 4);
-velNas = vecnorm(main.NASData(1:main.NASApogeeIndex, 5:7), 2, 2);
-
-tPitot = main.PITOT(:, 1);
-vPitot = main.PITOT(:, 3);
-
-% Air density
-[~, vSon, ~, rho] = atmosisa(altNas + alt0);
-
-% Mach number
-mach = velNas./vSon;
-
-% Normalize data
-timeRef = linspace(max([tNas(1), tAcc(1), tPitot(1)]), min([tNas(end), tAcc(end), tPitot(end)]) , NPOINTS);
-indexEndThrust = find(timeRef > timeThrust(end), 1, 'first');
-
-accXBRef = interp1(tAcc, accXB, timeRef);
-
-ThrustRef = zeros(1, NPOINTS);
-ThrustRef(1:indexEndThrust) = interp1(timeThrust, Thrust, timeRef(1:indexEndThrust));
-
-mRocket = (mMotor(end) + MSTRUCT) * ones(1, NPOINTS);
-mRocket(1:indexEndThrust) = interp1(timeThrust, mMotor, timeRef(1:indexEndThrust)) + MSTRUCT;
-
-rhoRef = interp1(tNas, rho, timeRef);
-
-velRef = interp1(tNas, movmean(velNas, 10), timeRef);
-
-velPitRef = interp1(tPitot, vPitot, timeRef);
-
-altRef = interp1(tNas, altNas, timeRef);
-
-% Compute CA
-CA = (ThrustRef - mRocket.*accXBRef) ./ (0.5 .* rhoRef .* velRef.^2 * AREF);
-
-CA_pitot = (ThrustRef - mRocket.*accXBRef) ./ (0.5 .* rhoRef .* velPitRef.^2 * AREF);
-
-tCO = tCO -1.8519;
-indexCO = find(ThrustRef == 0, 1, 'first');
-
-mRocketNoPwr = mRocket(indexCO:end);
-CANoPwr = CA(indexCO:end);
-CA_100ARB = zeros(1, length(CANoPwr));
-CA_0ARB = zeros(1, length(CANoPwr));
-
-for i = 1:length(CANoPwr)
- [~, indexMach] = nearestValSorted(MachRef, mach(i));
- [~, indexAlt] = nearestValSorted(AltRef, altRef(i));
- CA_100ARB(i) = CA_MD(indexMach, indexAlt, 3);
- CA_0ARB(i) = CA_MD(indexMach, indexAlt, 1);
end
-figure
-plot(timeRef(indexCO:end),CANoPwr);
-hold on
-plot(timeRef(indexCO:end), CA_pitot(indexCO:end));
-plot(timeRef(indexCO:end), CA_100ARB);
-plot(timeRef(indexCO:end), CA_0ARB);
-xlabel('Time [s]');
-ylabel('CA [-]');
-title('CA from telemetry');
-grid on
-legend('NAS', 'PITOT', '100% ARB', '0% ARB', 'Location','southoutside', 'Orientation', 'horizontal');
-
-
diff --git a/RoccarasoFlight/postprocessing/AFD/mainPostprocess.m b/RoccarasoFlight/postprocessing/AFD/mainPostprocess.m
index 8ce83195d867c22c18b6c0ed67be45b58b09dd0f..3616d5a4580deab3333a7d604fcc6a5dae1556b5 100644
--- a/RoccarasoFlight/postprocessing/AFD/mainPostprocess.m
+++ b/RoccarasoFlight/postprocessing/AFD/mainPostprocess.m
@@ -9,7 +9,7 @@ alt0 = 1414;
%% adding path and recall data
addpath("..\commonFunctions\");
-addpath("..\..\..\..\msa-toolkit\commonFunctions\miscellaneous\");
+addpath(genpath("..\..\..\..\msa-toolkit\commonFunctions\"));
main = importData('main');
%% Trajectory
@@ -21,8 +21,47 @@ end
%% CA Estimation
if opt.flagCA
- estimationCA;
+ res = estimationCA(main,alt0);
+
end
+%% TEST Simulation
+
+if opt.flagCA
+
+ indexCO = res.events.shutdonw;
+ indexOpeningARB = res.events.ARB_opening;
+
+ configSimPostp;
+
+ % CONTINUE
+ Thrust = res.Thrust(1:indexCO);
+ Thrust(end) = 0;
+ timeThrust = res.time(2:indexCO);
+ timeThrust = [0, timeThrust];
+
+ settings.motor.expM = res.mass.propMass(1:indexCO);
+ settings.motor.Pe = reshapeVector(timeThrust, settings.motor.expTime, settings.motor.Pe);
+ settings.mTotalTime = res.mass.rocketMass(1:indexCO);
+ settings.ms = res.mass.rocketMass(end);
+ settings.xcg = reshapeVector(timeThrust, settings.motor.expTime, settings.xcg);
+ settings.Ixx = reshapeVector(timeThrust, settings.motor.expTime, settings.Ixx);
+ settings.Iyy = reshapeVector(timeThrust, settings.motor.expTime, settings.Iyy);
+ settings.Izz = reshapeVector(timeThrust, settings.motor.expTime, settings.Izz);
+ settings.I = [settings.Ixx; settings.Iyy; settings.Izz];
+
+ %%% Inertias derivatives (avoiding calculations in ascent.m)
+ settings.Idot = diff(settings.I')'./diff(timeThrust);
+ settings.Idot(:, end+1) = settings.Idot(:, end);
+
+ settings.State.xcgTime = settings.State.xcgTime * timeThrust(end)/settings.State.xcgTime(end);
+
+ settings.motor.expTime = timeThrust;
+ settings.motor.expThrust = Thrust;
+ settings.tb = settings.motor.expTime(end);
+ settings.tControl = settings.tb;
+
+ [Tf, Yf, Ta, Ya, bound_value] = stdRun(settings);
+end
diff --git a/RoccarasoFlight/postprocessing/AFD/sensorPitot.mat b/RoccarasoFlight/postprocessing/AFD/sensorPitot.mat
new file mode 100644
index 0000000000000000000000000000000000000000..3519293201d763415582910d6409c41a8273a728
Binary files /dev/null and b/RoccarasoFlight/postprocessing/AFD/sensorPitot.mat differ
diff --git a/RoccarasoFlight/postprocessing/RPS/cleanData/ROC2/engineFlightData.mat b/RoccarasoFlight/postprocessing/RPS/cleanData/ROC2/engineFlightData.mat
index cfd6deb231917c80950f74c2dbfb8ad51dc084f0..aef81bab5ef6ee0f4df703b96c10110556d95a46 100644
Binary files a/RoccarasoFlight/postprocessing/RPS/cleanData/ROC2/engineFlightData.mat and b/RoccarasoFlight/postprocessing/RPS/cleanData/ROC2/engineFlightData.mat differ
diff --git a/RoccarasoFlight/postprocessing/RPS/configDataAnalyzer.m b/RoccarasoFlight/postprocessing/RPS/configDataAnalyzer.m
index 5e45f0fb52fdfdf25262ec17e7b29681e7982524..441dcc97598c2c3f464aa411ec2b11bf716684e7 100644
--- a/RoccarasoFlight/postprocessing/RPS/configDataAnalyzer.m
+++ b/RoccarasoFlight/postprocessing/RPS/configDataAnalyzer.m
@@ -1,10 +1,9 @@
%% FLIGHT and ACQUISITION SETTINGS
-settings.launchSite = 'Roccaraso'; % [-] Launch site name
settings.testName = 'ROC2'; % [-] Test name
settings.frPT = 1000; % [Hz] PTs sampling rate
settings.frSrv = 10; % [Hz] SRVs sampling rate
-settings.fileCSVpath = strcat('../../../', settings.launchSite, '/data/MOTORE/');
+settings.fileCSVpath = '../../data/MOTORE/';
% Load mission parameters
run(strcat('engineData/missions/', settings.testName));
@@ -16,6 +15,20 @@ settings.engineName = 'Furia';
% Retrieve engine data
run(strcat('engineData/engines/', settings.engineName, '/', settings.engineName));
+% Load CEA database
+load('engineData/engines/CEA_db.mat');
+engine.CEA_db = CEA_db;
+clear CEA_db
+
+%% SIMULATION SETTINGS
+settings.interpMthd = 1;
+
+%% CONSTANTS
+settings.g0 = 9.806; % [m/s2] Acceleration of Gravity, constant
+settings.R = 8314.51; % [kJ/mol-K]
+settings.z = 1410; % [m] Altitude above-sea-level
+[~, ~, settings.Pamb, ~] = atmosisa(settings.z); % [Pa] Ambient initial pressure
+
%% ODE SETTINGS
% !! Do not change unless you know what you are doing !!
settings.odeEngineSettings = odeset('RelTol', 1e-5, 'AbsTol', 1e-5, 'Events', @eventBurnOut);
@@ -28,5 +41,4 @@ settings.bisection.xTol = 1e-3;
settings.bisection.absTol = 1e-3;
%% PLOT SETTINGS
-settings.plots = true; % [-] True if you want to see the plots
-
+settings.plots = true; % [-] True if you want to see the plots
\ No newline at end of file
diff --git a/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m b/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m
index b3f9c00e32f4fa8641288d81ef355cb38d87864e..ca73d482882adb77bf26974dfaa359276ef2373b 100644
--- a/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m
+++ b/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m
@@ -30,9 +30,20 @@ end
%% DATA PROCESSING AND ANALYSIS
engine.acquisition = flightData;
-engine = tankEvacuation(engine, settings);
+% Computing the mass flow rate of oxidizer
+engine = tankEvacuation(engine);
+
+% Simulating the engine
+engineSim = combustionSimulation(engine, settings);
+
+% Computing actual performance
+engineFly = performanceReconstruction(engine, engineSim);
%% PLOTTING
+fprintf('DATA ANALYSIS CONCLUDED: \n\n')
+fprintf('\t - Total Impulse simulated: %.2f\n', engineSim.performances.Itot)
+fprintf('\t - Total Impulse computed: %.2f\n', engineFly.Itot)
+
if settings.plots
f1 = figure;
yyaxis left
@@ -45,6 +56,21 @@ if settings.plots
stairs(flightData.timeSrvMain, flightData.srvMain)
xlabel('Time [s]'); ylabel('SRV Opening [-]'); title('Tank pressure values')
legend('Top tank', 'Bottom tank', 'Combustion chamber')
+ axis padded
+
+ f2 = figure;
+ plot(engineSim.time, engineSim.mfr.mp)
+ grid minor
+ ylabel('Mass Flow [kg/s]'); xlabel('Time [s]'); title('Propellant MFR')
+ axis padded
+
+ f3 = figure;
+ plot(engineSim.time, engineSim.performances.T)
+ hold on; grid minor
+ plot(engineFly.time, engineFly.T)
+ ylabel('Thrust [N]'); xlabel('Time [s]'); title('Thrust compare')
+ legend('Simulated', 'Computed')
+ axis padded
end
diff --git a/RoccarasoFlight/postprocessing/RPS/engineData/engines/CEA_db.mat b/RoccarasoFlight/postprocessing/RPS/engineData/engines/CEA_db.mat
index 2712efdac18d9d0c55b4166d865a2d70d69286e8..92de427b9ecb77f7f42a48eb7eda47fdef275435 100644
Binary files a/RoccarasoFlight/postprocessing/RPS/engineData/engines/CEA_db.mat and b/RoccarasoFlight/postprocessing/RPS/engineData/engines/CEA_db.mat differ
diff --git a/RoccarasoFlight/postprocessing/RPS/engineData/engines/Furia/Furia.m b/RoccarasoFlight/postprocessing/RPS/engineData/engines/Furia/Furia.m
index 09f094e2b292763bb2615892ff7fae2645ea4120..7190b74c04cf1b46d10ba6f40da6562379955075 100644
--- a/RoccarasoFlight/postprocessing/RPS/engineData/engines/Furia/Furia.m
+++ b/RoccarasoFlight/postprocessing/RPS/engineData/engines/Furia/Furia.m
@@ -44,10 +44,12 @@ engine.data.d_pipe = 1.25e-2; % [m] Diameter of the feedline pipes
engine.data.dInj = 1.5e-3; % [m] Diameter of injection holes
engine.data.Aox = engine.data.Ninj* ...
pi/4*engine.data.dInj.^2; % [m^2] Injection area
-
+engine.data.alphaFit = [-11.7668580891917, 21.5544619916252, -8.85332436329970]; % [-] Scaling factor between top tank pressure and injector pressure
+engine.data.alphaDeg = length(engine.data.alphaFit) - 1;
+
%%% EFFICIENCIES
-engine.data.etaCT = 0.96;
-engine.data.etaCs = 0.85;
+engine.data.etacT = 0.97;
+engine.data.etaCs = 0.84;
%%% COMBUSTION CHAMBER
engine.data.rt = 0.0149; % [m] Throat radius
@@ -56,6 +58,17 @@ engine.data.lambda = 0.985; % [-] Loss due to nozzle divergence
engine.data.Rc_CM = 0.5;
engine.data.R_pre = 1.15; % [-] Ratio between pre-chamber length and port radius
engine.data.R_post = 2.35; % [-] ratio between post-chamber length and port radius
+engine.data.Pc0 = 10; % [bar] Initial pressure value
+
+%%% PRE CC
+engine.data.Dpre = 66e-3;
+engine.data.Lpre = 29.3e-3;
+engine.data.Vpre = pi/4*engine.data.Dpre^2*engine.data.Lpre;
+
+%%% POST CC
+engine.data.Dpost = 62e-3;
+engine.data.Lpost = 60e-3;
+engine.data.Vpost = pi/4*engine.data.Dpost^2*engine.data.Lpost;
%%% NOZZLE
engine.data.zopt = 0; % [m] Altitude of optimal expansion
diff --git a/RoccarasoFlight/postprocessing/RPS/src/Ap2RegressionCoor.m b/RoccarasoFlight/postprocessing/RPS/src/Ap2RegressionCoor.m
new file mode 100644
index 0000000000000000000000000000000000000000..4324600f74fb39d8559fb2f0c19dff1ca67cd459
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/RPS/src/Ap2RegressionCoor.m
@@ -0,0 +1,30 @@
+function [r, dr, Ab] = Ap2RegressionCoor(Ap, dAp, engine)
+%{
+Ap2RegressionCoor - this function takes the port area of the grain and
+gives back the regression coordinate and its derivative. This is done
+accordin to the grain type specified in the structure engine.
+
+CALLED SCRIPTS:
+
+CALLED FUNCTIONS:
+
+CALLED DATA FILES:
+
+REVISIONS:
+- 0, 12/09/2022, Release, Alberto Boffi, Paolo Gnata, PRP
+Department
+
+Copyright © 2022, Skyward Experimental Rocketry, PRP department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+%% INPUT
+L = engine.data.Lgrain;
+
+%% CALCULATIONS
+r = sqrt(Ap./pi);
+dr = dAp./(2*pi.*r);
+Ab = L*2*pi.*r;
\ No newline at end of file
diff --git a/RoccarasoFlight/postprocessing/RPS/src/combustionSimulation.m b/RoccarasoFlight/postprocessing/RPS/src/combustionSimulation.m
new file mode 100644
index 0000000000000000000000000000000000000000..19cad597d0dd54d74ed290e02fa2abc746045092
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/RPS/src/combustionSimulation.m
@@ -0,0 +1,136 @@
+function engineSim = combustionSimulation(engine, settings)
+%{
+combustionCompare - this function simulates combustion chamber behavior
+using estimated/measured mox curve and compares it with those resulting
+from combustion with most reliable a, n coefficients retrieved until now
+
+CALLED SCRIPTS: /
+
+CALLED FUNCTIONS: /
+
+CALLED DATA FILES: /
+
+REVISIONS:
+-1, 20/04/2023, Release, Alberto Boffi, PRP Department
+
+Copyright © 2023, Skyward Experimental Rocketry, PRP Department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+%% INPUT
+
+%%% GRAIN
+L = engine.data.Lgrain;
+
+%%% ENGINE
+epsilon = engine.data.epsilon;
+At = pi * engine.data.rt^2;
+lambda = engine.data.lambda;
+etacT = engine.data.etacT;
+
+%%% PROPELLANT
+rhoFu = engine.data.fu.rhof;
+
+%%% ENGINE GUESSES
+Ap0 = pi*engine.data.r0^2;
+Pc0 = engine.data.Pc0*1e5;
+
+%%% SETTINGS
+fsolveOptions = settings.fsolveSettings;
+tSpan = engine.acquisition.cutData.timePTop - engine.acquisition.cutData.timePTop(1);
+Pamb = settings.Pamb * 1e-5;
+switch settings.interpMthd
+ case 1
+ interpFun = 'interp1';
+ case 3
+ interpFun = 'spline';
+ otherwise
+ error(['\nThe chosen interpolation method cannot be used in this code.\n ...' ...
+ 'Look at configTestSimulator.m for the allowed options'])
+end
+
+%%% CONSTANTS
+g0 = settings.g0;
+
+
+%% COMBUSTION SIMULATION
+% Defining the time
+engineSim.time = engine.acquisition.cutData.timePTop - engine.acquisition.cutData.timePTop(1);
+engineSim.mfr.mox = engine.performances.mox;
+
+% Performing Marxman simulation
+[solution] = ode45(@(t, Ap) marxmanRegressionOde(t, Ap, engine, engineSim, interpFun), tSpan, Ap0);
+
+% Evaluating solution derivative
+[Ap, dAp] = deval(tSpan, solution);
+engineSim.cc.Ap = Ap';
+engineSim.cc.dAp = dAp';
+clear solution
+
+% Converting into radial coordinate
+[~, dr, ~] = Ap2RegressionCoor(Ap, dAp, engine);
+
+
+%% CHAMBER MASS BALANCE & PERFORMANCES
+% Computing fuel mass flow rate
+engineSim.mfr.mf = L*dAp'*rhoFu;
+
+% Computing oxidizer-to-fuel ratio
+OF = engineSim.mfr.mox./engineSim.mfr.mf;
+OF(1) = OF(2);
+engineSim.performances.OF = OF;
+
+% Integration of mass balance equation to retrieve pressure
+[solution] = ode45(@(t, Pc) massBalance(t, Pc, engine, engineSim, settings, interpFun), tSpan, Pc0);
+
+% Evaluating solution derivative
+[Pc, ~] = deval(tSpan, solution);
+
+% Conversion from Pa to bar
+Pc = Pc'*1e-5;
+clear solution
+
+
+%% PERFORMANCES
+% Gamma
+y = interp2(engine.CEA_db.OF_input, engine.CEA_db.P_input, engine.CEA_db.gamma', OF, Pc, 'spline');
+
+% Exit pressure
+PeFun = @(Pe, Pc, y) -1./epsilon + ((y + 1)./2).^(1./(y - 1)) .* ((Pe/Pc).^(1./y)) .* sqrt((y + 1)./(y - 1).*(1 - (Pe./Pc).^((y - 1)./y)));
+Pe = [];
+for ii = 1:length(Pc)
+ jj = Pc(ii);
+ kk = y(ii);
+ Pe = [Pe fsolve(@(Pe) PeFun(Pe, jj, kk), Pamb, fsolveOptions)];
+end
+
+% Turning the array
+Pe = Pe';
+
+% Thrust coefficient
+cT = sqrt(2*y.^2./(y - 1).*(2./(y + 1)).^((y + 1)./(y - 1))).*sqrt(1 - (Pe./Pc).^((y - 1)./y)) * lambda + (Pe - Pamb)./Pc.*epsilon;
+
+% Thrust
+T = cT.*Pc*1e5*At*etacT;
+
+% Specific impulse
+Isp = T./engineSim.mfr.mox/g0;
+
+% Total impulse
+Itot = trapz(engineSim.time, T);
+
+
+%% OUTPUT
+engineSim.cc.dr = dr;
+engineSim.cc.Pc = Pc;
+engineSim.performances.Pe = Pe;
+engineSim.performances.T = T;
+engineSim.performances.cT = cT;
+engineSim.performances.Isp = Isp;
+engineSim.performances.Itot = Itot;
+engineSim.mfr.mp = engineSim.mfr.mox + engineSim.mfr.mf;
+
+
diff --git a/RoccarasoFlight/postprocessing/RPS/src/dataLoader.m b/RoccarasoFlight/postprocessing/RPS/src/dataLoader.m
index e884a7dd4414bcc1dd3180d3f5ab78be24f728a8..253274bc6057b1794967114771f7c2cc552b8f9e 100644
--- a/RoccarasoFlight/postprocessing/RPS/src/dataLoader.m
+++ b/RoccarasoFlight/postprocessing/RPS/src/dataLoader.m
@@ -53,7 +53,7 @@ if strcmp(settings.testName, 'ROC1')
elseif strcmp(settings.testName, 'ROC2')
% Defining the burning time
- tb = 3.3;
+ tb = 4.0;
% Defining the first index
ip = 8481 * settings.frPT;
@@ -75,6 +75,17 @@ elseif strcmp(settings.testName, 'ROC2')
flightData.PEng = rawData.PEng(ip:jp);
flightData.srvMain = rawData.srvMain(is:js);
flightData.srvVent = rawData.srvVent(is:js);
+
+
+ % Defining manual cutting indeces
+ [~, ic] = min(abs(0.980 - flightData.timePTop));
+ [~, jc] = min(abs(0.980 + tb - flightData.timePTop));
+
+ flightData.cutData.timePTop = flightData.timePTop(ic:jc);
+ flightData.cutData.timePEng = flightData.timePEng(ic:jc);
+
+ flightData.cutData.PTankTop = flightData.PTankTop(ic:jc);
+ flightData.cutData.PEng = flightData.PEng(ic:jc);
end
% Fixing pressures
diff --git a/RoccarasoFlight/postprocessing/RPS/src/marxmanRegressionOde.m b/RoccarasoFlight/postprocessing/RPS/src/marxmanRegressionOde.m
new file mode 100644
index 0000000000000000000000000000000000000000..e39c70dcf2181a2c8460f30cc181f439448c20a3
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/RPS/src/marxmanRegressionOde.m
@@ -0,0 +1,11 @@
+function dAp = marxmanRegressionOde(t, Ap, engine, engineSim, interpFun)
+
+%%% INPUT
+a = engine.data.fu.a;
+n = engine.data.fu.n;
+
+%%% INTERPOLATION
+mox = eval(strcat(interpFun, '(engineSim.time, engineSim.mfr.mox, t)'));
+
+% Port area derivative
+dAp = 2*a*sqrt(pi.*Ap).*(mox./Ap).^n;
\ No newline at end of file
diff --git a/RoccarasoFlight/postprocessing/RPS/src/massBalance.m b/RoccarasoFlight/postprocessing/RPS/src/massBalance.m
new file mode 100644
index 0000000000000000000000000000000000000000..cccae0b178926426dd03cfc3d989501114729b20
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/RPS/src/massBalance.m
@@ -0,0 +1,46 @@
+function dPc = massBalance(t, Pc, engine, engineSim, settings, interpFun)
+
+%% INPUT
+
+%%% GRAIN
+Vpre = engine.data.Vpre;
+Vpost = engine.data.Vpost;
+L = engine.data.Lgrain;
+etacS = engine.data.etaCs;
+
+%%% NOZZLE
+At = pi * engine.data.rt^2;
+
+%%% CONSTANTS
+R = settings.R;
+
+
+%% INTERPOLATION
+mox = eval(strcat(interpFun, '(engineSim.time, engineSim.mfr.mox, t)'));
+mf = eval(strcat(interpFun, '(engineSim.time, engineSim.mfr.mf, t)'));
+dAp = eval(strcat(interpFun, '(engineSim.time, engineSim.cc.dAp, t)'));
+Ap = eval(strcat(interpFun, '(engineSim.time, engineSim.cc.Ap, t)'));
+
+if t == 0
+ OF = engineSim.performances.OF(2);
+else
+ OF = mox/mf;
+end
+
+% CEA thermodynamics parameters
+% CC temperature
+Tc = interp2(engine.CEA_db.OF_input, engine.CEA_db.P_input, engine.CEA_db.T_comb', OF, Pc*1e-5, 'spline');
+
+% Gamma
+y = interp2(engine.CEA_db.OF_input, engine.CEA_db.P_input, engine.CEA_db.gamma', OF, Pc*1e-5, 'spline');
+
+% Molecular mass
+mwc = interp2(engine.CEA_db.OF_input, engine.CEA_db.P_input, engine.CEA_db.mw_comb', OF, Pc*1e-5, 'spline');
+
+% Characteristic velocity
+cS = sqrt(R./mwc.*Tc)./(sqrt(y .* (2./(y + 1)).^((y + 1)./(y - 1))));
+
+
+%% EQUATION
+% chamber pressure
+dPc = (mox + mf - Pc*At/etacS/cS - Pc*dAp*L/R/mwc/Tc/etacS^2) * (R/mwc*Tc*etacS^2/(Ap*L+ Vpre + Vpost));
\ No newline at end of file
diff --git a/RoccarasoFlight/postprocessing/RPS/src/performanceReconstruction.m b/RoccarasoFlight/postprocessing/RPS/src/performanceReconstruction.m
new file mode 100644
index 0000000000000000000000000000000000000000..7ceb117fbee5db8d8ea6b30678dda827bede730c
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/RPS/src/performanceReconstruction.m
@@ -0,0 +1,51 @@
+function engineFly = performanceReconstruction(engine, engineSim)
+%{
+combustionCompare - this function simulates combustion chamber behavior
+using estimated/measured mox curve and compares it with those resulting
+from combustion with most reliable a, n coefficients retrieved until now
+
+CALLED SCRIPTS: /
+
+CALLED FUNCTIONS: /
+
+CALLED DATA FILES: /
+
+REVISIONS:
+-1, 24/09/2023, Release, Paolo Gnata, PRP Department
+
+Copyright © 2023, Skyward Experimental Rocketry, PRP Department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+%% INPUT
+
+%%% DATA
+time = engine.acquisition.cutData.timePEng - engine.acquisition.cutData.timePEng(1);
+Pc = engine.acquisition.cutData.PEng;
+
+%%% PARAMETERS
+cT = engineSim.performances.cT;
+etaCT = engine.data.etacT;
+At = pi * engine.data.rt^2;
+
+
+%% RESHAPING
+% Interpolating the thrust coefficient
+cT = interp1(engineSim.time, cT, time);
+
+
+%% PERFORMANCES
+% Computing the thrust
+T = Pc .* cT * At * etaCT * 1e5;
+
+% Computing the total impulse
+Itot = trapz(time, T);
+
+
+%% OUTPUT
+engineFly.time = time;
+engineFly.T = T;
+engineFly.Itot = Itot;
\ No newline at end of file
diff --git a/RoccarasoFlight/postprocessing/RPS/src/tankEvacuation.m b/RoccarasoFlight/postprocessing/RPS/src/tankEvacuation.m
index 8e90e499b58e0fb049d3800c89be0a395b724464..89e9769dacf413a1737cc0b58d072067c8a0c694 100644
--- a/RoccarasoFlight/postprocessing/RPS/src/tankEvacuation.m
+++ b/RoccarasoFlight/postprocessing/RPS/src/tankEvacuation.m
@@ -1,4 +1,4 @@
-function engine = tankEvacuation(engine, settings)
+function engine = tankEvacuation(engine)
%{
tankEvacuation - this function outputs oxydizer mass flow rate usinìg NHNE
model
@@ -23,11 +23,12 @@ SPDX-License-Identifier: GPL-3.0-or-later
d_pipe = engine.data.d_pipe;
N2O = engine.data.ox;
-pTank = engine.acquisition.PTankTop;
+pTank = engine.acquisition.cutData.PTankTop;
Aox = engine.data.Aox;
cd = engine.data.cd_inj;
-alpha = 1;
-time = engine.acquisition.timePTop;
+alphaFit = engine.data.alphaFit;
+
+alphaFun = @(x) sum(alphaFit .* x.^(engine.data.alphaDeg:-1:0), 2);
%% CALCULATIONS
@@ -40,7 +41,9 @@ v = 10;
for i = 1 : length(p1)
rho1(i) = interp1(N2O.P, N2O.rhoL, p1(i));
- p1(i) = pTank(i) - (alpha/(2*Aox.^2*rho1(i)).*mox(i).^2)*1e-5;
+ if i > 1
+ p1(i) = pTank(i) - alphaFun(mox(i - 1)) .* (0.5*(mox(i - 1)/Aox).^2)./rho1(i) * 1e-5;
+ end
s1(i) = interp1(N2O.P, N2O.sL, p1(i));
h1(i) = interp1(N2O.P, N2O.hL, p1(i));
p2(i) = p1(i) * N2O.eta_cr(p1(i));
@@ -55,15 +58,15 @@ for i = 1 : length(p1)
rhoV2(i) = interp1(N2O.P, N2O.rhoV, p2(i));
rho2(i) = 1/(x2(i) * 1/rhoV2(i) + (1-x2(i))*1/rhoL2(i));
k(i) = sqrt((p1(i)-p2(i))./(pTank(i)-p2(i)));
- mox(i) = Aox * cd * (k(i)*sqrt(2*abs(p1(i)-p2(i))*1e5*rho1(i)) + ...
+ mox(i + 1) = Aox * cd * (k(i)*sqrt(2*abs(p1(i)-p2(i))*1e5*rho1(i)) + ...
rho2(i)*sqrt(2*abs(h1(i)-h2(i))))/(1+k(i));
- v(i+1) = mox(i)./rho1(i)/(pi*d_pipe^2/4);
+ v(i + 1) = mox(i)./rho1(i)/(pi*d_pipe^2/4);
end
%% OUTPUT
% ozydizer mass flow rate
-engine.performances.mox = mox;
+engine.performances.mox = mox(1:end - 1);