function f = reentryForce( r, v, m, A, Cd, T )
% Total force on a satellite re-entrying the atmosphere
% Input:
%   r(3,1) : position relative to Earth centre (m)
%   v(3,1) : velocity (m/s)
%   m      : satellite mass (kg)
%   A      : satellite cross sectional area (m^2)
%   Cd     : satellite drag coefficient
%   T      : atmosphere temperature (K)
% Output:
%   f(3,1) : total force, including gravitation and drag (N)
% JMS, Jan.2020

% Set parameters
G = 6.674e-11;  % gravitational constant (N*m^2/kg^2)
M = 5.972e24;   % Earth mass (kg)
R = 6.371e6;    % Earth radius (m)
w = 7.292e-5;   % Earth rotation angular velocity (rad/s)
P0 = 1.01e5;    % atmspheric pressure at surface (N/m^2)
ma = 28.97e-3;  % atmosphere molar mass (kg/mol)
Rg = 8.314;     % gas constant (J/K/mol)

% Find atmosphere density
g = G*M/R^2;    % gravitational acceleration (m/s^2)
h = norm(r)-R;  % satellite height over surface (m)
rho = ma*P0/(Rg*T)*exp(-ma*g*h/Rg/T);  % air density

% Find drag force
va = cross([0;0;w],r);          % velocity of Earth's atmosphere
vd = v-va;                      % drag velocity
f = -0.5*rho*Cd*A*norm(vd)*vd;  % drag force

% Add gravitational force
f = f-G*M*m/norm(r)^3*r;

end
