function area=numsimpson(option)
%NUMINT  Shows Simpson's rule for a function
%USE:    Type  numsimpson  at the MATLAB prompt

%Peter Dunn
%26 April 1999

if nargin==0,

	figure('Name','Simpson''s Rule', 'tag','tagSimpson',...
   		 'NumberTitle','off');

   enterframe = uicontrol('Style','Frame','Units','normalized',...
            'Position',[0.02 0.02 0.79, 0.24]);

   funtext = uicontrol('Style','text','Units','normalized',...
             'FontSize',14,...
             'Position',[0.03 0.12 0.16,0.08],'String','Function:',...
             'HorizontalAlignment','center');
   fun = uicontrol('Style','edit','Units','normalized',...
         'FontSize',14,...
         'ToolTip','The function to integrate',...
         'Position',[0.19 0.15 0.47 0.08],...
         'Callback','numsimpson(5);',...
         'String','sin(x.^2)','tag','tagSimpsonFun' );

   lowertext = uicontrol('Style','text','Units','normalized',...
         'FontSize',14,...
         'Position',[0.15 0.15 0.47 0.08],...
         'Position',[0.03 0.03 0.11 0.08],'String','Min x:' );
   lower = uicontrol('Style','edit','Units','normalized',...
         'FontSize',14,...
         'ToolTip','The lower integration limit',...
         'Position',[0.15 0.03 0.08 0.08],...
         'Callback','numsimpson(6);',...
         'String','0','tag','tagSimpsonLower');

   uppertext = uicontrol('Style','text','Units','normalized',...
         'FontSize',14,...
         'Position',[0.24 0.03 0.11 0.08],'String','Max x:');
   upper = uicontrol('Style','edit','Units','normalized',...
         'FontSize',14,...
         'ToolTip','The upper integration limit',...
         'Position',[0.36 0.03 0.08 0.08],...
         'Callback','numsimpson(6);',...
         'String','2','tag','tagSimpsonUpper');

   numptstext = uicontrol('Style','text','Units','normalized',...
         'FontSize',14,...
         'Position',[0.45 0.03 0.12 0.08],'String','Panels:');
   numpts = uicontrol('Style','edit','Units','normalized',...
         'FontSize',14,...
         'ToolTip','Number of panels to use',...
         'Callback','numsimpson(10);',...
         'Position',[0.58 0.03 0.08 0.08],'String','6',...
         'tag','tagSimpsonNum');

   speedtextup = uicontrol('Style','text','Units','normalized',...
         'FontSize',14,...
         'Position',[0.70 0.20 0.10 0.04],'String','50');
   speedtextlo = uicontrol('Style','text','Units','normalized',...
         'FontSize',14,...
         'ToolTip','Changes the number of panels to use',...
         'Position',[0.70 0.03 0.10 0.04],'String','2');
   speed = uicontrol('Style','slider','Units','normalized',...
         'Position',[0.72 0.08 0.06 0.12],'Min',2,'Max',50,...
         'Callback','numsimpson(11);',...
         'SliderStep',[2/48 4/48],  'Value',6,'tag','tagSimpsonSlide');

   gobtn = uicontrol('Style','pushbutton','Units','normalized',...
             'Position',[0.83 0.14 0.12 0.12],...
             'FontSize',14,...
             'ToolTip','Approximates the integral',...
             'String','GO','Callback','numsimpson(4); ');
   quitbtn = uicontrol('Style','pushbutton','Units','normalized',...
             'Position',[0.83 0.02 0.12 0.12],...
             'FontSize',14,...
             'ToolTip','Closes the program',...
             'String','QUIT','Callback','delete(gcf);');

   uimenu('Label','DEMONSTRATIONS', 'tag','tagSimpsonDemo');
   uimenu( findobj('tag', 'tagSimpsonDemo'), ...
           'Callback','numsimpson(20);',...
           'Label','Difficult Regions: sqrt(x)');
   uimenu( findobj('tag', 'tagSimpsonDemo'), ...
           'Callback','numsimpson(21);',...
           'Label','Difficult Regions: exp( -abs(x) )');
   uimenu( findobj('tag', 'tagSimpsonDemo'), ...
           'Callback','numsimpson(22);',...
           'Label','Difficult Regions: 1 ./ ( 1 + x.^6 )');
   uimenu( findobj('tag', 'tagSimpsonDemo'), ...
           'Callback','numsimpson(23);',...
           'Label','Difficult Regions: besselj(0.2, x)');
   uimenu( findobj('tag', 'tagSimpsonDemo'), ...
           'Callback','numsimpson(24);',...
           'Label','Difficult Regions: sin( 1 ./ x )');
   uimenu( findobj('tag', 'tagSimpsonDemo'), ...
           'Callback','numsimpson(25);',...
           'Label','Exact Answer: x.^3 - x - 1');
   uimenu( findobj('tag', 'tagSimpsonDemo'), ...
           'Callback','numsimpson(26);',...
           'Label','Exact Answer! ( sin(x) )^2');

   axes('Position',[0.10 0.35 0.85, 0.6]);
   axis( [ 0 2 0 1 ] );
   
   set( gca,'FontSize',16);

   zoom;

   %Flag that nothing has been printed yet:
   set( findobj('tag','tagSimpson'), 'UserData',0);

   %numsimpson(4);


