clc;
clear all;
R0 = [8000.2 -6900.5];      %[km]
V0 = [0.450 8.9000];        %[km/s]
% Compute the position and velocity vectors when the change in satellite
% true anomaly equals to 100 [deg].
dfi = 100;                  %[deg]
% Magnitude of R0 and V0 vectors
r0 = norm(R0);              %[km]
v0 = norm(V0);              %[km/s]
% Radial component of V0
vr0 = R0*V0'/r0;             %[km/s]
% Magnitude of the constant angular momentum
h = r0*(v0^2-vr0^2)^0.5;    %[km^2/s]
mu      = 398600;           % Earth’s gravitational parameter [km^3/s^2]
r = h^2/mu*1/(1+(h^2/mu/r0 -1)*cosd(dfi)- h*vr0/mu*sind(dfi));
% Lagrange coefficients in terms of the change in true anomaly
f = 1 -  mu*r/h^2*(1-cosd(dfi));
g = r*r0/h*sind(dfi);
fd = mu/h*(1 - cosd(dfi))/sind(dfi)*(mu/h^2*(1 - cosd(dfi)) -1/r0 - 1/r);
gd = 1 - mu*r0/h^2*(1 - cosd(dfi));
% Resultant position and velocity vectors
R = f*R0 + g*V0;
V = fd*R0 + gd*V0;
fprintf('R = %4.2f*i + %4.2f*j [km] \n',R(1),R(2));
fprintf('V = %4.4f*i + %4.4f*j [km/s] \n',V(1),V(2));

R = 3634.07*i + 6101.27*j [km] 
V = -7.6623*i + 7.5831*j [km/s]