function aren_ode45,
% Aren_ode45
%
% This m-file is part of the material of the lecture
%   Theory and Numerics of Ordinary Differential Equations
% at the Math Department of Halle University, Germany.
% Winter term 2004/05
%
% Author :      Prof. Dr. M. Arnold, arnold@mathematik.uni-halle.de
% Version of :  Nov 17, 2004
%
% Matlab integrator ode45 applied to benchmark "Arenstorf orbit"
%
% Parameters:
%   none
%
% Examples:
%   aren_ode45;    %  print -dpng ../png/aren_ode45.png

% -> control data plot
graphwidth  = 28.0;
graphheight = 15.0;
vfac        =  1.0;
vshift      =  0.0;
hfac        =  1.0;
hshift      =  0.0;
papershift  =  0.0;

ActFontSize   = 14;
ActLineWidth  =  2;
ActMarkerSize = 12;
RatioFlag     =  0;

border      = 1.0;

%
% -> Bestimme Groesse des Bildschirms und Skalierungsfaktor fuer die Ausgabe
root = 0;
set ( root, 'Units', 'centimeters' );
screensize = get ( root, 'ScreenSize' );
left   = screensize(1);
bottom = screensize(2);
width  = screensize(3);
height = screensize(4);
set ( root, 'Units', 'pixels' );
if RatioFlag==0,
  screenscale = 1.0;
elseif RatioFlag==1,
  screenscale = (width-2*border) / graphwidth;
end;

%
% -> Definiere den Bildschirmausschnitt
if length ( get ( root, 'Children' ) ) == 0,
  h = figure;
else
  h = gcf;
  clf;
end;
set ( h, 'Units',    'centimeters' );
set ( h, 'Position', ...
  [ left+width-screenscale*graphwidth-border height-2*border-screenscale*graphheight ...
    screenscale*graphwidth                   screenscale*graphheight ] );
set ( h, 'Units',    'pixels' );

%
% -> Skalierung der Ausgabe in PostScript
%
set ( gcf, 'PaperUnits',    'centimeters' );
set ( gcf, 'PaperPosition', [ papershift papershift graphwidth graphheight ] );

% -> data
y0 = [ 0.994                    ...
     , 0.0                      ...
     , 0.0                      ...
     , -2.0015851063790825      ...
     ]';
te = 17.06521656015796255889;
mu = 0.012277471;

ntref = 1001;   
tsref = linspace ( 0, te, ntref )';

rtolref = 1.0e-8;
atolref = 1.0e-8;

% -> control parameters output
y1min = -1.5;
y1max =  1.5;
y2min = -1.5;
y2max =  1.5;

xmin = -1;
xmax = 2*te + 1;
ymin = -1.5;
ymax = 1.5;
hmin = 1.0e-4;
hmax = 1.0;

% -> evaluate reference solution
optref = odeset ( 'RelTol', rtolref, 'AbsTol', atolref );
[ tref, yref ] = ode45 ( @aren_rhs, tsref, y0, optref, mu );

% -> evaluate numerical solution
optdop = odeset ( 'RelTol', 1e-4, 'AbsTol', 1e-4, 'Stat', 'on', 'Refine', 1 );
[ tdop, ydop ] = ode45 ( @aren_rhs, [ 0 te ], y0, optdop, mu );
colout = 'b';

[ tstep, ystep ] = ode45 ( @aren_rhs, [ 0 2*te ], y0, optdop, mu );

% -> output
clf,
subplot ( 1, 2, 1 );
plot ( -mu , 0, 'om', 'Markersize', 7, 'MarkerFaceColor', 'm' );
hold on,
plot ( 1-mu, 0, 'om', 'Markersize', 4, 'MarkerFaceColor', 'm' );
plot ( [ y1min y1max ], [ 0 0 ], ':k' );
plot ( [ 0 0 ], [ y2min y2max ], ':k' );

plot ( yref(:,1), yref(:,3), '-k' );
hplotdop = plot ( ydop(:,1), ydop(:,3), colout, 'LineWidth', ActLineWidth );

xlabel ( 'y_1' );
ylabel ( 'y_2' );
title ( 'Arenstorf-Orbit: Integration mit ode45 (TOL=1e-4)' );

% -> scaling of axis etc.
axis ( [ y1min y1max y2min y2max ] );
set ( gca, 'DataAspectRatio'   , [ 1 1 1 ] );

subplot ( 2, 2, 2 );
plot ( tstep, ystep, 'LineWidth', ActLineWidth );
axis ( [ xmin xmax ymin ymax ] );
xlabel ( 't' );
ylabel ( 'y' );
title ( 'Loesung  ( 2 Perioden )' );

subplot ( 2, 2, 4 );
semilogy ( tstep, diff([tstep(1);tstep]), '-b'   ...
         , tstep, diff([tstep(1);tstep]), '+b'   ...
         );
axis ( [ xmin xmax hmin hmax ] );
xlabel ( 't' );
ylabel ( 'h_m' );
title ( 'Schrittweitenfolge  ode45  (TOL=1e-4)' );

%
% -> Korrektur der Groesse und Lage der Koordinatensysteme
allax = get ( gcf, 'Children' )';
for ax = allax,
  rect    = get ( ax, 'Position' );
  rect(1) = rect(1) + hshift * (1-hfac) * rect(3);
  rect(3) = hfac * rect(3);
  rect(2) = rect(2) + vshift * (1-vfac) * rect(4);
  rect(4) = vfac * rect(4);
  set ( ax, 'Position', rect );
end;

%
% -> Aenderung der Schriftgroesse etc.
for hf = allax,
  set ( hf, 'FontSize', [ActFontSize] );
  set ( get ( hf, 'title' ), 'FontSize', [ActFontSize] );
  set ( get ( hf, 'xlabel' ), 'FontSize', [ActFontSize] );
  set ( get ( hf, 'ylabel' ), 'FontSize', [ActFontSize] );
  allch = get ( hf, 'Children' )';
  % for hch = allch,
  %   if strcmp(get(hch,'Type'),'line'),
  %     set ( hch, 'LineWidth',  ActLineWidth );
  %     set ( hch, 'MarkerSize', ActMarkerSize );
  %   end;
  % end;
end;
