clear all;

p=8;
a=-0.38673;
e=4.656854;
alpha=1/a;
mu=pi/3;
r0=sqrt(2);

%p=2.42;
%a=-5.5;
%e=1.2;
%alpha=1/a;
%mu=pi/3;
%r0=1.1;


n=5;
%
%
F=acosh(  (e+cos(mu)) / (1+e*cos(mu)) )
t=sqrt( (-a)^3 )*(e*sinh(F)-F)

X=input('X?  >'); %X=0.383109;
y=alpha*X^2;
C= (1- cosh(sqrt(-y))) / y
S= (sinh(sqrt(-y))-sqrt(-y)) / (sqrt((-y)^3))

FUN=  (1-r0*alpha) * X^3 * S + r0*X; 
FUN= FUN - t

FUNder=(1-r0*alpha)*X^2*C+r0
FUNderder=(1-r0*alpha)*X*( 1- alpha*X^2*S )

if(FUNder<0)
    sign=-1;
else
    sign=1;
end
term=abs( ( (n-1)*FUNder )^2 - n*(n-1)*FUN*FUNderder);                         

newX= X - (   (n*FUN) / ( FUNder + ( (1*sign)* sqrt(term)  )))

f=1- (newX^2/r0)*C
g=t-(newX^3)*S
r=f^2+f^2+(2*g)^2;
r=sqrt(r)

fdot=(newX/(r*r0))*(  alpha*newX^2*S - 1)
gdot=1- (newX^2/r)*C

FINAL= f*gdot - fdot*g