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