% reentry.m
% Simulate a satellite reentry
% JMS, Jan.2020

% Set parameters
clear
m = 5.5e3;       % satellite mass (kg)
A = 12.0;        % satellite cross sectional area (m^2)
Cd = 0.4;        % satellite drag coefficient
h0 = 150e3;      % initial height (m)
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)
T0 = 247;        % atmospheric temperature (K)
DT = 1;          % uncertainty in atmospheric temperature (K)
dt = 10;         % integration time interval (s)
tmax = 2e5;      % max simulation time (s)

% Find initial position and velocity
r0 = R+h0;          % initial distance to Earth centre
v0 = sqrt(G*M/r0);  % inital orbital velocity (norm)
r0 = [r0;0;0];      % initial position vector
v0 = [0;v0;0];      % initial velocity vector

% Loop on temperature
for iT = -1:+1
  T = T0+iT*DT;     % atmospheric temperature
  
  % Find trajectory
  myforce = @(r,v,t)reentryForce(r,v,m,A,Cd,T);
  r = rungeKutta( myforce, m, r0, v0, dt, tmax );
  r = squeeze(r);
  
  % Remove trajectory below surface
  nt = find(sum(r.^2)<R^2,1);   % first position below surface
  r = r(:,1:nt);

  % Find time, angle, and distance travelled
  t = (0:nt-1)*dt;              % time
  h = sqrt(sum(r.^2))-R;        % height above surface
  phi = atan2(r(2,:),r(1,:));   % angle in inertial system
  it = 1+find(diff(phi)<0);     % times when phi decreases
  for j=1:numel(it)             % add 2*pi every time that phi decreases
    phi(it(j):nt) = phi(it(j):nt)+2*pi;
  end
  phi = phi-w*t;                % angle relative to Earth surface
  fprintf('T = %3.0f K. Fall after %5.2f h, %5.2f orbits, %7.0f km\n', ...
          T, t(nt)/3600,phi(nt)/(2*pi),phi(nt)*R/1e3)

  % Plot height as a function of orbit number
  figure(1)
  plot(phi/2/pi,h/1e3), hold on
  grid on
  xlabel('orbit number')
  ylabel('height (km)')

end % for iT
