function fseries(action)

% FSERIES is a function which allows the user to observe the convergence of
%          Fourier series.  It is started with the command fseries at the
%          MATLAB prompt.  This causes the FSERIES Setup window to appear.
%          In this window the user can enter the function that is to be
%          approximated.  Since it is often the case that the function of
%          interest is discontinuous, it is possible to give different
%          definitions on up to four contiguous intervals.  However, it is
%          important to choose the number of intervals used first.
%
%          When the function is entered, clicking the Proceed button will
%          cause MATLAB to compute a large number of Fourier coefficients,
%          and to open the FSERIES Display window.  In this window the
%          function will be plotted, along with the Fourier approximation
%          of order 0.  The order of the approximation can be changed to
%          anything <= 1000 by changing the number in the edit box, and
%          hitting the Return or Enter keys.
%
%          For many purposes, it is useful to use the MATLAB command ZOOM
%          with FSERIES.  Use HELP ZOOM to get information about ZOOM.  the
%          click and drag option for ZOOM is probably the most useful in
%          this case.

%  Copywright (c) John C. Polking, Rice University
%  Last modified:  July 3, 1995

if nargin <1
   action ='initialize';
end


if strcmp(action,'initialize')
   
   % First we make sure that there is no other copy of FSERIES
   % running, since this causes problems.
   
   figs = findobj('tag','fseries');
   
   if (~isempty(figs))
      error('Only one copy of fseries can be open at one time.')
   end
   
   % Build the Setup window.
   
   N = 4096;		% Number of data points.
   M = 2;			% Default degree of approx.
   
   texth =20;   		% Height of text boxes.
   varw = 40;		% Length of variable boxes.
   equalw =20;		% Length of equals.
   eqlength = 250;		% Length of right hand sides.
   space = 30;
   nw = 60;
   left =10;		% Left margin.
   frw = 10+eqlength+space+2*nw+2*equalw+varw;
   
   NN = 2;
   
   separation = texth+3;	% Separation between boxes.
   buttbott = 10;
   
   frbott = separation +2*buttbott;
   eqbott = frbott + buttbott + (fliplr(0:3))*separation;
   dirbott = eqbott(1)+separation;
   labelbott = dirbott + separation;
   frht = 2*buttbott + 6*separation;
   
   defigwidth =2*left + frw;	% Width of the figure.
   defigureheight = frbott + frht + buttbott;		% Height of the figure.
   defigurebot = 40;				% Bottom of the figure.
   setpos = [20 defigurebot defigwidth defigureheight];
   
   
   fourset = figure('pos',setpos,...
      'resize','off','name','FSERIES Setup','numb','off',...
      'visible','off','tag','fseries');
   
   
   frame = uicontrol('style','frame',...
      'pos',[left,frbott,frw,frht],...
      'visible','off');
   
   label=uicontrol('style','text',...
      'pos',[left+5 labelbott frw-10 texth],...
      'horizon','left',...
      'string','The function may be defined differently on up to four contiguous intervals.',...
      'visible','off');
   
   dirlength = eqlength;
   directions=uicontrol('style','text',...
      'pos',[left+5 dirbott dirlength texth],...
      'horizon','left',...
      'string','Choose the number of intervals.','visible','off');
   
   radleft =left + dirlength + 10;
   radwidth = 40;
   
   rad(1) = uicontrol('style','radio',...
       'pos',[radleft dirbott radwidth texth],...
      'string','1','value',0,'visible','off');
   
   rad(2) = uicontrol('style','radio',...
       'pos',[radleft+radwidth dirbott radwidth texth],...
      'string','2','value',2,'max',2,'visible','off');
   
   rad(3) = uicontrol('style','radio',...
       'pos',[radleft+2*radwidth dirbott radwidth texth],...
      'string','3','value',0,'max',3,'visible','off');
   
   rad(4) = uicontrol('style','radio',...
       'pos',[radleft+3*radwidth dirbott radwidth texth],...
      'string','4','value',0,'max',4,'visible','off');
   
   for i=1:4
      set(rad(i),'user',[rad(:,[1:(i-1),(i+1):4])]);
   end
   
   set(rad,'call','fseries(''callrad'')');
   
   
   Xname = str2mat(' x + 1 ', ' x - 1 ','','');
   lefts = str2mat(' -1 ',' 0 ','','');
   rights = str2mat(' 0 ',' 1 ','','');
   
   eqcall = '';
   leftcall = '';
   varcall = '';
   rightcall = '';
   
   forl = left +5 + eqlength;
   leftl = forl + space;
   leftr = leftl + nw;
   varl = leftr +equalw;
   rightl = varl + varw;
   rightr = rightl +equalw;
   
   for kk = 1:4
      eqn(kk) = uicontrol('pos',[left+5 eqbott(kk) eqlength texth],...
         'style','edit','tag',int2str(kk),...
         'horizon','left','string',deblank(Xname(kk,:)),...
         'call','fseries(''eqcall'')','visible','off');
      
      forr(kk) = uicontrol('pos',[forl eqbott(kk) space texth],...
         'style','text','tag',int2str(kk),...
         'horizon','center','string','for',...
         'visible','off');
      
      leftw(kk) = uicontrol('pos',[leftl eqbott(kk) nw texth],...
         'style','edit','tag',int2str(kk),...
         'horizon','center','string',deblank(lefts(kk,:)),...
         'call','fseries(''leftcall'')','visible','off');
      
      lefteq(kk) = uicontrol('pos',[leftr eqbott(kk) equalw texth],...
         'style','text','tag',int2str(kk),...
         'horizon','center','string','<',...
         'visible','off');
      
      if kk ==1
         varr(kk) = uicontrol('pos',[varl eqbott(kk) varw texth],...
            'style','edit','tag',int2str(kk),...
            'horizon','center','string',' x ',...
            'call','fseries(''varcall'')','visible','off');
      else
         varr(kk) = uicontrol('pos',[varl eqbott(kk) varw texth],...
            'style','text','tag',int2str(kk),...
            'horizon','left','string',' x ',...
            'visible','off');
      end
      
      righteq(kk) = uicontrol('pos',[rightl eqbott(kk) equalw texth],...
         'style','text','tag',int2str(kk),...
         'horizon','center','string','<',...
         'visible','off');
      
      rightw(kk) = uicontrol('pos',[rightr eqbott(kk) nw texth],...
         'style','edit','tag',int2str(kk),...
         'horizon','center','string',deblank(rights(kk,:)),...
         'call','fseries(''rightcall'')','visible','off');
      
      
   end
   
   butt(1) = uicontrol('style','push',...
      'pos',[(defigwidth/4-35),buttbott,70,texth],...
      'string','Quit',...
      'call','fseries(''quit'')','visible','off');
   
   butt(2) = 0;
   
   butt(3) = uicontrol('style','push','enable','on',...
      'pos',[(3*defigwidth/4-35),buttbott,70,texth],...
      'string','Proceed',...
      'call','fseries(''proceed'')','visible','off');
   
   butt(4) = 0;
     
   ud = [eqn,leftw,rightw,varr,butt,frame,NN,separation,N,0];
   
   set(fourset,'resize','on');
   hh = get(fourset,'child');
   set(hh,'units','normal');
   global FSSETPOS
   
   if all(size(FSSETPOS) == [1,4])
      set(fourset,'pos',FSSETPOS);
   end
   
   
   set(fourset,'user',ud);
   set(fourset,'vis','on');
   set(hh,'vis','on');
   for kk = NN+1:4
      set(findobj(hh,'tag',int2str(kk)),'vis','off');
   end	
   
