diff --git a/RoccarasoFlight/postprocessing/RCS/descentPostProcess.m b/RoccarasoFlight/postprocessing/RCS/descentPostProcess.m
index f36a7dd5907b447381a45e90822e3fc2c244f213..2fe4fa7f0f3392f00039ab40bfa6c67ff6c68ba5 100644
--- a/RoccarasoFlight/postprocessing/RCS/descentPostProcess.m
+++ b/RoccarasoFlight/postprocessing/RCS/descentPostProcess.m
@@ -1,6 +1,4 @@
-clear; 
-clc; 
-close all; 
+clear; clc; close all; 
 
 %% adding path and recall data
 addpath("..\commonFunctions\"); 
@@ -55,6 +53,11 @@ para.descent.altitudetime = [para.descent.pressureRCSalt(:,1); para.descent.stat
 para.descent.altitude = [para.descent.pressureRCSalt(:,2); ...
      para.descent.static1alt(:,2); para.descent.static2alt(:,2); -para.descent.NASalt(:,2) + h0];
 
+para.descent.NASvd = main.NASData(main.NASData(:,1)>desc(1) & main.NASData(:,1)<desc(2),[1,7]);
+
+[para.descent.altitudetime, idx_sort] = sort(para.descent.altitudetime);
+para.descent.altitude = para.descent.altitude(idx_sort,:);
+
 [~, ~, ~, para.descent.rhoVec] = atmosisa(para.descent.altitude);
 
 
@@ -68,19 +71,46 @@ plot(para.descent.NASalt(:,1), -para.descent.NASalt(:,2)+h0)
 legend('RCS analogic', 'STATIC 1', 'STATIC 2', 'NAS')
 grid on
 
-para.descent.cfit = fit(para.descent.altitudetime, para.descent.altitude, 'poly1');
+% Model for speed vs altitude
+
+para.descent.cfit = fit(para.descent.altitudetime, para.descent.altitude, 'poly2');
+para.descent.cfit2 = fit(para.descent.altitudetime, para.descent.altitude, 'poly1');
+
 xlim(desc)
 plot(para.descent.cfit);
 
+para.descent.v_interp = differentiate(para.descent.cfit, para.descent.altitudetime);
+figure
+plot(para.descent.altitudetime, para.descent.v_interp);
+hold on
+plot(para.descent.NASvd(:,1), -para.descent.NASvd(:,2))
+xlabel('Time [s]')
+ylabel('Vz [m/2]')
 %% Descent Cd evaluation
 para.descent.refMass  = M;
 para.descent.paraS    = S;
 para.descent.avgSpeed = -para.descent.cfit.p1;
-para.descent.rho      = [min(para.descent.rhoVec), max(para.descent.rhoVec)];
 para.descent.Cd       = 2 * para.descent.refMass * 9.81 ./ ...
-                       (para.descent.paraS * para.descent.avgSpeed^2 *para.descent.rho);
+                       (para.descent.paraS * para.descent.v_interp.^2 .*rhofun(para.descent.cfit(para.descent.altitudetime)));
 
+para.descent.Cd2       = 2 * para.descent.refMass * 9.81 ./ ...
+                       (para.descent.paraS * para.descent.cfit2.p1.^2 .*rhofun(para.descent.cfit(para.descent.altitudetime)));
+
+figure
+plot(para.descent.altitudetime, para.descent.Cd )
 xlabel('Time [s]')
-ylabel('Altitude agl [m]')
+ylabel('Cd [-]')
+
+hold on
+plot(para.descent.altitudetime, para.descent.Cd2 )
+legend('poly2', 'poly1')
+fprintf('%s Cd = %3.2f~%3.2f\n', name, para.descent.Cd)
+fprintf('%s Vz = %3.2f m/s\n', name, para.descent.avgSpeed)
+
+end
+
+function rho = rhofun(h)
+
+[~, ~, ~, rho] = atmosisa(h);
 
 end
\ No newline at end of file