function f = satelliteForce( m, r, v )

% Force on a satellite including gravity and atmospheric drag
% Intup:
%   m       : satellite mass (kg)
%   r(3,1)  : satellite position (m)
%   v(3,1)  : satellite velocity (m/s)
% Output:
%   f(3,1)  : force on satellite (N)
% J.M.Soler, Apr.2018

M = 5.972e24;  % Earth mass (kg)
R = 6.371e6;   % Earth radius (m)
h = 6.0e5;     % atmosphere thickness (m)
G = 6.672e-11; % gravitational constant (N*m^2/kg^2)
mur = 0.15;    % drag coeff (kg/s)

rmod = sqrt(sum(r.^2));
rmod = repmat(rmod,3,1);
mu = (1-(rmod-R)/h)*mur;
mu = max(mu,0);
f = -G*M*m*r./rmod.^3 - mu.*v;

end

