%%
%% 3.2 Action field of the Kepler/Coulomb problem (Computational example)
% We consider the dynamics of the Kepler problem
% Energy: E={\vec p \over 2m}-{G_m\over r}
% Parameter: G_m=1,\ m=1,\ p_\varphi=1,\ E=-{1\over 5}
%

clear all
syms m Gm p_phi Xma r phi positive
syms E negative
Par=[E==-1/5 m==1 Gm==1 p_phi==1 Xma==4]

%
%% 1 Coordinates and metric for the x-y-plane
% Polar-coordinates: $$\vec r = \left( {\matrix{
%    {r\cos \varphi }  \cr
%    {r\sin \varphi }  \cr
%
%  } } \right)$$

vr=[r*cos(phi);r*sin(phi)]

% 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(vr,[r,phi])
G=simplify(J'*J)
G_i=G^(-1)

%
%% 2 Energy
% Momentum: $$\vec p = \left( {\matrix{
%    {{p_r}}  \cr
%    {{p_\varphi }}  \cr
%
%  } } \right)$$,    Energy: E={1\over 2m}
% {\vec p}\prime G_i\vec p

syms E p_r p_phi
p_v=[p_r;p_phi]
e2=E==(transpose(p_v)*G_i*p_v)/2/m-Gm/r

%
% Separation of the variables
% Everything that depends on \varphi goes to the left, everything that depends on r to the right

solve(e2,p_phi,'IgnoreProperties', true')
const=p_phi==ans(1)

%
%
%% 3 Effective potential
% V_{eff}(r)=E|_{p_r=0}

V_eff=sube(E,[e2 p_r==0])
V_effp=sube(V_eff,Par)

%
%% 4 Perihelion and aphelion
% p_r(r)=0\ \to \ r_{ap},\ r_{pe}
% Limits at p_r=0

r_ap=solve(subs(e2,p_r,0),r,'IgnoreProperties', true')
r_ap=[0 1;1 0]*sube(r_ap,Par)
vpa(r_ap,3)

%% 5 Graphics

fplot(sube(V_eff,Par),'b')
hold on

sube(E,Par)
plot([0,10],[ans ans],'r')
plot(r_ap(1)*[1 1],0.5*[-1 1],'g')
plot(r_ap(2)*[1 1],0.5*[-1 1],'g')
plot([0 10],[0 0],':k')
axis([0 10 -.5 .5])
xlabel('Distance r'),ylabel('V_eff')
legend('Veff','energy','aphel','perihel')

%
%% 6 Momentum field
% \[\left( {\begin{array}{*{20}{c}}
%   {{p_x}} \\
%   {{p_y}}
% \end{array}} \right) = \vec p = {\vec {g^i}}{p_i} = {\vec {g_j}}{g^{ji}}{p_i} = J{G^{ - 1}}\left( {\begin{array}{*{20}{c}}
%   {{p_r}} \\
%   {{p_\varphi }}
% \end{array}} \right)\]

syms x y

% Radial momentum

e7=p_r==solve(e2,p_r)

% Transformation into Cartesian coordinates

double(sube(Xma,Par))
[X Y]=meshgrid(linspace(-ans,ans,12));
xyP=[sin(phi)==y/r cos(phi)==x/r r==sqrt(x^2+y^2)]
p_v7=J*G_i*p_v

% inwards

figure
p_v7p=sube(p_v7,[xyP e7(2) xyP Par]);
VectorField(p_v7p,1,X,Y,'incoming')

% to the outside

p_v7p=sube(p_v7,[xyP e7(1) xyP Par]);
VectorField(p_v7p,2,X,Y,'outgoing')

%
%% 7 Action field radial
% S_r=\int p_rdr-\int p_rdr|_{r=r_{pe}} (includes gauge)

pr=sube(p_r,[e7(1),Par]);

S_rf1=GetS_r(e7,r_ap,Par,2);
S_rf2=GetS_r(e7,r_ap,Par,1);
vpa(S_rf1(r_ap(2)),5)

Plo_S_r(r_ap,S_rf1,1,'incoming')
Plo_S_r(r_ap,S_rf2,2,'outgoing')

%
%% 8a Action field S(r,phi)
% S_{r,\varphi}=S_r+S_\varphi\  {\rm mit}\
% S_\varphi=\int p_\varphi d\varphi

S_p=p_phi*phi
S_pf=matlabFunction(sube(S_p,Par));

% Graphics of S

Plo_S_rphi(S_rf1,S_pf,r_ap,1,'incoming')
Plo_S_rphi(S_rf2,S_pf,r_ap,2,'outgoing')

%
%% 8b Action field -  Plot symbolic

Sp=sube(p_phi*phi,Par);
figure
txt={'incoming' 'outgoing'};
for n=1:2
    pr=sube(p_r,[e7(3-n),Par]);
    int(pr,r);
    Sr=ans-sube(ans,r==r_ap(1));
    
    subplot(1,2,n)
    fsurf(r*cos(phi),r*sin(phi),real(Sr)+Sp,double([0 2*pi r_ap']),'MeshDensity',15)
    title(['S(r,phi)  ' txt{n}])
    xlabel('Position x')
    ylabel('Position y')
    zlabel('S_r(r)')
    axis square
end

%
%% 9 Animation of S(r,phi,t) \rightarrow'Kepler_S.gif'
% S=S_r+S_\varphi-Et

N=15;
[rr,phi_p]=meshgrid(double(linspace(r_ap(1),r_ap(2),N)),linspace(0,10*pi,5*N));
f=@(t)surf(rr.*cos(phi_p),rr.*sin(phi_p),...
    real(S_rf1(rr))+S_pf(phi_p)-6*pi*t*ones(size(rr)))
figure
fanimator(f,'AnimationRange',[0 1],'FrameRate',75)
title('S(r,phi,t)')
xlabel('Position x')
ylabel('Position y')
zlabel('S(r,phi,t)')
axis equal
axis([-4 4 -4 4 0 2*pi])

%

playAnimation

%

writeAnimation('Kepler_S.gif','Loopcount',inf)

% _________________________________________________________________
%
%% Functions
% Plot Vector field for momenta

function VectorField(p_vP,n,X,Y,txt)
u=matlabFunction(p_vP(1));
v=matlabFunction(p_vP(2));
subplot(1,2,n)
quiver(X,Y,real(u(X,Y)),real(v(X,Y)),'b')
axis square
xlabel('x-direction')
ylabel('y-direction')
title(['Momentum field ' txt])
end

% _________________________________________________________________
% MatLab function for actionS_r

function S_rf=GetS_r(e7,r_ap,Par,n)
syms p_r r positive
pr=sube(p_r,[e7(n),Par]);
int(pr);
S_r=simplify(ans-sube(ans,r==r_ap(1)));
S_rf=matlabFunction(S_r);
end

% _________________________________________________________________
% Plot S(r)

function Plo_S_r(r_ap,S_rf,n,txt)
rr=double(linspace(r_ap(1),r_ap(2),10));
subplot(1,2,n)
plot(rr,real(S_rf(rr)))
title(['S_r(r)  ' txt])
xlabel('Position r')
ylabel('S_r(r)')
axis square
end

% _________________________________________________________________
% 3D-Plot S(r,phi)

function Plo_S_rphi(S_rf,S_phi,r_ap,n,txt)
N=15;
[rr,phi_p]=meshgrid(double(linspace(r_ap(1),r_ap(2),N)),linspace(0,2*pi,N));
subplot(1,2,n)
surfc(rr.*cos(phi_p),rr.*sin(phi_p),real(S_rf(rr))+S_phi(phi_p))
title(['S(r,phi)  ' txt])
xlabel('Position x')
ylabel('Position y')
zlabel('S_r(r)')
axis square
end