elseif nargin==1,

   if option == 3, %NUMBER REGIONS

      np = get( findobj(gcf,'tag','tagSimpsonNum'),'String');
      np = floor(str2num(np));

      if np<1,
         watchoff;
         herr=errordlg('Need at least two points!','Error!');
         set(herr,'WindowStyle','modal');
         return;
      end;

      if np>50,
         watchoff;
         herr=errordlg(...
               'More than 50 points will slow MATLAB down too much!',...
               'Error!');
         set(herr,'WindowStyle','modal');
         set( findobj(gcf,'tag','tagSimpsonNum'),'String','50');
         np = 50;
         return;
      end;

      if ( floor(np/2) ~= (np/2) ),
         watchoff;
         herr=errordlg(...
               'Need an *even* number of points!',...
               'Error!');
         set(herr,'WindowStyle','modal');
         return;
      end;
  

   elseif option == 4, %FUNCTION

      prt = get( findobj('tag','tagSimpson'),'UserData');
      set( findobj('tag','tagSimpson'),'UserData',1);

      if ( prt==0 ),

         %Set up Command Window output
         numsimpson(5);

      end;

      watchon;

      delete(gca);

      numsimpson(3);

      f = get( findobj(gcf,'tag','tagSimpsonFun'),'String');
      
      xlow = get( findobj(gcf,'tag','tagSimpsonLower'), 'String');
      xupp = get( findobj(gcf,'tag','tagSimpsonUpper'), 'String');
      numpt = get( findobj(gcf,'tag','tagSimpsonNum'), 'String');

      %Plot the function
      xl = str2num(xlow);
      xu = str2num(xupp);
      np = str2num(numpt);

      axes('Position',[0.10 0.35 0.85, 0.6]);

      x = linspace(xl,xu,1000);

      prob = 0;
      eval(['fy=',f,';'],'fy=NaN;');
      fx = x;

      if any(isnan(fy)),
         watchoff;
         herr=errordlg(...
              ['Couldn''t evaluate the function!  Check the function ',...
              'and the evaluation limits.'],'Error evaluating function');
         set(herr,'WindowStyle','modal');
         return;
      end;

      hold on;

      plot( fx, fy, 'g','LineWidth',4 );
      plot( [min(fx),max(fx)],[0 0],'k-');
      hold on;

      %Determine interval info
      h = ( xu - xl ) / np ;
      xx = [ xl : h : xu ];
      x = xx;

      approxarea = 0;

      for i = 1:2:np-1,

         x1 = x(i);
         x2 = x(i+1);
         x3 = x(i+2);
         x = x1;
         y1 = eval(f);
         x = x2;
         y2 = eval(f);
         x = x3;
         y3 = eval(f);
         x = xx;

         [area, hndl] = plotquad( x1, x2, x3, y1, y2, y3, f );
         approxarea = approxarea + area;

      end;

      %LEGEND is slow and causes the figure to be redrawn... speed up?
      legend( [hndl],['Integral = ',num2str(approxarea)],0);

      set( gca,'FontSize',16);

      fprintf( '%5i   %15.7f\n', np, approxarea);

      watchoff;

   elseif option==5, %Header for output

      %Set up Command Window output

      f = get( findobj('tag','tagSimpsonFun'),'String');
      lo = get( findobj('tag','tagSimpsonLower'),'String');
      hi = get( findobj('tag','tagSimpsonUpper'),'String');
      disp(' ');
      disp( [ 'Integrating ',f,' between ', ...
              num2str(lo),' and ',num2str(hi) ] );
      disp(' ');
      disp('Number of    Approximate');
      disp(' Panels        Integral');
      disp('~~~~~~~~~    ~~~~~~~~~~~~');

   elseif option==6,

      %Flag need for new header

      set( findobj('tag','tagSimpson'), 'UserData',0);

   elseif option==10, %Number Panels changed

      val = get( findobj('tag','tagSimpsonNum'), 'String');
      %Convert to number
      val = val{1};
      val = eval(val);
      val = floor(val);

      if val > 50,
         val = 50;
         set( findobj('tag','tagSimpsonNum'), 'String','50');
      elseif val <= 0,
         val = 2;
         set( findobj('tag','tagSimpsonNum'), 'String','2');
      else
         set( findobj('tag','tagSimpsonNum'), 'String',num2str(val));
      end;

      set( findobj('tag','tagSimpsonSlide'), 'Value', val );

   elseif option==11, %Slider moves

      val = get( findobj('tag','tagSimpsonSlide'), 'Value' );
      val = floor(val);
      %Make even
      if ( (val/2) ~= floor(val/2) ),
         val = val+1;
      end
      set( findobj('tag','tagSimpsonNum'), 'String',num2str(val) );

      numsimpson(4);

   elseif option==20,

      numsimpson(6);

      set( findobj('tag','tagSimpsonFun'), ...
           'String','sqrt( x )' );
      set( findobj('tag','tagSimpsonLower'),...
           'String','0');
      set( findobj('tag','tagSimpsonUpper'),...
           'String','4');
      set( findobj('tag','tagSimpsonSlide'),...
           'Value',2);
      set( findobj('tag','tagSimpsonNum'),...
           'String','2');

      numsimpson(4);

   elseif option==21,

      numsimpson(6);

      set( findobj('tag','tagSimpsonFun'), ...
           'String','exp( -abs(x) )' );
      set( findobj('tag','tagSimpsonLower'),...
           'String','-2');
      set( findobj('tag','tagSimpsonUpper'),...
           'String','2');
      set( findobj('tag','tagSimpsonSlide'),...
           'Value',2);
      set( findobj('tag','tagSimpsonNum'),...
           'String','2');

      numsimpson(4);

   elseif option==22,

      numsimpson(6);

      set( findobj('tag','tagSimpsonFun'), ...
           'String','1 ./ ( 1 + x.^6 )' );
      set( findobj('tag','tagSimpsonLower'),...
           'String','-3');
      set( findobj('tag','tagSimpsonUpper'),...
           'String','4');
      set( findobj('tag','tagSimpsonSlide'),...
           'Value',2);
      set( findobj('tag','tagSimpsonNum'),...
           'String','2');

      numsimpson(4);

   elseif option==23,

      numsimpson(6);

      set( findobj('tag','tagSimpsonFun'), ...
           'String','besselj( 0.2, x )' );
      set( findobj('tag','tagSimpsonLower'),...
           'String','0');
      set( findobj('tag','tagSimpsonUpper'),...
           'String','1');
      set( findobj('tag','tagSimpsonSlide'),...
           'Value',2);
      set( findobj('tag','tagSimpsonNum'),...
           'String','2');

      numsimpson(4);

   elseif option==24,

      numsimpson(6);

      set( findobj('tag','tagSimpsonFun'), ...
           'String','sin( 1 ./ x )' );
      set( findobj('tag','tagSimpsonLower'),...
           'String','0.1');
      set( findobj('tag','tagSimpsonUpper'),...
           'String','1');
      set( findobj('tag','tagSimpsonSlide'),...
           'Value',2);
      set( findobj('tag','tagSimpsonNum'),...
           'String','2');

      numsimpson(4);

   elseif option==25,

      numsimpson(6);

      set( findobj('tag','tagSimpsonFun'), ...
           'String','x.^3 - x - 1' );
      set( findobj('tag','tagSimpsonLower'),...
           'String','-1');
      set( findobj('tag','tagSimpsonUpper'),...
           'String','2');
      set( findobj('tag','tagSimpsonSlide'),...
           'Value',2);
      set( findobj('tag','tagSimpsonNum'),...
           'String','2');

      numsimpson(4);

   elseif option == 26,

      % Thanks to Richard Oakey (Auckland University of Technology)
      % for this demonstration
      
     	numsimpson(6);

      set( findobj('tag','tagSimpsonFun'), ...
           'String','sin(x) .^ 2' );
      set( findobj('tag','tagSimpsonLower'),...
           'String','0');
      set( findobj('tag','tagSimpsonUpper'),...
           'String','2*pi');
      set( findobj('tag','tagSimpsonSlide'),...
           'Value',2);
      set( findobj('tag','tagSimpsonNum'),...
           'String','4');

      numsimpson(4);
   end;

