%
%   Lösung des Radial-Anteils x der DGL für die 
%     Bewegung von Planeten um die Sonne: 
%     d^2x/dt^2 = L^2 / (m^2 * x^3) - k / x^2
%
%   Erde - ungefähre Daten in km und Tagen d: 
%    m = 5.974*10^24   [kg]  - Erdmasse 
%    r = 149.6*10^6    [km]  - mittlerer Abstand (1 AE) 
%    L = m*r^2*2*pi/365   [kg km^2/d] = m*r^2*(2*pi/Jahr) - ungef. Drehimpuls  
%    k = 1.327*10^20 * (24*60*60)^2/10^9   [km^3/d^2] = G * mS  
%           G = 6.67*10-11 [m^3/kg s^2], mS = 1.989*10^30  [kg]  
%    x0 = 0.983*149.6*10^6   [km] = 0.983 AE - Perihel
%    v0 = 0                  [km/d]          - mittlere Radius-Änderung = 0 
% 
%    Perihel = 0.983 AE 
%    Aphel   = 1.017 AE 
%    Aphel - Perihel = 5.1*10^6 km 
% 
%   >> PlanetDGL( m, L, k, x0, v0, 365 )
%
function PlanetDGL( m, L, k, x0, v0, tm ) % 1/2 P  -> 9 + 5 P
   
   w0 = [ x0, v0 ];                         % 1/2 P
   tspan = [0,tm];                          % 1/2 P
   
   [t,w] = ode45( @dglFnc, tspan, w0, [], L, m, k );     % 2 P

   x = w( :, 1 );                           % 1/2 P
   v = w( :, 2 );                           % 1/2 P  
   subplot( 1,2,1 )
   plot( t, x );                            % 1 P 
   subplot( 1,2,2 )
   plot( t, v );                            % 1 P 
   
   maxR = max( x );                         % 1 P
   minR = min( x );   	                    % 1 P
   diff = maxR - minR;
   fprintf( 'Perihel-Aphel-Abstand: %g \n', diff );  % 1/2 P 
   
   AE = 149.6*10^6;
   Aphel_AE   = maxR / AE   % in AE 
   Perihel_AE = minR / AE 

end

function dy_dt = dglFnc( t, y, L, m, k ) 
   dy_dt(1,1) = y(2);                          % 1 P
   dy_dt(2,1) = L^2/(m^2*y(1)^3) - k/y(1)^2;   % 4 P;
end
