%%
%
%% 4.2 Keplerian orbits (Computational example)
% A body with mass m moves in the
%
% Gravitational potentialV(r)={\alpha \over r}
% The movement takes place in the x-y plane, so that it can be described with
% Polar coordinates$$\vec r = \left( {\matrix{
%    {r\cos \varphi }  \cr
%    {r\sin \varphi }  \cr
%
%  } } \right)$$
%
%

clear all
syms alpha m p_r p_phi r phi positive
syms E
assume(E<0)
Par=[E==-2/25 m==1 alpha==1 p_phi==1]
r_v=[r*cos(phi);r*sin(phi)]

%
%% 1 Ground vectors and metric
% Jacobian J, metric G and contravariant metric G_i
% J=\left( {
% \partial x_i \over\partial q_j
% }\right), \
% G=J\prime J, \
% G_i=G^{-1}

J=jacobian(r_v,[r,phi])
G=simplify(transpose(J)*J)
G_i=G^(-1)

%
%% 2 Hamiltonian
% H={1\over 2m}(p_r,p_\varphi)G_i
% \left({p_r\atop p_\varphi}\right)
% +V(r)

p_v=[p_r;p_phi]
H=1/(2*m)*(p_v'*G_i*p_v)-alpha/r;
H=expand(H)

%
%% 3 Radial momentum
% E=H(p_r,...)\ \to\ p_r(E,...)

E==H
p_rs=solve(ans,p_r,'IgnoreProperties', true')

%
%% 4 Aphelion and perihelion (reversal points)
% p_r(r)=0\ \to \ r_{ap},\ r_{pe}

r_ape=solve(p_rs(2),r,'IgnoreProperties', true);
[r_ap_,ind]=sort((sube(r_ape,Par))); %Sorting
r_ap=r_ape(ind)
double(sube(r_ap,Par))

%
%% 5 Keplerian orbits
% Azimuth angle \varphi(r)+{\pi\over 2}=
% \int d\varphi=\int{m\dot r\over G^{-1}_{rr}p_r(r)}\dot \varphi dt=
% \int{p_\varphi\over r^2p_r(r)}dr
% Therefore: r(\varphi)={p_\varphi^2\over m\alpha+A\cos(\varphi), \left({A\equiv\sqrt{\alpha^2m^2+2Emp_\varphi}}\right)
% Integrand p_\varphi\over r^2p_r(r)

IG=p_phi/(r^2*p_rs(1))

% Integration and substraction of the constant of integration (integal value in the perihel)

int(IG,r)
phi_int=ans-sube(ans,r==r_ap(1))

% Inverting to the form r(\varphi)={p\over a+\epsilon\cos(\varphi)}  (Kepler-orbits)

cos(phi+pi)==simplify(cos(phi_int),'Steps',100)
er=isolate(ans,r)
erP=r==simplify(sube(rhs(er),Par))

% Graphics

fplot(sube(r,erP),[0 2*pi])
xlabel('\phi')
ylabel('r(\phi)')

%
%% 6 Period
% T=2\int\limits_{t(r_{pe})}^{t(r_{ah})}dt=
% 2\int\limits_{t(r_{pe})}^{t(r_{ah})}{m\dot r\over G^{-1}_{rr}p_r(r)}dt=
% 2\int\limits_{r_{pe}}^{r_{ah}}{m\over p_r(r)}dr

T=2*m*int(1/p_rs(1),r,r_ap(1),r_ap(2))

% with parameters

TP=real(sube(T,Par))
TP=double(TP)

%
%% 7 Numerical solution
% Equation of motion: \frac{d}{d t} \phi = \frac{\partial}{\partial p_\phi}H(t)

syms phi(t)
dphi_dt=diff(H,p_phi)
sube(dphi_dt,[Par,erP])
sube(ans,phi==phi(t))
odeFunction(ans,phi(t))
[tn,phi_s] = ode45(ans,linspace(0,TP,500),0);

%
%% 8 Graphics

subplot(1,2,1)
plot(tn,phi_s,[0 100],2*pi*[1 1],':k')
xlabel('Time t')
ylabel('Angle \phi')
title('Azimuth \phi(t)')
axis square
subplot(1,2,2)
r_s=matlabFunction(sube(r,erP));
rx_s=matlabFunction(r_v(1));
ry_s=matlabFunction(r_v(2));

x=rx_s(phi_s,r_s(phi_s));
y=ry_s(phi_s,r_s(phi_s));
plot(x,y,'r')
hold on
rectangle('Position',0.5*[-1 -1 2 2],'Curvature',[1 1],'FaceColor','y')
k=8;
plot(x(k),y(k),'ob','MarkerSize',5);
xlabel('x-direction')
ylabel('y-direction')
title('Keplerian orbit')
axis equal
axis([-2 13 -3 3])
hold off

%
%% 9 Animation \rightarrow 'Keplerian_Orbit.gif'

figure
hold on
plot(x,y,'g')

rectangle('Position',0.4*[-1 -1 2 2],'Curvature',[1 1],'FaceColor','y','EdgeColor','r');
PM=rectangle('Position',0.25*[-1 -1 2 2]+[x(1) y(1) 0 0],'Curvature',[1 1],'FaceColor','b');
axis equal
axis off
axis([-.7 12.5 -3 3])

clear im
im{1}=frame2im(getframe);
cnt=1;
for n=3:4:size(x,1)
    PM.Position=0.25*[-1 -1 2 2]+[x(n) y(n) 0 0];
    im{cnt}=frame2im(getframe);
    cnt=cnt+1;
end

%

FiNa='Keplerian_Orbit.gif';
n=1;
[A,map]=rgb2ind(im{n},16);
imwrite(A,map,FiNa,'gif','LoopCount',Inf,'DelayTime',0);
for n=2:size(im,2)
    [A,map]=rgb2ind(im{n},16);
    imwrite(A,map,FiNa,'gif','WriteMode','append','DelayTime',0);
end