From c964fbaa43de88fdd322ff25ae89f412027ce43a Mon Sep 17 00:00:00 2001
From: davideperani <91607731+davideperani@users.noreply.github.com>
Date: Mon, 25 Sep 2023 18:38:32 +0200
Subject: [PATCH] Added functions to estimate parachute opening force

---
 .../postprocessing/RCS/descentPostProcess.m   |   6 +-
 ...rceReductionFactorFromBallisticParameter.m | 231 ++++++++++++++++++
 ...uteForceReductionFactorFromCanopyLoading.m |  77 ++++++
 3 files changed, 313 insertions(+), 1 deletion(-)
 create mode 100644 RoccarasoFlight/postprocessing/commonFunctions/ParachuteOpeningForceModels/computeForceReductionFactorFromBallisticParameter.m
 create mode 100644 RoccarasoFlight/postprocessing/commonFunctions/ParachuteOpeningForceModels/computeForceReductionFactorFromCanopyLoading.m

diff --git a/RoccarasoFlight/postprocessing/RCS/descentPostProcess.m b/RoccarasoFlight/postprocessing/RCS/descentPostProcess.m
index a29c393..ecd2961 100644
--- a/RoccarasoFlight/postprocessing/RCS/descentPostProcess.m
+++ b/RoccarasoFlight/postprocessing/RCS/descentPostProcess.m
@@ -63,7 +63,11 @@ A_dcm = (qs.^2 - pagemtimes(qq, 'transpose', qq, 'none')) .* eye(3) + ...
 
 
 % A_dcm = repmat(eye(3), [1 1 10]);
-animateOrientation(A_dcm, time_qq);
+% animateOrientation(A_dcm, time_qq);
+
+%% Opening force estimation
+
+
 
 function para = descentPost(main, desc, h0, name, S, M)
 
