function [ tout, yout ] = aren_rk ( imeth ),
% Aren_RK
%
% 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 2006/07
%
% Author :      Prof. Dr. M. Arnold, martin.arnold@mathematik.uni-halle.de
% Version of :  Nov 24, 2006
%
% Explicit Runge-Kutta methods applied to benchmark "Arenstorf orbit"
%
% Parameters:
%   imeth  (input)  : control parameter method (empty input and vector valued input accepted)
%                       =1 .. explicit Euler,        h=0.0005
%                       =2 .. explicit Euler,        h=0.00002
%                       =3 .. classical Runge-Kutta, h=0.002
%                       =4 .. DOPRI5 by ode45,       TOL=1e-4, output at t_n
%                       =5 .. DOPRI5 by ode45,       TOL=1e-4, dense output
%   tout   (output) : numerical solution
%   yout   (output) : numerical solution
%
% Examples:
%   [ tout, yout ] = aren_rk ( [] );           % print -dpng ../png/aren_rk0.png
%   [ tout, yout ] = aren_rk ( 1 );            % print -dpng ../png/aren_rk1.png
%   [ tout, yout ] = aren_rk ( 2 );            % print -dpng ../png/aren_rk2.png
%   [ tout, yout ] = aren_rk ( [ 1 3 ] );      % print -dpng ../png/aren_rk3.png
%   [ tout, yout ] = aren_rk ( 4 );            % print -dpng ../png/aren_rk4.png
%   [ tout, yout ] = aren_rk ( 5 );            % print -dpng ../png/aren_rk5.png

% -> control data plot
graphwidth  = 19.0;
graphheight = 12.0;
vfac        =  0.99;
vshift      =  0.8;
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 = -2.0;
y1max =  2.5;
y2min = -1.5;
y2max =  1.5;

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

% -> evaluate numerical solution for all methods in parameter imeth
nrun = length(imeth);

for irun=1+(1:nrun),

    imact = imeth(irun-1);
    
    if ( (imact==1) | (imact==2)  ),                      % Euler

        if imact==1,
            heul = 5.0e-4;
            hstr = '5.0e-4';
        elseif imact==2,
            heul = 2.0e-5;
            hstr = '2.0e-5';
        end;
        neul = round ( te / heul );
        teul = linspace ( 0, te, neul+1 );
        yeul = zeros ( neul+1, 4 );
        yeul(1,:) = y0';
        for ieul=1:neul
            yn = yeul(ieul,:)';
            k1 = aren_rhs ( teul(ieul), yn, mu );
            ynp = yn + heul*k1;
            yeul(ieul+1,:) = ynp';
        end;
        tout{irun} = teul;
        yout{irun} = yeul;
        colout{irun} = 'r';
        mstr{irun} = [ ' h = ', hstr, ',  Euler' ];

    elseif (imact==3),                                    % classical Runge-Kutta
        hrk = 2.0e-3;
        nrk = round ( te / hrk );
        trk = linspace ( 0, te, nrk+1 );
        yrk = zeros ( nrk+1, 4 );
        yrk(1,:) = y0';
        for irk=1:nrk
            yn  = yrk(irk,:)';
            yn1 = yn;
            k1  = aren_rhs ( trk(irk)      , yn1, mu );
            yn2 = yn + hrk/2 * k1;
            k2  = aren_rhs ( trk(irk)+hrk/2, yn2, mu );
            yn3 = yn + hrk/2 * k2;
            k3  = aren_rhs ( trk(irk)+hrk/2, yn3, mu );
            yn4 = yn + hrk   * k3;
            k4  = aren_rhs ( trk(irk)+hrk  , yn4, mu );
            ynp = yn + hrk/6 * ( k1 + 2*k2 + 2*k3 + k4 );
            yrk(irk+1,:) = ynp';
        end; 
        tout{irun} = trk;
        yout{irun} = yrk;
        colout{irun} = 'g';
        mstr{irun} = ' h = 2.0e-3,  RK';

    elseif ( (imact==4) | (imact==5) ),

        if imact==4,
            optdop = odeset ( 'RelTol', 1e-4, 'AbsTol', 1e-4, 'Stat', 'on', 'Refine', 1 );
        elseif imact==5,
            optdop = odeset ( 'RelTol', 1e-4, 'AbsTol', 1e-4, 'Stat', 'on' );
        end;
        [ tdop, ydop ] = ode45 ( @aren_rhs, [ 0 te ], y0, optdop, mu );
        tout{irun} = tdop;
        yout{irun} = ydop;
        colout{irun} = 'b';
        mstr{irun} = ' TOL = 1.0e-4,  ode45';

    end;

end;

% -> output
clf,
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' );

hrun = zeros(nrun,1);
for irun=1+(1:nrun),
    hrun(irun-1) = plot ( yout{irun}(:,1), yout{irun}(:,3), colout{irun}, 'LineWidth', ActLineWidth );
    lstr{irun-1} = mstr{irun};
end;   
   
hold off,

xlabel ( 'y_1' );
ylabel ( 'y_2' );
title ( 'Arenstorf-Orbit' );

% -> legend
if nrun>0,
    hlegend = legend ( hrun, lstr );
end;

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

%
% -> 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;

% -> enforce correct position of legend
if nrun>0,
    set ( hlegend, 'Location', 'Best' );
end;
