function wvchar(action);
%WVCHAR - Characteristics of the Wave Equation
%
%         <DESCRIPTION>
%         wvchar is a graphical interface for illustrating the use of 
%         characterics in construction the d'Alembert solution of the 
%         Initial Boundar Value Problem 
%         
%               PDE: u_tt = c^2 u_xx,  0 < x < 1, t > 0
%
%               BC:  u(0,t)=0=u(1,t) or u_x(0,t)=0=u_(1,t)
%
%               IC:  u(x,0) = f(x), u_t(x,0) = 0.
%
%         <TO RUN>
%         Type wvchar at the MATLAB prompt with no options.
%
%         <TO USE>
%         Once started, the following fields may be edited:
%
%            Time Interval      - Plot xt plane for 0 < t < Final Time 
%            Speed              - The wave speed c in the wave equation
%            Boundary Condition - Radio buttons to select u(0,t)=u(1,t)=0 or
%                                 u_x(0,t)=u_x(1,t)=0 as the boundary conditions.
%            Initial condition  - The function of x defining the initial 
%                                 condition u(x,0).
%
%            Restart Button - Redraw plot area after setting the time interval,
%                             speed, boundary condition, and initial condition.
%         
%          Once set up, clicking the point (x,t) in the plot area constructs
%          characteristic lines which extend back to the x-axis. The value of 
%          u(x,0) at the starting point of these characteristics is used to
%          construct the value of u(x,t) at (x,t) by d'Alembert's formula. 

%Created:  5/95 by R. Knobel

if nargin<1,
   action='start';
end;

% Default X interval is [0,1]

MinX = 0;
MaxX = 1;

% Default initial time is t=0

MinT = 0;

% If this is not the inital setup, get stored graphics handles

if ~strcmp(action,'start'),
   handles=get(gcf,'UserData');
   hTmax = handles(1);
   hChar1= handles(2);
   hChar2= handles(3);
   hClick= handles(4);
   hBC1  = handles(5);
   hBC2  = handles(6);
   hRoot1= handles(7);
   hRoot2= handles(8);
   hFunction = handles(9);
   hU     = handles(10);
   hSpeed = handles(11);
   hAxes  = handles(12);
end;

% If this is the initial startup, create the figure, controls, and 
% and initial plot:

