function func = aren_mov ( outdir );
% Aren_Mov
%
% 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
%
% Movie output for three body problem ("Arenstorf orbit")
%
% Parameter:
%   outdir (input)  : directory for output of animation figures
%
% Examples:
%   aren_mov ( 'h:\tmp\output' );    % print -dpng ../png/aren_fig99.png
%
% Movie in mpg format:
%   The Matlab command
%     aren_mov ( 'output' );
%   results in files .\output\fig00.eps, ..., .\output\fig99.eps
%   that are summarized to a movie using the ImageMagick package from
%     http://www.imagemagick.org
%   by the command
%     cd output
%     convert *.eps aren5.mpg
%   Because of known problems in the generation of mpg files with the
%   Windows version of convert the Linux or UNIX versions should be used.
%
% Width and height of the mpg video are smaller by about 1/4 than the
% original Matlab output. Resize the movie by factor 133% when importing
% the mpg video to PowerPoint.

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

% -> Euler on coarse grid, slide version
heulc = 5.0e-4;
neulc = round ( te / heulc );
teulc = linspace ( 0, te, neulc+1 );
yeulc = zeros ( neulc+1, 4 );
yeulc(1,:) = y0';
for ieul=1:neulc
    yn = yeulc(ieul,:)';
    k1 = aren_rhs ( teulc(ieul), yn, mu );
    ynp = yn + heulc*k1;
    yeulc(ieul+1,:) = ynp';
end;

% -> Matlab solver ode45 with TOL=1e-4, dense output, movie version
nout   = 100;
eval ( [ '! mkdir ', outdir ] );
eval ( [ '! del ', outdir, '\fig*.eps' ] );

optdop = odeset ( 'RelTol', 1e-4, 'AbsTol', 1e-4, 'Stat', 'on' );
[ tdop, ydop ] = ode45 ( @aren_rhs, [ 0 te ], y0, optdop, mu );

ntact = length(tdop);
nffac  = 535 / length(tdop);
for iout=1:min(ntact,nout);

    % -> output: reference solution etc.  
    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' );

    % -> movie
    if iout<min(ntact,nout),
        irange = 1:iout*ntact/min(ntact,nout);
    else
        irange = 1:ntact;
    end; 
    hplotdop = plot ( ydop(irange,1), ydop(irange,3), 'b', 'LineWidth', ActLineWidth );
    crange = 1:nffac*max(irange);
    hploteul = plot ( yeulc(crange,1), yeulc(crange,3), 'r', 'LineWidth', ActLineWidth  );
    if iout>=min(ntact,nout)-1,
        plot ( yeulc(:,1), yeulc(:,3), 'r', 'LineStyle', '--', 'LineWidth', 1  );
    end;
   
    hold off,

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

    % -> legend
    hlegend = legend ( [ hplotdop hploteul ]         ...
                     , ' TOL = 1.0e-4,  ode45'       ...
                     , ' h = 5.0e-4,  Euler'         ...
                     );

    % -> 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' )';
    end;

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

    % -> output to file
    set ( gcf, 'PaperUnits',    'centimeters' );
    set ( gcf, 'PaperPosition', [ 0 0 graphwidth graphheight ] );
    drawnow,
    eval ( [ 'print -depsc ', outdir, sprintf('\\fig%02i.eps',iout-1) ] );
   
end;