%% Beispiel zur Gravitationsabstrahlung
% 

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016
% Vers.1.0 02.08.2016

close all, clear, clc
%% Parameter und Konstanten
G = 6.67384e-11;  % Gravitationskonstante
c = 299792458;    % m/s
m2s = 1.983e30;   % Masse der Sonne in kg
m2 = 3.e6*m2s;    % Masse des schwarzen Lochs
%                   m1 = 10*m2;       % Masse des Sterns  ACHTUNG
%                                       Tippfehler im Buch
m1 = 10*m2s;      % Masse des Sterns 
AE = 1.496e11;    % Mittlere Entfernung Erde-Sonne in m = 1 astronomische Einheit (1AE)
%%
AE = 1*AE;
epsilon = 0.2056;
mu = m1*m2/(m1+m2);    % Reduzierte Masse
Vorf = G*m1*m2;
%% Startwerte
r0 = (1+epsilon)*AE;   % Anfangswert: Aphel
theta0 = 0;      
pr0=0;
%
E0 = -Vorf/(2*r0) * (1+epsilon); % Energie
pt0 = sqrt((E0*r0+Vorf)*2*mu*r0);
%% Differentialgleichung loesen
y0 = [r0;theta0;pr0;pt0];
tmax=sqrt(4*pi^2/(G*(m1+m2)))*r0^(3/2)*6;      % *n Perioden
options = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',360);
tic
[t,y] = ode45(@(t,y) keplerDGL(t,y,mu,Vorf),[0,tmax],y0,options);
toc
%% Visualisieren
x1=y(:,1).*sin(y(:,2));
x2=y(:,1).*cos(y(:,2));
plot(x1,x2),hold on, plot(x1(1),x2(1),'pr'),plot(0,0,'yp'),axis equal, shg
%% Leistung: Gravitationswelle
dE = 8*G^4*m1^2*m2^2*(m1+m2)/(15*r0^5*c^5*(1-epsilon^2)^5) * (1+epsilon.*cos(y(:,2))).^4 ...
      .*(12*(1+epsilon.*cos(y(:,2))).^2 +epsilon^2.*sin(y(:,2)).^2);
dE_mittel = -32*G^4*m1^2*m2^2*(m1+m2)/(5*r0^5*c^5*(1-epsilon^2)^(7/2)) * ...
             (1 + 73/24*epsilon^2 + 37/96*epsilon^4) 
figure
plot(t,dE),shg
%%
% t und r0 skalieren
tneu = t.*t(end)./(3*(t+t(end)));
r0neu = r0*sqrt((t(end)./(t+t(end))));
dEneu = 8*G^4*m1^2*m2^2*(m1+m2)./(15*r0neu.^5*c^5.*(1-epsilon^2)^5) .* (1+epsilon.*cos(y(:,2))).^4 ...
      .*(12*(1+epsilon.*cos(y(:,2))).^2 +epsilon^2.*sin(y(:,2)).^2);
figure
plot(tneu,dEneu),shg
figure
x1n=y(:,1).*r0neu.*sin(y(:,2))/r0;
x2n=y(:,1).*r0neu.*cos(y(:,2))/r0;
plot(x1n,x2n),hold on,axis equal, shg
%% Buchvisualisierung
figure
h(1)=subplot(2,1,1);
hp1=plot(t,dE,'k'); axis tight
hp1.LineWidth=2; ylabel('-dE'), xlabel('t')
%
h(2)=subplot(2,1,2);
hp2 = plot(tneu,dEneu,'k');axis tight, hold on
hp2.LineWidth=2; ylabel('-dE')
ha = axes('Position',[0.275,0.25,0.22,0.25]);
plot(x1n,x2n,'k'), axis equal
ha.XTickLabel={};
ha.YTickLabel={};
set(h,'FontSize',12);
