asympt (cbmu, rsc, vsc)

function [shat, iberr] = asympt (cbmu, rsc, vsc)

% c3-scaled outgoing unit asymptote vector

% input

%  cbmu  = central body gravitational constant (km**3/sec**2)
%  rsc   = spacecraft position vector (kilometers)
%  vsc   = spacecraft velocity vector (kilometers/second)
%  iberr = error ==> orbit is not hyperbolic (c3 < 0)

% output

%  shat = c3-scaled outgoing unit asymptote vector

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

zmu = 1.0d0 / cbmu;

iberr = 0;

% position magnitude (kilometers)

rmag = norm(rsc);

% velocity magnitude (kilometers/second)

vmag = norm(vsc);

% twice specific orbital energy (km/sec)**2

c33 = vmag * vmag - 2.0d0 * cbmu / rmag;

if (c33 < 0.0d0)
    
   shat = 0.0;
   
   iberr = 1;
   
   return
   
end

% angular momentum vector

hv = cross(rsc, vsc);

hmag = norm(hv);

vxh = cross(vsc, hv);

% compute unit position and eccentricity vectors

for i = 1: 3
    
    rhat(i) = rsc(i) / rmag;

    ev(i) = zmu * vxh(i) - rhat(i);
    
end

fac1 = 1.0d0 / (1.0d0 + c33 * hmag * hmag / cbmu / cbmu);

fac2 = sqrt(c33) / cbmu;

hxe = cross(hv, ev);

% compute components of c3-scaled unit asymptote vector

shat(1) = fac1 * (fac2 * hxe(1) - ev(1));

shat(2) = fac1 * (fac2 * hxe(2) - ev(2));

shat(3) = fac1 * (fac2 * hxe(3) - ev(3));