elseif strcmp(action,'proceed')
 
   ud = get(gcf,'user');
   eqn = ud(1:4);
   leftw = ud(5:8);
   rightw = ud(9:12);
   var = ud(13);
   var = get(var,'string');
   butt = ud(17:20);
   NN = ud(22);
   N = ud(24);
   fourdisp = ud(25);
   fourset = gcf;
   M = 0;
   mmmm = int2str(M);
   
   a = [];
   b = [];
   for kk = 1:NN
      aaa = str2num(get(leftw(kk),'string'));
      bbb = str2num(get(rightw(kk),'string'));
      ffcn = get(eqn(kk),'string');
      l = length(ffcn);
      if ((l==0)|(isempty(aaa))|(isempty(bbb)))
         fseries('alert');
         return;
      end
      a = [a,aaa];
      b = [b,bbb];
      eval([var,'=aaa;']);
      A(kk) = eval(ffcn);
      eval([var,'=bbb;']);
      B(kk) = eval(ffcn);
      
      for (k=fliplr(find((ffcn=='^')|(ffcn=='*')|(ffcn=='/'))))
         ffcn = [ffcn(1:k-1) '.' ffcn(k:l)];
         l = l+1;
      end
      if kk ==1
         fffcn = ffcn;
      else
         fffcn = str2mat(fffcn,ffcn);
      end
   end
   mmerr = abs(B(NN) - A(1));
   for kk = 1:NN-1
      mmerr = max(mmerr, abs(B(kk)-A(kk+1)));
   end
   mmerr = mmerr/2;
   aa = a(1);
   bb = b(NN);
   xx = aa +(0:N)*(bb-aa)/N;
   wf = NaN*ones(1,N+1);
   fm = wf;
   F1 = deblank(fffcn(1,:));
   FN = deblank(fffcn(NN,:));
   eval([var,'=xx(1);']);
   %	x = xx(1);
   fm(1) = eval(F1);
   eval([var,'=xx(N+1);']);
   %	x = xx(N+1);
   fm(1) = fm(1) + eval(FN);
   fm(1) = fm(1)/2;
   fm(N+1) = fm(1);
   FF = '';
   for kk = 1:NN
      kkk = find((a(kk)<xx)&(xx<b(kk)));
      FL = FF;
      FF = deblank(fffcn(kk,:));
      eval([var,'=xx(kkk);']);
      %		x = xx(kkk);
      wff= eval(FF);
      if length(wff) == 1
         wf(kkk) = wff*ones(size(kkk));
      else
         wf(kkk) = wff;
      end
      fm(kkk) = wf(kkk);
      kmin = min(kkk);
      if isnan(fm(kmin - 1))
         eval([var,'=xx(kmin-1);']);
         %			x = xx(kmin - 1);
         fm(kmin - 1) = (eval(FL) + eval(FF))/2;
      end
      if kmin > 2
         wf(kmin-1) = NaN;
      end
   end
   
   if (fourdisp == 0)
      
      ss = get(0,'screensize');
      sw = ss(3);
      setpos = get(fourset,'pos');
      setl = setpos(3);
      setw = setpos(3);
      
      displeft = min(setl+setw, sw -520); 
      dispbott = 20; 
      dispwidth = 550; 
      dispht = 392;
      fourdisp = figure('pos',[displeft, dispbott, dispwidth, dispht],...
         'resize','off','name','FSERIES Display','numb','off',...
         'visible','off','tag','fseries');
      ud(25) = fourdisp;
      set(fourset,'user',ud);
      
      
      texth =20;   		% Height of text boxes.
      space = 40;
      buttbott = 10;
      
      daxbott = 80/dispht;
      dxw = 450;
      daxw = dxw/dispwidth;
      dxleft = (dispwidth - dxw)/2;
      daxleft = dxleft/dispwidth;
      daxht = 281.25/dispht;
      diax = axes('pos',[daxleft, daxbott, daxw, daxht],'vis','off');
      
      frame = uicontrol('pos',[dxleft,buttbott,dxw,texth+10],...
         'style','frame');
      
      degree = uicontrol('pos',[dxleft+5 buttbott+5 space texth],...
         'style','edit','tag','deg',...
         'horizon','center','string',mmmm,...
         'call','fseries(''approx'')','visible','off');
      
      degrlab = 'Enter the desired degree of approximation.  ( <= 1000 )';
      deglabw = dxw - space -20 -5;
      deglabel = uicontrol('pos',[dxleft+space+20 buttbott+5 deglabw texth],...
         'style','text',...
         'horizon','left','string',degrlab,...
         'visible','off');
      hh = get(fourdisp,'child');
      set(findobj(fourdisp,'type','axes'),'units','normal');
      set(findobj(fourdisp,'type','uicontrol'),'units','normal');
      set(fourdisp,'resize','on');
      global FSDISPPOS
      if all(size(FSDISPPOS) == [1,4])
         set(fourdisp,'pos',FSDISPPOS);
      end
      
      set(fourdisp,'vis','on');
      set(hh,'vis','on');
      grid
      set(diax,'vis','on'); %,'xgrid','on','ygrid','on');
      
   else
      figure(fourdisp);
      hh = get(fourdisp,'user');
      degree = hh(4);
      set(degree,'string',mmmm);
      diax = gca;
   end
   
   aa = min(xx);  bb = max(xx);
   cc = min(wf(finite(wf)));
   dd = max(wf(finite(wf)));
   ee = 0.05*(max((dd-cc),1));
   cc = cc-ee;
   dd = dd+ee;
   
   hact = plot(xx,wf);
   axis([aa,bb,cc,dd]);
   set(hact,'user',fm);
   v = fft(fm(1:N));
   set(diax,'user',v);
   
   w = v;
   w(M+2:N-M) = zeros(1,N-2*M-1);
   sfm = zeros(1,N+1);
   sfm(1:N) = real(ifft(w));
   sfm(N+1) = sfm(1);
   cc = min(cc,min(sfm));
   dd = max(dd,max(sfm));
   hold on 
   happ = plot(xx,sfm,'--r');
   hold off
   axis([aa,bb,cc,dd]);
   set(fourdisp,'user',[hact,happ,mmerr,degree]);
   s2=['The Fourier series approximation of degree ',mmmm,'.'];
   title(s2);
   ferr=abs(fm-sfm);
   maxerr=max(max(ferr),mmerr);
   step = (xx(N+1)-xx(1))/N;
   L2err=norm(ferr)*sqrt(step);
   l2e=num2str(L2err);
   merr=num2str(maxerr);
   s3=['The maximum error is ',merr,'.    The L2 error is ',l2e,'.'];
   xlabel(s3);
   hh = get(fourdisp,'child');
   %set(hh,'vis','on');
      
   
