%
%% 6.5 Schwarzschild metric and black holes (Computational example)
% The equivalence principle combines the continuity of the momentum gradient or the electromagnetic potential with the metric.
% The Schwarzschild metric can be determined by using a suitable approach.
% In consequence, seen from the outside, in the field of a point-like gravitational source, the clocks slow down and the space lines are stretched radially.

clear all
syms r theta phi a0 q mu omega ct e c r_s real
syms A(r) B(r)
Par=[a0==1 q==1 mu==1 omega==1 r_s==1]

%
%% 1 Schwarzschild metric
% Approach: (ds)^2=A(r)(dr)^2+
% r^2(d\vartheta)^2+
% r^2\sin^2\vartheta(d\varphi)^2-B(r)(dct)^2
% Metric:   $G = \left( {\begin{array}{*{20}{c}}
%   {A(r)}&0&0&0 \\
%   0&{{r^2}}&0&0 \\
%   0&0&{{r^2}\sin {{(\varphi )}^2}}&0 \\
%   0&0&0&{ - B(r)}
% \end{array}} \right)$

G=[[A(r) 0 0 0];[0 r^2 0 0];[0 0 r^2*sin(theta)^2 0];[0 0 0 -B(r)]]
G_i=G^-1

%
%% 2 Christoffel symbols
% Definition:  \[{\Gamma _{\mu \nu }^\lambda  \equiv \vec g_{}^{\left( 4 \right)\lambda } \cdot \frac{{\partial \vec g_\mu ^{\left( 4 \right)}}}{{\partial {q^\nu }}}}\]
% or \[{\Gamma _{\mu \nu }^\lambda  = \frac{1}{2}{g^{\lambda \kappa }}\left( {\frac{{\partial {g_{\kappa \mu }}}}{{\partial {q^\nu }}} + \frac{{\partial {g_{\kappa \nu }}}}{{\partial {q^\mu }}} - \frac{{\partial {g_{\mu \nu }}}}{{\partial {q^\kappa }}}} \right)}\]
% or for a diagonal metric: \[\Gamma _{\mu \nu }^\lambda  = \frac{1}{2}{g^{\lambda \lambda }}\left( {\frac{{\partial {g_{\lambda \mu }}}}{{\partial {q^\nu }}} + \frac{{\partial {g_{\lambda \nu }}}}{{\partial {q^\mu }}} - \frac{{\partial {g_{\mu \nu }}}}{{\partial {q^\lambda }}}} \right)\]

Koord=[r,theta,phi,ct];
for lam=1:4
    dG{lam}=diff(G,Koord(lam));
end
Gamm=sym('g',[4,4,4]);
for lam=1:4, for mu=1:4, for nu=1:4
            Gamm(lam,mu,nu)=G_i(lam,lam)*(dG{nu}(lam,mu)+dG{mu}(lam,nu)-dG{lam}(mu,nu))/2;
        end, end, end

% Output of non-disappearing symbols

for lam=1:4, for mu=1:4, for nu=1:4
            if ~isequal(Gamm(lam,mu,nu),sym(0))
                txt=sprintf('Gamma_%d%d^%d',mu,nu,lam);
                disp(str2sym(txt)==Gamm(lam,mu,nu))
            end
        end, end, end

%
%% 3 Schwarzschild metric B(r), A(r)
% Equivalence principle:  $\frac{{\partial {p_\nu }}}{{\partial {q^\mu }}} = {m_0}\frac{{\partial {A_\nu }}}{{\partial {q^\mu }}} =  - {m_0}{g_{\mu \lambda }}\Gamma _{\kappa \nu }^\lambda {v^\kappa }$
% for $\nu  = 4,\;\mu  = 1,\;{A_4} = \frac{M_g}{c}\frac{G_g}{r}\; \to $  ${m_0}\frac{{MG}}{r} =  - {m_0}{G_{11}}\Gamma _{1\kappa }^4{v^\kappa }$
% with Schwarzschild radius r_s

syms m_0 M_g G_g r_s positive
v_v=[0;0;0;c]

A_4(r)=M_g*G_g/(r*c)
m_0*diff(A_4)==-m_0*G(1,1)*Gamm(1,:,4)*v_v
e1(r)=isolate(ans,diff(B,r))
e2=r_s==sube(rhs(e1(r)),r==1)
sube(e1,isolate(e2,M_g))
e3=B(r)==dsolve(ans,B(inf)==1)

%
% Radial motion on light cone: 0=(ds)^2=A(r)^2(dr)^2-B(r)^2(dct)^2\ \to \
% 1=\frac{dr}{dct}=\frac{B(r)}{A(r)}

A(r)*B(r)==1
eq3=[e3,isolate(sube(ans,e3),A(r))]

%
%% 4 Resting clock in the gravitational field
% (ds)^2=A(r)(dr)^2
% r^2(d\vartheta)^2+
% r^2\sin^2\vartheta(d\varphi)^2-B(r)(dct)^2
% For object at rest: (ds)^2=G_{44}(dct)^2=-(d\tau)^2
% Interpretation:  A local time interval d\tau at the event horizon appears infinitely long from the outside. Frequencies appear infinitely small from the outside.

syms dt
dtau=sqrt(sube(-G(4,4),eq3))*dt
sube(dtau,[dt==1,r_s==1])
fplot(ans,[0 10])
axis([0 10 0 1])
xlabel('Radial coordinate r [r_s]')
ylabel('d\tau [dt]')

%
%% 5 Black hole
% The Schwarzschild elongation can be represented by a curvature of the space in an additional space dimension H. The degree of bulge in H direction is determined by the theorem of Pythagoras.
%
% Interpretation: viewed from the outside, a local length dr is stretched to \sqrt{g_{rr}}dr. The path to (or from) the event horizon ultimately becomes infinitely large. Light cannot cover this way. Therefore, the event horizon appears black.
% A local observation (for example by an astronaut at the event horizon) is not affected.

syms dr dH
dr^2+dH^2==G(1,1)*dr^2
sube(ans,eq3)
isolate(ans,dH)
sube(dH,ans)
int(ans/dr);
H=sube(ans,r_s==1)

%

figure
set(gcf,'DefaultTextFontSize',18)
fplot(H,[0 10],'Color','b')
xlabel('Radial coordinate [r_s]')
hold on
Hn=matlabFunction(H);
plot([2 4 4],Hn([2 2 4]),'r')
ylabel('H(r) [r_s]')
text(1.9,2.5,'$\sqrt{g_{rr}} dr$','Interpreter','latex','Color','r',...
    'Rotation',40)
text(3,1.8,'$dr$','Color','r','Interpreter','latex', 'HorizontalAlignment','center')
text(4.1,mean(Hn([2 4])),'$dH$','Color','r','Interpreter','latex')

axis([0 10 0 6])
hold off

%
% The so-called Flamm'sche paraboloid is created by rotation-symmetrical representation of H(r)

[r,phi]=meshgrid(linspace(1.01,5,25),linspace(0,2*pi,25));
matlabFunction(H);
surfc(r.*cos(phi),r.*sin(phi),ans(r))
axis equal
view([-50.30 30.00])
jet_=jet;
colormap(jet_(end:-1:1,:))
axis([-5 5 -5 5 0 4])
zlabel('H(r) [r_s]')
title('Flame paraboloid')
view([-51.72 38.54])