testing.m

global Thrust O Rbody ScaleH Rho0 Aref Cd Mdot

Aref = 2.7770013648359075e-12
Cd = 0.80000000000000004
Mdot = 161107.22668495434
O = 0.072123427831452785
%O = 0
Rbody = 0.99986023590346207
Rho0 = 2.6471097575924006e+17
ScaleH = 0.0095986582646732359
Thrust = 201870.12975876412

y0 = [ -0.99424331272405853 -0.10713242047435829 -0.0016969345289107589 0.0073169237820196412 -0.071752443605568844 -7.0231902251673552e-07 138261.79742440602 ]
%y0 = [ -0.99424331272405853 -0.10713242047435829 -0.0016969345289107589 -0.000409833614471594 -4.42077934106611e-05 -7.02319022516736e-07 138261.79742440602 ]

opts = odeset('RelTol',1e-12,'AbsTol',1e-12);

[ t y ] = ode45(@dydt, [ 0 0.20307864529123865 ], y0, opts);

y(end,:);

[ -1.0047976508300847 -0.12352200260296523 -0.0009132257571778606 -0.12586565020619353 -0.094212823689453398 0.015014678306399685 105544.36008259695];

rf = y(end,1:3)';
vf = y(end,4:6)';

r = y(:,1:3)';
v = y(:,4:6)';

east = normc(cross(repmat([0 0 1]',1,size(r,2)), r));
north = normc(cross(r, east));
up = normc(r);

vcoriolis = cross(repmat([0 0 O]',1,size(r,2)), r);
vsurf = v - vcoriolis;
norm_vsurf = sqrt(sum(vsurf.^2,1));

vert = dot(vsurf,up);
horiz = vsurf - vert.*up;

heading = atan2d(dot(horiz,east), dot(horiz,north))
pitch = asind(vert./norm_vsurf)

function M = normc(M)
   n = sqrt(sum(M.^2,1));
   M = bsxfun(@rdivide,M,n);
end

function dy = dydt(t, y)
  global Thrust O Rbody ScaleH Rho0 Aref Cd Mdot

  r = y(1:3);
  v = y(4:6);
  m = y(7);

  a = Thrust / m;
  rm2 = dot(r,r);
  rm = sqrt(rm2);
  rm3 = rm2 * rm;

  vsurf = v - cross([0 0 O]', r);

  % point along the surface velocity vector
  u = vsurf/norm(vsurf);

  % height in ScaleH units
  h = (rm - Rbody) / ScaleH;
  rho = Rho0 * exp(-h);

  % drag accel
  d = -rho * Aref * Cd * norm(vsurf) * vsurf * 0.5 / m;

  dy = [
    v
    -r/rm3 + a*u + d;
    -Mdot;
  ];
end