elseif strcmp(action,'eqcall')
   
   ud = get(gcf,'user');
   butt = ud(17:20);
   set(butt(3),'enable','on');
   
elseif strcmp(action,'leftcall')
   
   ud = get(gcf,'user');
   butt = ud(17:20);
   rightw = ud(9:12);
   set(butt(3),'enable','on');
   numb = get(gco,'tag');
   N = str2num(numb);
   if N>1
      bb = get(gco,'string');
      nn = int2str(N-1);
      hh = findobj(rightw,'tag',nn);
      set(hh,'string',bb);
   end
   
elseif strcmp(action,'rightcall')
   
   ud = get(gcf,'user');
   butt = ud(17:20);
   leftw = ud(5:8);
   set(butt(3),'enable','on');
   numb = get(gco,'tag');
   N = str2num(numb);
   if N<4
      bb = get(gco,'string');
      nn = int2str(N+1);
      hh = findobj(leftw,'tag',nn);
      set(hh,'string',bb);
   end
   
elseif strcmp(action,'varcall')
   
   ud = get(gcf,'user');
   butt = ud(17:20);
   varr = ud(13:16);
   set(butt(3),'enable','on');
   xx = get(gco,'string');
   set(varr,'string',xx);
   
   
   
elseif strcmp(action,'approx')
   
   mmax = 1000;
   hm = gco;
   v=get(gca,'user');
   hh = get(gcf,'user');
   hact = hh(1);
   happ = hh(2);
   mmerr = hh(3);
   xx = get(hact,'xdata');
   fm = get(hact,'user');
   mmmm = get(hm,'string');
   M = str2num(mmmm);
   M = max(0,min(M,mmax));
   mmmm = num2str(M);
   set(hm,'string',mmmm);
   w = v;
   N = length(xx) - 1;
   w(M+2:N-M) = zeros(1,N-2*M-1);
   sfm = zeros(1,N+1);
   sfm(1:N) = real(ifft(w));
   sfm(N+1) = sfm(1);
   wf = get(hact,'ydata');
   cc = min(wf(finite(wf)));
   dd = max(wf(finite(wf)));
   ee = 0.05* max((dd-cc),1);
   dd = max(max(sfm),dd)+ee;
   cc = min(min(sfm),cc)-ee;
   set(gca,'ylim',[cc,dd]);
   
   hold on 
   set(happ,'ydata',sfm);
   hold off
   s2=['The Fourier series approximation of order ',mmmm,'.'];
   title(s2);
   ferr=abs(fm-sfm);
   maxerr=max(max(ferr),mmerr);
   step = (xx(N+1)-xx(1))/N;
   L2err=norm(ferr)*sqrt(step);
   l2e=num2str(L2err);
   merr=num2str(maxerr);
   s3=['The maximum error is ',merr,'.    The L2 error is ',l2e,'.'];
   xlabel(s3);
   outdata =[xx;wf;sfm];
   assignin('base','fsdata',outdata);
   
elseif strcmp(action,'callrad')
   
   me = get(gcf,'currentobject');
   vv = get(me,'value');
   mm = get(me,'max');
   if vv 
      hand = get(me,'user');
      set(hand(1:3),'value',0);
   end
   set(me,'value',mm);
   us = get(gcf,'user');
   NN = us(22);
   if mm < NN
      us(22)=mm;
      set(gcf,'user',us);
      hh = get(gcf,'child');
      %		set(hh,'vis','on');
      for kk = mm+1:NN
         set(findobj(hh,'tag',int2str(kk)),'vis','off');
      end
   else
      us(22)=mm;
      set(gcf,'user',us);
      hh = get(gcf,'child');
      %		set(hh,'vis','on');
      for kk = NN+1:mm
         set(findobj(hh,'tag',int2str(kk)),'vis','on');
      end
   end
   
elseif strcmp(action,'alert')
   
   ud = get(gcf,'user');
   NN = ud(22);
   astr = ['All of the data for the ' num2str(NN) ' intervals must be provided.'];
errh = errordlg(astr,'FSERIES data entry error');

   
elseif strcmp(action,'quit')
  
   h = findobj('tag','fseries');
   close(h);
end