CA postprocessing - WIP

RoccarasoFlight/postprocessing/AFD/estimationCA.m

+
+%% 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'); 
+
+
+
+
+

RoccarasoFlight/postprocessing/AFD/mainPostprocess.m

@@ -2,46 +2,27 @@ clear;
 clc; 
 close all; 
 
+%% Select options
+opt.flagTraj = false; 
+opt.flagCA = true; 
+alt0 = 1414; 
+
 %% adding path and recall data
 addpath("..\commonFunctions\"); 
+addpath("..\..\..\..\msa-toolkit\commonFunctions\miscellaneous\"); 
 main = importData('main'); 
 
+%% Trajectory
 
-%% 2D TRAJECTORY - GPS
-tGps = main.GPS(:, 1); 
-lat = main.GPS(:, 2); 
-lon = main.GPS(:, 3); 
-altGps = main.GPS(:, 4); 
-
-indApo = find(tGps > main.tApogee, 1, 'first');
-
-figure
-geoplot(lat, lon, 'LineWidth', 2); 
-geobasemap satellite
-hold on
-geoplot(lat(1), lon(1), 'Marker','o', 'MarkerSize', 10, 'LineWidth',2); 
-geoplot(lat(indApo), lon(indApo),  'Marker','square', 'MarkerSize', 10, 'LineWidth',2); 
-geoplot(lat(end), lon(end), 'Marker','x', 'MarkerSize', 10, 'LineWidth',2, 'Color', 'r'); 
-legend('Trajectory', 'Liftoff', 'Apogee', 'Landing'); 
-
-
-%% 3D TRAJECTORY - GPS
-alt = -main.NASData(:, 4); 
-tNas = main.NASData(:, 1); 
-altNasInterp = interp1(tNas, alt, tGps); 
-alt0 = atmospalt(main.RCS_PRESSURE(1, 2)); 
-
-uif = uifigure;
-g = geoglobe(uif);
-geoplot3(g,lat,lon,altNasInterp + alt0,'b','Linewidth',2.5);
-
-%% CA estimation 
-
-
-
-
+if opt.flagTraj 
+    trajectory; 
+end
 
+%% CA Estimation
 
+if opt.flagCA
+    estimationCA; 
+end
 
 
 

RoccarasoFlight/postprocessing/AFD/trajectory.m

+%% 2D TRAJECTORY - GPS
+tGps = main.GPS(:, 1); 
+lat = main.GPS(:, 2); 
+lon = main.GPS(:, 3); 
+altGps = main.GPS(:, 4); 
+
+indApo = find(tGps > main.tApogee, 1, 'first');
+
+figure
+geoplot(lat, lon, 'LineWidth', 2); 
+geobasemap satellite
+hold on
+geoplot(lat(1), lon(1), 'Marker','o', 'MarkerSize', 10, 'LineWidth',2); 
+geoplot(lat(indApo), lon(indApo),  'Marker','square', 'MarkerSize', 10, 'LineWidth',2); 
+geoplot(lat(end), lon(end), 'Marker','x', 'MarkerSize', 10, 'LineWidth',2, 'Color', 'r'); 
+legend('Trajectory', 'Liftoff', 'Apogee', 'Landing'); 
+
+
+%% 3D TRAJECTORY - GPS
+altNas = -main.NASData(:, 4); 
+tNas = main.NASData(:, 1); 
+altNasInterp = interp1(tNas, altNas, tGps); 
+
+uif = uifigure;
+g = geoglobe(uif);
+geoplot3(g,lat,lon,altNasInterp + alt0,'b','Linewidth',2.5);
\ No newline at end of file