if strcmp(action,'start'),
   
   set(gcf, ...
      'NumberTitle','off', ...
      'Name','Characteristics of the Wave Equation', ...
      'Backingstore','off',...
      'Units','normalized');
   
   bgcolor = get(gcf,'Color');
   textcolor=[0 0 0];
   if max(bgcolor)<.5,
      textcolor=[1 1 1];
   end;
   
   % Console Frame 
   
   top     =0.95;
   bottom  =0.05;
   left    =0.81;
   btnWid  =0.13;
   btnHt   =0.06;
   spacing =0.02;
   frmBorder=0.02;
   
   frmPos=[left-frmBorder bottom-frmBorder btnWid+3*frmBorder ...
         top-bottom+2*frmBorder]; 
   uicontrol( ...
      'Style','frame', ...
      'Units','normalized', ...
      'Position', frmPos, ...
      'ForegroundColor','white',...
      'BackgroundColor','black');
   
   % Time Interval
   
   uicontrol( ...
      'Style','text',...
      'Units','normalized',...
      'Position', [left top-.05 btnWid .045],...
      'BackgroundColor','black',...
      'ForegroundColor','red',...
      'String','Final Time t');
   MaxT = 2;
   hTmax=uicontrol( ...
      'Style','edit',...
      'Units','normalized',...
      'Position', [.85 top-.10 .06 .045],...
      'BackgroundColor',[1 1 1],...
      'String',num2str(MaxT));
   
   % Characteristic Speed c in u_tt = c^2 u_xx
   
   SpeedPos = .80;
   uicontrol(...
      'Style','text',...
      'String','Wave',...
      'Units','normalized',...
      'Position',[left SpeedPos-.04 btnWid .04],...
      'Background','black',...
      'Foreground','red');
   uicontrol(...
      'Style','text',...
      'String','Speed',...
      'Units','normalized',...
      'Position',[left SpeedPos-.08 btnWid .04],...
      'Background','black',...
      'Foreground','red');
   
   c = 0.5;
   uicontrol('Style','text',...
      'Units','normalized',...
      'Position',[.83 SpeedPos-.13 .04 .05],...
      'BackgroundColor','black',...
      'ForegroundColor','white',...
      'String','c =');
   hSpeed = uicontrol(...
      'Style','edit',...
      'Position',[.88 SpeedPos-.13 .07 .05],...
      'Units','normalized',...
      'BackgroundColor','white',...
      'UserData',c,...
      'String',num2str(c));
   
   % Boundary Condition Radio Control 
   
   uicontrol(...
      'Style','frame',...
      'Units','normalized',...
      'Position',[.81 .38 .15 .18],...
      'Background','black',...
      'Foreground','red');
   uicontrol(...
      'Style','text',...
      'String','Boundary',...
      'Units','normalized',...
      'Position',[.83 .51 .11 .04],...
      'Background','black',...
      'Foreground','red');
   uicontrol(...
      'Style','text',...
      'String','Conditions',...
      'Units','normalized',...
      'Position',[.83 .47 .11 .04],...
      'Background','black',...
      'Foreground','red');
   hBC1 = uicontrol(...
      'Style','radio',...
      'String','u=0',...
      'Units','normalized',...
      'Position',[left+.01 .43 btnWid .04],...
      'Background','black',...
      'Foreground','white',...
      'Value',1,...
      'CallBack','wvchar(''bc1'')');
   hBC2 = uicontrol(...
      'Style','radio',...
      'String','u''=0',...
      'Units','normalized',...
      'Position',[left+.01 .39 btnWid .04],...
      'Background','black',...
      'Foreground','white',...
      'CallBack','wvchar(''bc2'')'); 
   
   % The Restart button 
   
   uicontrol( ...
      'Style','push', ...
      'Units','normalized', ...
      'Position',[left .15  btnWid btnHt], ...
      'String','Restart', ...
      'Callback','wvchar(''restart'')');
   
   % The Close button 
   
   uicontrol('Style','Pushbutton', ...
      'Units','normalized',...
      'Position',[left .07 btnWid btnHt], ...
      'Callback','wvchar(''close'')',...
      'String','Close');
   
   % Fields to output starting points of characteristics and
   % the value of u(x,t) 
   
   hRoot1=uicontrol(...
      'Style','text',...
      'String','',...
      'Units','normalized',...
      'Background',bgcolor,...
      'Foreground',textcolor,...
      'Position',[.20 .16 .10 .04]);
   hRoot2=uicontrol(...
      'Style','text',...
      'String','',...
      'Units','normalized',...
      'Background',bgcolor,...
      'Foreground',textcolor,...
      'Position',[.55 .16 .10 .04]);
   
   uicontrol(...
      'Style','text',...
      'String','u(x,t) = ',...
      'Units','normalized',...
      'BackgroundColor',bgcolor,...
      'ForegroundColor',textcolor,... 
      'Position',[.25 .90 .12 .04]);
   hU=uicontrol(...
      'Style','text',...
      'String','',...
      'Units','normalized',...
      'Background',bgcolor,...
      'Foreground',textcolor,...
      'Position',[.37 .90 .10 .04]); 
   
   % Function field for initial condition u(x,0)
   
   uicontrol('Style','text',...
      'Units','normalized',...
      'Position',[.25 .02 .12 .05],...
      'BackgroundColor',bgcolor,...
      'ForegroundColor',textcolor,...
      'String','u(x,0) = ');
   
   f = 'x';
   hFunction = uicontrol(...
      'Style','edit',...
      'Units','normalized',...
      'Position',[.37 .02 .15 .05],...
      'BackgroundColor','white',...
      'String',f);
   
   % Set up initial plot are 
   
   hAxes = axes(...
      'Units','Normalized',...
      'Position',[.12 .28 .6 .6],...
      'XLim',[MinX MaxX],...
      'YLim',[MinT MaxT],...
      'Color','black',...
      'Box','on',...
      'ButtonDownFcn','wvchar(''click'')' ); 
   
   axis('square');
   hold on;
   
   ylabel('t');
   xlabel('x');
   
   % Divide xt-area into reference areas by plotting a background of
   % dashed characteristics. 
   
   Num_char = ceil( MaxT / ( (MaxX-MinX)/c ) );
   
   for n=1:Num_char,
      y1 = (n-1)* (MaxX - MinX)/c;
      y2 = y1 + (MaxX - MinX)/c;
      
      plot( [MinX MaxX],[y1 y2], 'b:', 'ButtonDownFcn','wvchar(''click'')')
      plot( [MinX MaxX],[y2 y1], 'b:', 'ButtonDownFcn','wvchar(''click'')')
   end;
   
   % Initial plot objects 
   
   hChar1 = plot([0], [0],'y','Erasemode','xor',...
      'ButtonDownFcn','wvchar(''click'')');
   hChar2 = plot([0], [0], 'y', 'Erasemode','xor',...
      'ButtonDownFcn','wvchar(''click'')');
   hClick = plot([0], [0], 'g*', 'Erasemode','xor', ...
      'ButtonDownFcn','wvchar(''click'')');
   drawnow;
   
   % store graphics handles for later reference
   
   handles=[hTmax hChar1 hChar2 hClick hBC1 hBC2 ...
         hRoot1 hRoot2 hFunction hU hSpeed hAxes];
      set(gcf,'UserData',handles); 
      
      % Refresh the plot area - redraw characteristics
      
   elseif strcmp(action,'restart'),
      
      MaxT = str2num( get(hTmax,'String') );
      
      set(hAxes,'XLim',[MinX MaxX]);
      set(hAxes,'YLim',[MinT MaxT]);
      
      cla;
      c = str2num(get(hSpeed,'string'));
      set(hSpeed,'UserData',c);
      Num_char = ceil( MaxT / ( (MaxX-MinX)/c ) );
      
      for n=1:Num_char,
         y1 = (n-1)* (MaxX - MinX)/c;
         y2 = y1 + (MaxX - MinX)/c;
         
         plot( [MinX MaxX],[y1 y2], 'b:', 'ButtonDownFcn','wvchar(''click'')')
         plot( [MinX MaxX],[y2 y1], 'b:', 'ButtonDownFcn','wvchar(''click'')')
      end;
      
      hChar1 = plot([0], [0], 'Erasemode','xor',...
         'ButtonDownFcn','wvchar(''click'')');
      hChar2 = plot([0], [0], 'Erasemode','xor',...
         'ButtonDownFcn','wvchar(''click'')');
      hClick = plot([0], [0], 'g*', 'Erasemode','xor', ...
         'ButtonDownFcn','wvchar(''click'')');
      drawnow;
      
      handles=[hTmax hChar1 hChar2 hClick hBC1 hBC2 ...
            hRoot1 hRoot2 hFunction hU hSpeed hAxes];
         set(gcf,'UserData',handles); 
         
         % Toggle the boundary conditions
         
      elseif strcmp(action(1:2),'bc'),
         
         if strcmp(action,'bc1'),
            set(hBC1,'Value',1)
            set(hBC2,'Value',0)
         else,
            set(hBC1,'Value',0)
            set(hBC2,'Value',1)
         end;
         
         % Draw characteristics from point of click to x-axis.  Compute
         % and display solution u(x,t)
         
      elseif strcmp(action,'click'),
         
         MaxT = str2num( get(hTmax,'String') );
         
         c = str2num(get(hSpeed,'string'))';
         c_old = get(hSpeed,'UserData');
         
         % If characteristic speed c has changed, update reference characteristics
         
         if ~(c == c_old),
            % clear background of reference characteristics 
            cla;        
            
            set(hAxes,'XLim',[MinX MaxX],'YLim',[MinT MaxT]);
            
            % redraw reference characteristics back to t=0
            
            Num_char = ceil( MaxT / ( (MaxX-MinX)/c ) );
            
            for n=1:Num_char,
               y1 = (n-1)* (MaxX - MinX)/c;
               y2 = y1 + (MaxX - MinX)/c;
               plot( [MinX MaxX],[y1 y2], 'b:','ButtonDownFcn','wvchar(''click'')')
               plot( [MinX MaxX],[y2 y1], 'b:','ButtonDownFcn','wvchar(''click'')')
            end;
            
            set(hSpeed,'UserData',c);    % redeclare char1,2 due to cla delete
            hChar1 = plot([0], [0], 'Erasemode','xor',...
               'ButtonDownFcn','wvchar(''click'')');
            hChar2 = plot([0], [0], 'Erasemode','xor',...
               'ButtonDownFcn','wvchar(''click'')');
            hClick = plot([0], [0], 'g*', 'Erasemode','xor',...
               'ButtonDownFcn','wvchar(''click'')');
            drawnow;
            
            handles=[hTmax hChar1 hChar2 hClick hBC1 hBC2 ...
                  hRoot1 hRoot2 hFunction hU hSpeed hAxes];
               set(gcf,'UserData',handles); 
            end;
            
            L = MaxX - MinX;
            
            % Get clicked point (x,t), extend characteristics from (x,t) back to
            % x-axis, compute solution u(x,t) based on initial values u(x,0).
            
            pt = get(gca,'currentpoint');
            x1 = pt(1,1);
            t1 = pt(1,2);
            
            set(hClick,'XData',x1,'YData',t1)
            
            %  follow characteristics starting left
            
            t(1) = t1;
            x(1) = x1;
            count = 1;
            
            while t(count) > 0,
               if x(count)==MinX, 
                  x(count+1) = MaxX;
                  t(count+1) = t(count) - L/c;
               elseif x(count)==MaxX,
                  x(count+1) = MinX;
                  t(count+1) = t(count) - L/c;
               else,
                  x(count+1) = MinX;
                  t(count+1) = MinX/c + t(count) - x(count)/c;
               end;
               count = count + 1;
            end;
            
            set(hChar1,'XData',x,'YData',t) 
            
            if x(count-1)==MaxX,
               y_int= t(count);
            elseif x(count-1)==MinX,
               y_int= t(count-1);
            else,
               y_int= -x(1)/c + t(1);
            end;     
            
            PorN=(-1)^count;
            slope=PorN*(1/c);
            p=[slope y_int];
            r=roots(p);
            
            t(1) = t1;
            x(1) = x1;
            count = 1;
            
            while t(count) > 0,
               if x(count)==MinX, 
                  x(count+1) = MaxX;
                  t(count+1) = t(count) - L/c;
               elseif x(count)==MaxX,
                  x(count+1) = MinX;
                  t(count+1) = t(count) - L/c;
               else,
                  x(count+1) = MaxX;
                  t(count+1) = -MaxX/c + t(count) + x(count)/c;
               end;
               count = count + 1;
            end;
            
            set(hChar2,'XData',x,'YData',t) 
            
            if x(count-1)==MaxX,
               y_int= t(count);
            elseif x(count-1)==MinX,
               y_int= t(count-1);
            else,
               y_int= x(1)/c + t(1);
            end;
            
            NorP=(-1)^count;
            slope=(-1)^(count+1)*(1/c);
            p=[slope y_int];
            r1=roots(p);
            
            f = get(hFunction,'String');
            x=r;
            y=eval(f);
            x=r1;
            y1=eval(f);   
            
            if get(hBC1,'Value')==1,
               u=.5*(y*PorN + y1*NorP);
               if (min(r,r1)==r),
                  set(hRoot1,'string',num2str(y*PorN));
                  set(hRoot2,'string',num2str(y1*NorP)); 
               else,
                  set(hRoot2,'string',num2str(y*PorN));
                  set(hRoot1,'string',num2str(y1*NorP)); 
               end;
            else,
               u=.5*(y + y1);
               if (min(r,r1)==r),
                  set(hRoot1,'string',num2str(y));
                  set(hRoot2,'string',num2str(y1)); 
               else,
                  set(hRoot2,'string',num2str(y));
                  set(hRoot1,'string',num2str(y1)); 
               end;
            end;
            
            set(hU,'string',num2str(u));
            
         elseif strcmp(action,'close'),
            
            close;
            
         end;