function func = demoosc ( iplot );
% DemoOsc
%
% 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 :  Dec 6, 2004
%
% Demonstration linear mechanic Oscillations
%
% Parameters:
%   iplot  (input)  : control flag "plot"
%                       =0 .. resonance
%                       =1 .. periodic case
%                       =2 .. aperiodic case
%                       =3 .. aperiodic limit case
%
% Examples:
%   demoosc ( 0 );   %   print -dpng ../png/demoosc0.png
%   demoosc ( 1 );   %   print -dpng ../png/demoosc1.png
%   demoosc ( 2 );   %   print -dpng ../png/demoosc2.png
%   demoosc ( 3 );   %   print -dpng ../png/demoosc3.png

graphwidth  = 14.0;
graphheight =  8.0;
fac         =  0.9;
shift       =  0.4;
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+2.0*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
tmin =  0.0;
tmax =  1.0;
xmin = -1.0;
xmax =  1.0;

ommin =  0.0;
ommax =  1.0;
Vmin  =  0.0;
Vmax  = 30.0;

% -> plot
tt   = linspace ( tmin, tmax, 1001 );
omom = linspace ( ommin, ommax, 1001 );
if iplot==0,
    om0   = 0.4;
    VV0   = 1 ./ sqrt ( ( om0^2 - omom.^2 ).^2 + 4*0.025^2*omom.^2 );
    VV1   = 1 ./ sqrt ( ( om0^2 - omom.^2 ).^2 + 4*0.050^2*omom.^2 );
    VV2   = 1 ./ sqrt ( ( om0^2 - omom.^2 ).^2 + 4*0.100^2*omom.^2 );
    VV3   = 1 ./ sqrt ( ( om0^2 - omom.^2 ).^2 + 4*0.150^2*omom.^2 );
    VV4   = 1 ./ sqrt ( ( om0^2 - omom.^2 ).^2 + 4*0.200^2*omom.^2 );
elseif iplot==1.,
    alpha =  2;
    om0   = 24;
    om1   = sqrt ( om0^2 - alpha^2 );
    xx    = 0.9*exp(-alpha*tt).*cos(om1*tt+pi/6);
    xxmax = exp(-alpha*tt);
elseif iplot==2,
    alpha = 4.0;
    om0   = 3.6;
    beta  = sqrt ( alpha^2 - om0^2 );
    xx0   =  0.7*exp((-alpha+beta)*tt);
    xx1   = -0.9*exp((-alpha-beta)*tt);
    xx2   = xx0 + xx1;
    xx3   = 2.1*xx0 + xx1;
    xx4   = 2.0*xx0 + 0.6*xx1;
elseif iplot==3,
    alpha = 4.0;
    om0   = alpha;
    xx0   =  0.9*exp(-alpha*tt);
    xx1   = ( 1  +  tt).*xx0;
    xx2   = ( 0.9-10*tt).*xx0;
    xx3   = (-1  + 8*tt).*xx0;
end;

plot ( [ tmin tmax ], [ 0 0 ], '-k' );
hold on;
if iplot==0,
    plot ( om0*[1 1], [ Vmin Vmax ], '--r', 'LineWidth', 2 );
    plot ( omom,    VV0, '-b', 'LineWidth', ActLineWidth );
    plot ( omom,    VV1, '-b', 'LineWidth', ActLineWidth );
    plot ( omom,    VV2, '-b', 'LineWidth', ActLineWidth );
    plot ( omom,    VV3, '-b', 'LineWidth', ActLineWidth );
    plot ( omom,    VV4, '-b', 'LineWidth', ActLineWidth );
elseif iplot==1,
    plot ( tt,   xxmax, '--r', 'LineWidth', 2 );
    plot ( tt,  -xxmax, '--r', 'LineWidth', 2 );
    plot ( tt,   xx,     '-b',  'LineWidth', ActLineWidth );
elseif iplot==2,
    plot ( tt,      xx1, '-b', 'LineWidth', ActLineWidth );
    plot ( tt,      xx2, '-b', 'LineWidth', ActLineWidth );
    plot ( tt, -2.5*xx2, '-b', 'LineWidth', ActLineWidth );
    plot ( tt,      xx3, '-b', 'LineWidth', ActLineWidth );
    plot ( tt,      xx4, '-b', 'LineWidth', ActLineWidth );
elseif iplot==3,
    plot ( tt,      xx1, '-b', 'LineWidth', ActLineWidth );
    plot ( tt,      xx2, '-b', 'LineWidth', ActLineWidth );
    plot ( tt,      xx3, '-b', 'LineWidth', ActLineWidth );
end;
hold off;

if iplot==0,
    axis ( [ ommin ommax Vmin Vmax ] );
    xlabel ( '\omega' );
    ylabel ( 'V_\alpha(\omega)' );
    title ( [ 'Amplitudenverstaerkung, \omega_0 = ', sprintf('%3.1f',om0) ] );
elseif ( (iplot==1) | (iplot==2) | (iplot==3) ),
    axis ( [ tmin tmax xmin xmax ] );
    xlabel ( 't' );
    ylabel ( 'x_h' );
    title ( [ '\alpha = ', sprintf('%3.1f ,  ',alpha), '\omega_0 = ', sprintf('%3.1f',om0) ] );
end;

%
% -> Korrektur der Groesse und Lage der Koordinatensysteme
allax = get ( gcf, 'Children' )';
for ax = allax,
  rect    = get ( ax, 'Position' );
  rect(2) = rect(2) + shift * (1-fac) * rect(4);
  rect(4) = fac * 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;