diff --git a/RoccarasoFlight/postprocessing/commonFunctions/ParachuteOpeningForceModels/computeForceReductionFactorFromBallisticParameter.m b/RoccarasoFlight/postprocessing/commonFunctions/ParachuteOpeningForceModels/computeForceReductionFactorFromBallisticParameter.m
new file mode 100644
index 0000000..24c1493
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/commonFunctions/ParachuteOpeningForceModels/computeForceReductionFactorFromBallisticParameter.m
@@ -0,0 +1,231 @@
+function FRF = computeForceReductionFactorFromBallisticParameter(A, n)
+
+    %%%%%%%%%%%%%% INFO %%%%%%%%%%%%%
+    
+    % Function to compute the Force Reduction Factor from the ballistic
+    % parameter, equivalent to using the graph on fig. 5.51 of La Bibbia
+    %
+    % INPUT:
+    % A [1] Ballistic Parameter, must be between 0.01 and 1000, otherwise
+    % the function returns an error
+    % n [1] exponent of drag area vs time. Can be either 1, 2 or 1/2 (0.5) 
+    % otherwise the function returns an error
+    %
+    % OUTPUT:
+    % FRF [1] = Force Reduction Factor
+    %
+    % CONTRIBUTORS:
+    % Tommaso Mauriello
+    %
+    % VERSION:
+    % 9 nov 2020, Fourth version
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    A_min = 0.01;
+    A_max = 1000;
+    
+    errorMessageOutOfBounds = ['The input Ballistic Parameter A (%d) is out of the graph',...
+												'bounds. It must be between 0.01 and 1000.'];
+    errorMessageWrongN = 'The requested n index (%d) is not valid. It must be 1, 2 or 0.5.';
+    
+    if A < A_min || A > A_max
+        % if the input Ballistic Parameter A is out of the graph bounds
+        % then return error
+        error(errorMessageOutOfBounds, A);
+    else
+        switch n
+            case 1
+                % if evalin('caller', 'exist(''FRF_A_n_1_curve'',''var'')')
+                %     % if the FRF_A_curve has already been created it's saved in
+                %     % the workspace. We use it instead of computing it once again
+                %     savedCurve = evalin('caller', 'FRF_A_n_1_curve');
+                % 
+                %     % we return the value of FRF
+                %     FRF = savedCurve(A);
+                % else
+                    % it's the first time that the FRF_PCF_curve is needed. We 
+                    % get it by interpolating the data obtained from the graph 
+                    % and then we save it in the workspace for further usage
+
+                    % data obtained from La Bibbia fig. 5.51, curve for n = 1
+                    A_data_n_1 = [0.01, 0.011270759, 0.012329285, 0.013357006, 0.014192211, ...
+						0.015189805, 0.016336641, 0.017274183, 0.018309839, 0.019407636, ...
+						0.020923608, 0.022285936, 0.023622127, 0.0250384, 0.02719134, ...
+						0.029601226, 0.033581491, 0.037729032, 0.040973288, 0.045148972, ...
+						0.047855955, 0.050848421, 0.054290733, 0.058106995, 0.062342461, ...
+						0.067867424, 0.072111222, 0.077180042, 0.082405255, 0.088626656, ...
+						0.096247268, 0.101032694, 0.109189228, 0.11714812, 0.123871122, ...
+						0.13225672, 0.140527185, 0.152981682, 0.162547111, 0.173973269, ...
+						0.191237662, 0.208692039, 0.225539798, 0.239643573, 0.257735764, ...
+						0.277194555, 0.298846433, 0.324543738, 0.347356469, 0.366400631, ...
+						0.387427465, 0.406689369, 0.436334405, 0.467003898, 0.502262865, ...
+						0.534965691, 0.576753004, 0.615797749, 0.647984231, 0.688504902, ...
+						0.766069417, 0.817932502, 0.866971916, 0.941521382, 1.015064333, ...
+						1.075924225, 1.154349222, 1.238492263, 1.32233511, 1.432554313, ...
+						1.590063152, 1.689484446, 1.852635856, 2.021711541, 2.233124572,...
+						2.401701751, 2.558062553, 2.737852125, 2.951681429, 3.159131874,...
+						3.405859244, 3.636394917, 3.882530119, 4.145330673, 4.404490468, ...
+						4.748461814, 5.020913071, 5.334819769, 5.751453094, 6.155661708, ...
+						6.524659739, 7.017166357, 7.528576382, 8.077226926, 8.58220292, ...
+						9.388057159, 10, 40, 1000]';
+                    FRF_data_n_1 = [0.044144116, 0.046906657, 0.049002192, 0.05119117, ...
+						0.052704682, 0.054527041, 0.056824631, 0.058646726, 0.060234291, ...
+						0.062165769, 0.064785286, 0.066700777, 0.06883961, 0.070874854, ...
+						0.073861539, 0.077536487, 0.082789711, 0.087970619, 0.092123686, ...
+						0.097178101, 0.100294224, 0.103510357, 0.106829803, 0.110792127, ...
+						0.114623062, 0.11945341, 0.12328392, 0.127546683, 0.132598623, ...
+						0.137516949, 0.143311949, 0.147190898, 0.153765934, 0.159082803, ...
+						0.164183828, 0.169038355, 0.175307551, 0.184029901, 0.188553697, ...
+						0.196022204, 0.204283502, 0.21444772, 0.223484223, 0.231772675, ...
+						0.23978704, 0.249285296, 0.25916001, 0.271394878, 0.280778843, ...
+						0.289079578, 0.297625962, 0.304201871, 0.316251407, 0.325602491, ...
+						0.33932231, 0.350203649, 0.364960391, 0.37666421, 0.385921999, ...
+						0.400234849, 0.423221013, 0.438917934, 0.451894544, 0.473228353, ...
+						0.491973921, 0.507749605, 0.525306407, 0.544790521, 0.562261264, ...
+						0.583118056, 0.604753671, 0.619619637, 0.64104956, 0.66321952, ...
+						0.679531947, 0.694551811, 0.708180878, 0.722078615, 0.738039527, ...
+						0.750699572, 0.765433677, 0.776676208, 0.786174035, 0.797721197, ...
+						0.805518179, 0.815371487, 0.823339581, 0.833406626, 0.845650412, ...
+						0.85599244, 0.864358233, 0.874930543, 0.887782854, 0.894290683, ...
+						0.903031541, 0.911867156, 0.918549968, 0.98, 1]';
+
+                    % we interpolate using 'pchipinterp' (Piecewise Cubic Hermite 7
+                    % Interpolating Polynomial) since it looks better than 'smoothingspline'
+                    interpolatedCurve = fit(A_data_n_1, FRF_data_n_1, 'pchipinterp');
+
+                    % we save the curve object in the workspace for further usage
+                    % assignin('caller', 'FRF_A_n_1_curve', interpolatedCurve);
+
+                    % we return the value of FRF
+                    FRF = interpolatedCurve(A);
+                % end
+            case 2
+                % if evalin('caller', 'exist(''FRF_A_n_2_curve'',''var'')')
+                %     % if the FRF_A_curve has already been created it's saved in
+                %     % the workspace. We use it instead of computing it once again
+                %     savedCurve = evalin('caller', 'FRF_A_n_2_curve');
+                % 
+                %     % we return the value of FRF
+                %     FRF = savedCurve(A);
+                % else
+                    % it's the first time that the FRF_PCF_curve is needed. We 
+                    % get it by interpolating the data obtained from the graph 
+                    % and then we save it in the workspace for further usage
+
+                    % data obtained from La Bibbia fig. 5.51, curve for n = 2
+                    A_data_n_2 = [0.01, 0.011130674, 0.011898702, 0.012675246, 0.01352295, ...
+						0.014317013, 0.016825234, 0.017919619, 0.018925181, 0.020431782, ...
+						0.021758606, 0.023189558, 0.024864972, 0.026661561, 0.028701934, ...
+						0.030870314, 0.032968083, 0.035890116, 0.038134752, 0.040274559, ...
+						0.042325435, 0.045523046, 0.048176292, 0.05139804, 0.055043372, ...
+						0.059379339, 0.063284615, 0.068780331, 0.073871615, 0.080663463, ...
+						0.087017618, 0.092041134, 0.098242012, 0.104702883, 0.113294043, ...
+						0.125982261, 0.134675807, 0.142233154, 0.151587339, 0.161556716, ...
+						0.17428285, 0.185744833, 0.207804985, 0.223496266, 0.239156198, ...
+						0.257738975, 0.273027679, 0.289223286, 0.304866904, 0.322592066, ...
+						0.345899465, 0.368647493, 0.397687202, 0.429014475, 0.457227089, ...
+						0.487295783, 0.509975841, 0.545163809, 0.581014571, 0.631527836, ...
+						0.690658392, 0.751872947, 0.818510465, 0.885660179, 0.94964643, ...
+						1.018250614, 1.095124364, 1.181379986, 1.243880326, 1.308129232, ...
+						1.395597368, 1.476918875, 1.575666424, 1.66748056, 1.784415453, ...
+						1.901742972, 2.008430599, 2.146990295, 2.27431396, 2.458716584, ...
+						2.623107629, 2.881136272, 3.05199236, 3.232985657, 3.456004146, ...
+						3.660944989, 3.961246992, 4.296621994, 4.540330935, 4.809572293, ...
+						5.125758482, 5.4627398, 5.839562325, 6.284157152, 6.774908483, ...
+						7.220309412, 7.671685259, 8.102001474, 9.118993962, 10, 40, 1000]';
+                    FRF_data_n_2 = [0.021215028, 0.022817684, 0.023880278, 0.025083749, ...
+						0.026332378, 0.027624349, 0.031212646, 0.032840164, 0.034369347, ...
+						0.036109668, 0.037874089, 0.039637801, 0.041609741, 0.044079027, ...
+						0.04678896, 0.049517293, 0.052084568, 0.055175765, 0.058096452, ...
+						0.060433984, 0.062622328, 0.066273927, 0.06957274, 0.072447286, ...
+						0.076108691, 0.080381293, 0.084636305, 0.090149118, 0.095268682, ...
+						0.101536986, 0.107237085, 0.112176231, 0.117813746, 0.123300069, ...
+						0.131014351, 0.141339797, 0.149273876, 0.156224728, 0.163996398, ...
+						0.172154683, 0.181268524, 0.190286041, 0.20841997, 0.220788667, ...
+						0.23246786, 0.24552724, 0.255406391, 0.265683044, 0.275554249, ...
+						0.288366128, 0.303631538, 0.317770981, 0.335610059, 0.354450589, ...
+						0.368713217, 0.384714792, 0.395366992, 0.412525739, 0.429125234, ...
+						0.451554776, 0.475767101, 0.504003677, 0.530687246, 0.550372581, ...
+						0.576003431, 0.597367818, 0.617649168, 0.64055959, 0.65431614, ...
+						0.666075733, 0.682460257, 0.69924583, 0.710671683, 0.728151136, ...
+						0.741048948, 0.756966815, 0.770882564, 0.789833968, 0.799490486, ...
+						0.809041718, 0.822261664, 0.834237756, 0.841879932, 0.852172763, ...
+						0.862594191, 0.867859994, 0.883821884, 0.891726323, 0.897369788, ...
+						0.902847883, 0.905610548, 0.91114089, 0.919490488, 0.928672786, ...
+						0.936424915, 0.942143432, 0.945025333, 0.947914031, 0.959523673, ...
+						0.968325926, 0.99, 1]';
+
+                    % we interpolate using 'pchipinterp' (Piecewise Cubic Hermite 7
+                    % Interpolating Polynomial) since it looks better than 'smoothingspline'
+                    interpolatedCurve = fit(A_data_n_2, FRF_data_n_2, 'pchipinterp');
+
+                    % we save the curve object in the workspace for further usage
+                    % assignin('caller', 'FRF_A_n_2_curve', interpolatedCurve);
+
+                    % we return the value of FRF
+                    FRF = interpolatedCurve(A);
+                % end
+            case 0.5
+                % if evalin('caller', 'exist(''FRF_A_n_05_curve'',''var'')')
+                %     % if the FRF_A_curve has already been created it's saved in
+                %     % the workspace. We use it instead of computing it once again
+                %     savedCurve = evalin('caller', 'FRF_A_n_05_curve');
+                % 
+                %     % we return the value of FRF
+                %     FRF = savedCurve(A);
+                % else
+                %     % it's the first time that the FRF_PCF_curve is needed. We 
+                %     % get it by interpolating the data obtained from the graph 
+                %     % and then we save it in the workspace for further usage
+                % 
+                %     % data obtained from La Bibbia fig. 5.51, curve for n = 0.5
+                    A_data_n_05 = [0.01, 0.010951767, 0.011807143, 0.012667726, 0.013642837, ...
+						0.014555066, 0.015418268, 0.016502002, 0.017533817, 0.018949181, ...
+						0.020445024, 0.023860471, 0.026742116, 0.03093212, 0.033673281, ...
+						0.036607438, 0.042710446, 0.049044166, 0.054834313, 0.060275537, ...
+						0.067065536, 0.073008511, 0.076438598, 0.081549322, 0.088216597, ...
+						0.096501096, 0.106077197, 0.117170728, 0.129110057, 0.151161957, ...
+						0.164857254, 0.197874321, 0.221771738, 0.25098006, 0.271240009, ...
+						0.293135035, 0.333268233, 0.352687785, 0.390289253, 0.410688345, ...
+						0.444917743, 0.480833604, 0.502293675, 0.531117042, 0.576780275, ...
+						0.623338313, 0.66714956, 0.703723996, 0.762376863, 0.806123682, ...
+						0.877557878, 0.959973789, 1.094717055, 1.176960702, 1.268889132, ...
+						1.36136906, 1.511047244, 1.689436949, 1.92589631, 2.127292535, ...
+						2.293439385, 2.498550763, 2.791413604, 3.220934983, 3.653961014, ...
+						4.236701475, 4.545477448, 5.106789707, 5.599903094, 6.081351691, ...
+						6.636298543, 7.189373432, 7.826425701, 8.458146888, 9.140834995, ...
+						10, 20, 1000]';
+                    FRF_data_n_05 = [0.095955515, 0.09927411, 0.101962065, 0.104722621, ...
+						0.106787857, 0.109414688, 0.111279201, 0.114569561, 0.117385893, ...
+						0.120564347, 0.12355117, 0.13062536, 0.135472456, 0.142561339, ...
+						0.147491532, 0.151286924, 0.159800009, 0.167753932, 0.175249612, ...
+						0.180871495, 0.18895298, 0.19406975, 0.197649803, 0.20087945, ...
+						0.206219738, 0.212834757, 0.220730881, 0.22947645, 0.235691582, ...
+						0.249837661, 0.258223827, 0.275995725, 0.286237058, 0.300482752, ...
+						0.309368634, 0.317745399, 0.329262665, 0.337361104, 0.348469752, ...
+						0.356169335, 0.365813632, 0.377546421, 0.381234312, 0.389658505, ...
+						0.399240095, 0.409056595, 0.418097362, 0.42733576, 0.438907091, ...
+						0.447518543, 0.459637163, 0.477846798, 0.505823651, 0.52023635, ...
+						0.536921502, 0.547458816, 0.565022831, 0.587412078, 0.61069267, ...
+						0.627233185, 0.642655024, 0.66085731, 0.682891018, 0.709958701, ...
+						0.730966068, 0.758098961, 0.771103755, 0.793917673, 0.809504854, ...
+						0.82740036, 0.843643478, 0.858118668, 0.874964086, 0.889975916, ...
+						0.900863101, 0.91855075, 0.98, 1]';
+
+                    % we interpolate using 'pchipinterp' (Piecewise Cubic Hermite 7
+                    % Interpolating Polynomial) since it looks better than 'smoothingspline'
+                    interpolatedCurve = fit(A_data_n_05, FRF_data_n_05, 'pchipinterp');
+
+                    % we save the curve object in the workspace for further usage
+                    % assignin('caller', 'FRF_A_n_05_curve', interpolatedCurve);
+
+                    % we return the value of FRF
+                    FRF = interpolatedCurve(A);
+                % end
+            otherwise
+                % if n is different from 1, 2 and 0.5 then return error
+                error(errorMessageWrongN, n);
+        end
+    end   
+end
diff --git a/RoccarasoFlight/postprocessing/commonFunctions/ParachuteOpeningForceModels/computeForceReductionFactorFromCanopyLoading.m b/RoccarasoFlight/postprocessing/commonFunctions/ParachuteOpeningForceModels/computeForceReductionFactorFromCanopyLoading.m
new file mode 100644
index 0000000..2f7bb74
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/commonFunctions/ParachuteOpeningForceModels/computeForceReductionFactorFromCanopyLoading.m
@@ -0,0 +1,77 @@
+function FRF = computeForceReductionFactorFromCanopyLoading(PCL)
+
+    %%%%%%%%%%%%%% INFO %%%%%%%%%%%%%
+    
+    % Function to compute the Force Reduction Factor from the Canopy
+    % Loading, equivalent to using the graph on fig. 5.48 of La Bibbia
+    
+    % INPUT:
+    % PCL [1] = Parachute Canopy Loading, must be between 0 and 100 
+    % otherwise the function returns an error
+    % 
+    % OUTPUT:
+    % FRF [1] = Force Reduction Factor
+    %
+    % CONTRIBUTORS:
+    % Tommaso Mauriello
+    %
+    % VERSION:
+    % 9 nov 2020, Fourth version
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    PCL_min = 0;
+    PCL_max = 100;
+    
+    errorMessageOutOfBounds = ['The input Parachute Canopy Loading PCL (%d) is out of'...
+											'the graph bounds. It must be between 0 and 100.'];
+
+    if PCL < PCL_min || PCL > PCL_max
+         %if the input Canopy Loading PCF is out of the graph bounds then returns error
+         %formattedError = sprintf(errorMessageOutOfBounds, PCL);
+         error(errorMessageOutOfBounds, PCL);
+    else
+        if evalin('caller', 'exist(''FRF_PCL_curve'', ''var'')')
+            % if the FRF_PCF_curve has already been created it's saved in
+            % the workspace. We use it instead of computing it once again
+            savedCurve = evalin('caller', 'FRF_PCL_curve');
+            FRF = savedCurve(PCL);
+        else
+            % it's the first time that the FRF_PCF_curve is needed. We 
+            % get it by interpolating the data obtained from the graph 
+            % and then we save it in the workspace for further usage
+            
+            %data obtained from La Bibbia fig. 5.48
+            PCL_data = [0, 0.2, 0.315238664, 0.353618869, 0.396643869, 0.469187078, 0.518862451, ...
+				0.55955745, 0.618800824, 0.667294804, 0.719589146, 0.775981673, 0.858088507, ...
+				0.927411952, 0.99163793, 1.075987155, 1.164534389, 1.251414299, 1.349138417, ...
+				1.454781569, 1.568697018, 1.711207738, 1.90463473, 2.01757444, 2.240172577, ...
+				2.518493085, 2.728899831, 3.01197155, 3.258187307, 3.504498865, 3.87060744, ...
+				4.358457656, 4.675440475, 4.965487707, 5.38024354, 6.731032391, 7.428725575, ...
+				8.295048363, 9.370986014, 10.84382219, 12.10145674, 13.25453732, 14.18324569, ...
+				15.39337664, 17.13135833, 18.7805021, 20.86587709, 23.60741993, 26.81084133, ...
+				31.11332399, 33.70945518, 38.14698112, 44.21401955, 56.83565422, 64.70349943, ...
+				69.52742428, 74.70418146, 80.68133451, 87.13535882, 93.48211331, 100]';
+            FRF_data = [0.020092574, 0.020092574, 0.020902911, 0.021456589, 0.026215821, ...
+				0.032543213, 0.039641814, 0.043210055, 0.050308656, 0.057388317, 0.064467978, ...
+				0.071547639, 0.082157661, 0.089471295, 0.097465569, 0.106908065, 0.115810715, ...
+				0.126145192, 0.13516989, 0.145760971, 0.156352053, 0.170305721, 0.185806806, ...
+				0.19533046, 0.212097039, 0.230409394, 0.246824866, 0.262786929, 0.278266572, ...
+				0.290392115, 0.309704495, 0.332690949, 0.34722972, 0.358957219, 0.37630726, ...
+				0.422270079, 0.442437695, 0.466352411, 0.495416078, 0.526509595, 0.550748537, ...
+				0.572309583, 0.587316492, 0.606069645, 0.628579987, 0.649209843, 0.671718923, ...
+				0.699156344, 0.726113031, 0.759058834, 0.775350621, 0.802546413, 0.833063827, ...
+				0.878901694, 0.902476415, 0.915614595, 0.927300202, 0.939512313, 0.952659353, ...
+				0.96206163, 0.968656595]';
+            
+            % we interpolate using 'pchipinterp' (Piecewise Cubic Hermite 7
+            % Interpolating Polynomial) since it looks better than 'smoothingspline'
+            interpolatedCurve = fit(PCL_data, FRF_data, 'pchipinterp');
+            
+            % we save the curve object in the workspace for further usage
+            assignin('caller', 'FRF_PCL_curve', interpolatedCurve);
+            
+            FRF = interpolatedCurve(PCL);
+        end
+    end
+end
\ No newline at end of file
-- 
GitLab