end;

return;


%%%%%%%%%%%%%%%%
%  SUBFUNCTION %
%%%%%%%%%%%%%%%%

function [area, h]=plotquad(x1, x2, x3, y1, y2, y3, f)
%PLOTS the parabola

%Find parabola parameters
X = [ x1.^2  x1  1 ; ...
      x2.^2  x2  1 ; ...
      x3.^2  x3  1 ];

pars = inv(X) * [y1; y2; y3] ;
%pars = [y1; y2; y3] \ X;
a = pars(1);
b = pars(2);
c = pars(3);

x = linspace( x1, x3, 20);
f = get( findobj(gcf,'tag','tagSimpsonFun'),'String');

suby = a * (x.^2) + b * x + c;
%suby = evalin('base',f);
%suby = eval(f);

h = plot( [x1 x x3], [0 suby 0] , 'r-',...
   'LineWidth', 2 );
plot( x1, y1, 'r.' , 'MarkerSize', 20);
plot( x2, y2, 'r.' , 'MarkerSize', 20);
plot( x3, y3, 'r.' , 'MarkerSize', 20);
%h = plot( [xlow xlow xupp xupp], [0 ht htr 0] , 'r-' );

fB = (a*x3^3)/3 + (b*x3^2)/2 + c*x3;
fA = (a*x1^3)/3 + (b*x1^2)/2 + c*x1;
area = fB - fA;

set( gca, 'FontSize', 14);

return;