function func = burgers ( iplot ),
% Burgers
%
% This m-file is part of the material of the lectures on
%   Numerical Solution of Partial Differential Equations
% at the Math Department of Halle University, Germany.
%
% Author :      Martin Arnold, martin.arnold@mathematik.uni-halle.de
% Version of :  Apr 11, 2006
%
% Plot output data for Burgers equation
%
% Parameters:
%   iplot  (input)  : control parameter output
%                       =1 ... plot characteristics
%                       =2 ... plot solution
%
% Examples:
%   burgers ( 1 );  %  print -dpng ../png/burgers1.png
%   burgers ( 2 );  %  print -dpng ../png/burgers2.png

% -> control data plot
graphwidth   = 12.0;
graphheight  = 12.0;
fac	     =  0.9;
shift	     =  0.5;
border       =  2.0;

ActFontSize  = 16;
ActLineWidth =  2;

%
% -> 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' );
screenscale = 1.0;

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

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

% control parameters output
nxx = 81;
xxmin = -4.5;
xxmax =  3.5;
xx = linspace ( xxmin, xxmax, nxx );

xmin = -1.5;
xmax =  2.5;
tmin =  0.0;
tmax =  2.0;

% compute the shock
tend = tmax * ones ( size(xx) );                                                    
xshock = [];                                                                        
tshock = [];                                                                        
nxxneg = length ( xx(xx<0) );                                                       
for ixx=nxxneg+1:nxx,                                                               
   xxact = xx(ixx);                                                                 
   if xxact==0,                                                                     
      tend(ixx) = 1.0;                                                              
      xshock = [ xshock; xxact+(pi/4-atan(xxact))*tend(ixx) ];                      
      tshock = [ tshock; tend(ixx) ];                                               
   elseif xxact>0,                                                                  
      ixxend = 0;                                                                   
      for ixxneg=1:nxxneg,                                                          
         if tend(ixxneg)==tmax,                                                     
            xxneg = xx(ixxneg);                                                     
            tcross = ( xxneg - xxact ) / ( atan(xxneg) - atan(xxact) );             
            if (tcross>0) & (tcross<tend(ixx)),                                     
               tend(ixx) = tcross;                                                  
               ixxend = ixxneg;                                                     
            end;                                                                    
         end;                                                                       
      end;                                                                          
      if ixxend>0,                                                                  
         tend(ixxend) = tend(ixx);                                                  
         xshock = [ xshock; xx(ixxend)+(pi/4-atan(xx(ixxend)))*tend(ixxend) ];      
         tshock = [ tshock; tend(ixxend) ];                                         
      end;                                                                          
   end;                                                                             
end;                                                                                

% plot
if iplot==1,   
   for ixx=1:nxx,
      xxact = xx(ixx);
      plot ( xxact+(pi/4-atan(xxact))*[0,tend(ixx)], [0,tend(ixx)] );
      hold on;
   end;
   plot ( xshock, tshock, '-r', 'LineWidth', 2 );
   hold off;
   axis ( [ xmin xmax tmin tmax ] );
   xlabel ( 'x' );
   ylabel ( 't' );
   title ( 'Charakteristiken fuer AB  u(x,0) = \pi/4 - atan(x)' );
elseif iplot==2,
   x0neg = -5.0:0.001:0.0;
   s0neg = pi/4 - atan(x0neg);
   x0pos = 0.0:0.001:5.0;
   s0pos = pi/4 - atan(x0pos);
   nt = 51;   nx = 101;
   tt = linspace ( tmin, tmax, nt );
   xx = linspace ( xmin, xmax, nx );
   u = zeros ( nx, nt );
   for ix=1:nx,
      xact = xx(ix);
      for it=1:nt,
         tact = tt(it);
         if tact>=4*xact/pi,
            [ auxd, ix0 ] = min ( abs(x0neg+s0neg*tact-xact) );
            x0act = x0neg(ix0);
         else
            [ auxd, ix0 ] = min ( abs(x0pos+s0pos*tact-xact) );
            x0act = x0pos(ix0);
         end;
         u(ix,it) = pi/4 - atan(x0act);
      end;
   end;
   mesh ( xx, tt, u');
   axis ( [ xmin xmax tmin tmax -pi/4 3*pi/4 ] );
   xlabel ( 'x' );
   ylabel ( 't' );
   zlabel ( 'u(x,t)' );
   title ( 'Loesung der Burgers-Gleichung  u_t+uu_x=0' );
   view ( -37.5+45.0, 45.0 );
end;

%
% -> Korrektur der Groesse und Lage der Koordinatensysteme, Aenderung der Schriftgroesse
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 );
  set ( ax, 'FontSize', ActFontSize );
  set ( get ( ax, 'XLabel' ), 'FontSize', ActFontSize );
  set ( get ( ax, 'YLabel' ), 'FontSize', ActFontSize );
  set ( get ( ax, 'ZLabel' ), 'FontSize', ActFontSize );
  set ( get ( ax, 'Title' ), 'FontSize', ActFontSize );
end;