%It provides an example to numerically solve Kepler's equation

clear all
clc
close all

%constants of the problem START%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
G = 6.6742*10^(-20); %km^3/(kg*s^2) gravitational constant
m1 = 5.972*10^24; %kg mass of Earth
m2 = 1000; %kg mass of S/C
mu = G*(m1+m2); %km^3/s^2 gravitational parameter
%constants of the problem END%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r0 = [8000; 0; 6000]; %km 
v0 = [0; 10; 0]; %km/s

h = cross(r0, v0); %km^2/s angular momentum vector
e = cross(v0, h)/mu - r0/norm(r0); %[-] eccentricity vector
enorm = norm(e); %printing norm of e vector
rp = (norm(h)^2/mu)*(1/(1+norm(e))); %km perigee
ra = (norm(h)^2/mu)*(1/(1-norm(e))); %km apogee
a = (rp+abs(ra))/2; %km semi-major axis
theta_inf = acos(-1/enorm)*180/pi

%numerically solving Kepler's equation for hyperbola
tdesired = 100000000000; %s time from perigee where we want true anomaly value
Mh = ((mu^2/norm(h)^3)*(enorm^2-1)^(3/2))*tdesired;
X0 = 0; %radians, initial guess for eccentric anomaly
eccentric_anomaly = fzero(@(F) keplerHyper(F, Mh, enorm), X0)
eccentric_anomaly_deg = (eccentric_anomaly/pi)*180
trueanomaly = 2*atan((sqrt((1+enorm)/(enorm-1)))*tanh(eccentric_anomaly/2));
trueanomaly_deg = (trueanomaly/pi)*180

