function Interpolate(x,y)
%
%INPUT ==>() or (x,y) Call with nothing or a discrete function defined by x and y
%OUTPUT==>Draws graph and lets you adds points and interpolate.
%
%You may declare a function first by typing x=-1:.001:1; and y=sin(5*x);
%on the command line first.  Or just run this program with no arguments.
%If you run it with no arguments, then after it starts, click the 
%button that says "New Fcn" to draw a function.
%After drawing a function either close that window or hit the "x" key.
%This returns you to the interpolation interface.  Click on your drawn
%function to add interpolation points.  When you are ready to see an
%interpolation, hit the button labeled "Poly".  To clear away a single
%point, click near it while holding down the 'ctrl' key.  To add a point
%not on the graph, hold down the shift key.  If you want your points to be
%the roots of Chebyshev polynomials, first fill in a number in bottom
%middle and then click "+Cheby".  "+Even" adds evenly spaced points.
%"PieceWs" interpolates with a linear piece wise function.  If you make the
%window full screen you will see all the the keyboard keys defined on the
%bottom of the window.  If you have any questions email me at
%cdeotte@ucsd.edu.  Enjoy.  I hope it helps you understand math better.
%

try %if no x,y was inputed
   y;
catch
    x=-1:.001:1;
    y=0*x;
    fprintf('\nYou can call Lagrange(x,y) after first declaring a discrete function.\n');
    fprintf('i.e. x=-1:.001:1; and then y=sin(5*x); or something.\n\n');
end;
hide2=0;
poly=0;
piece=0;
points=0;
index=0;
n=0;
t1=0;
t2=0;
tol=.0001;
init(1,1,0);

function init(erase_points,set_axis,hide2);
  if set_axis==0
      axis2=axis;
  end;
  oldvalue='empty';
  try
    oldvalue=get(t2,'String');
  end;
  clf;
  n=max(size(x));
  if erase_points
    points=zeros(1,3);
    index=0;
  end;
  if ~hide2
    plot(x,y,'HitTest','off');
  end;
  set(gca,'ButtonDownFcn',@down);
  set(gcf,'KeyPressFcn',@push4);
  ymx=max(y);
  ymn=min(y);
  ydf=ymx-ymn;
  if ydf==0 ydf=1; end;
  xdf=x(n)-x(1);
  if set_axis
    axis([x(1)-xdf/10 x(n)+xdf/10 ymn-.5*ydf ymx+.5*ydf]);
  else
    axis(axis2);
  end;
  b1=uicontrol('Style','PushButton','String','Poly','Position',[15 30 50 20],'CallBack',@push);
  b2=uicontrol('Style','PushButton','String','Clear','Position',[15 70 50 20],'CallBack',@push2);
  b3=uicontrol('Style','PushButton','String','New Fcn','Position',[15 90 50 20],'CallBack',{@push5,0});
  b4=uicontrol('Style','PushButton','String','PieceWs','Position',[15 50 50 20],'CallBack',@push6);
  t1=uicontrol('Style','Text','String','N=0, Error=','Position',[15 5 150 15]);
  t2=uicontrol('Style','Edit','String','2','Position',[288 5 50 15]);
  b3=uicontrol('Style','PushButton','String','+Cheby','Position',[200 5 75 15],'CallBack',{@push3,0});
  b4=uicontrol('Style','PushButton','String','+Even','Position',[350 5 75 15],'CallBack',{@push3,1});
  t3=uicontrol('Style','Text','String','Keys: p,q:poly; i/o:zoom +-; c:clr; h:hide; b:cheby; e:even. Mouse: click:+pnt; clk&ctrl:-pnt; click&shift:any pnt','Position',[430,5,550,15]);
  if ~strcmp(oldvalue,'empty')
      set(t2,'String',oldvalue);
  end;
  hold on;
end;
    
    function DrawPoints();
        for i=1:index
            plot(points(i,1),points(i,2),'.r');
        end;
    end;

    function down(src,event); %Mouse button was clicked
       X=get(gca,'CurrentPoint');
       p1=[X(1,1) X(1,2)];
       remove=0;
       %disp(get(gcf,'SelectionType'));
       if strcmp(get(gcf,'SelectionType'),'alt') && index>0 
         remove=1;
         min=1000;
         minindex=0;
         for i=1:index
             m=abs(points(i,1)-p1(1));
             if m<min
                 min=m;
                 minindex=i;
             end;
         end;
         plot(points(minindex,1),points(minindex,2),'.w');
         points(minindex,:)=points(index,:);
         points=points(1:index-1,:);
         index=index-1;
       end;
                 
       if ~remove
         p2=[0 0 0];
         i=1;
         while i~=0
            if x(i)>=p1(1) || i==n
                p2=[x(i) y(i) i];
                if strcmp(get(gcf,'SelectionType'),'extend')
                    p2=[x(i) p1(2) i];
                end;
                i=0;
            else
                i=i+1;
            end;
         end;
         duplicate=0;
         for i=1:index
             if points(index,1)==p2(1)
                 duplicate=1;
             end;
         end;
         if ~duplicate
           plot(p2(1),p2(2),'.r');
           index=index+1;
           points(index,:)=p2;
           set(t1,'String',strcat(strcat('N=',num2str(index)),' Error=NA'));
         end;       
       end;
    end;
    
    function push4(src,event) %A key was hit.
        try    
        if event.Character=='b'
            push3(0,0,0);
        end;
        if event.Character=='e'
            push3(0,0,1);
            return;
        end;
        if event.Character=='w'
            push6(0,0);
            return;
        end;
        if event.Character=='h'
            hide2=~hide2;
            if poly push(0,0);end;
            if piece push6(0,0);end;
            return;
        end;
        if event.Character=='c'
            init(1,1,0);
            return;
        end;
        if event.Character=='n'
            push5(0,0,0);
            return;
        end;
        if event.Character=='m'
            push5(0,0,1);
            return;
        end;
        if event.Character=='p' || event.Character=='q'
            push(0,0);
            return
        end;
        if event.Character=='o'
            axis2=axis;
            ydff=axis2(4)-axis2(3);
            axis2=axis2+[0 0 -ydff/2 ydff/2];
            axis(axis2);
        end;
        if event.Character=='i'
            axis2=axis;
            ydff=axis2(4)-axis2(3);
            axis2=axis2+[0 0 +ydff/4 -ydff/4];
            axis(axis2);
        end;
        catch
            return;
        end;
    end;
    function push5(src,event,edt) %new function was hit.
      if ~edt
        [x y]=getfunction();
      else
        [x y]=getfunction(0,x,y);
      end;
      hide2=0;
      init(1,1,0);
    end;
    function push(src,event) %poly button was hit.

        if poly || piece
            init(0,0,hide2);
            DrawPoints();
        end;
        a=points(:,1);
        A=[a.^0];
        for i=2:index
            A(:,i)=a.^(i-1);
        end;
        b=points(:,2);
        u=A\b;
        for i=1:n
            z(i)=0;
            for j=1:index
                z(i)=z(i)+u(j)*x(i)^(j-1);
            end;
        end;
        plot(x,z,'r');
        max=0;
        maxindex=1;
        for i=1:n
            m=abs(z(i)-y(i));
            if m>max
                max=m;
                maxindex=i;
            end;
        end;
        plot(x(maxindex),z(maxindex),'.k');
        plot(x(maxindex),y(maxindex),'.k');
        for i=1:n
           m=abs(z(i)-y(i));
           if abs(m-max)<tol
               plot(x(i),z(i),'.k');
               plot(x(i),y(i),'.k');
           end;
        end;
        set(t1,'String',strcat(strcat('N=',num2str(index)),strcat(' Error=',num2str(max))));
        poly=1;
        piece=0;
    end;
    
    function push6(src,event) %Piece wise was hit.
      if poly || piece
        init(0,0,hide2);
        DrawPoints();
      end;
      [t i]=sort(points(:,3));
      points=points(i,:);
      firstx=points(1,3);
      for i=1:firstx;
         z(i)=points(1,2);
      end;
      for i=2:index
          newy=points(i,2);
          oldy=points(i-1,2);
          dy=newy-oldy;
          oldx=points(i-1,3);
          dx=points(i,3)-oldx;
          for j=1:dx
              z(oldx+j)=oldy+j*dy/dx;
          end;
      end;
      lastx=points(index,3);
      for i=lastx+1:n
          z(i)=points(index,2);
      end;
      plot(x,z,'r');
      max=0;
        maxindex=1;
        for i=1:n
            m=abs(z(i)-y(i));
            if m>max
                max=m;
                maxindex=i;
            end;
        end;
        plot(x(maxindex),z(maxindex),'.k');
        plot(x(maxindex),y(maxindex),'.k');
        for i=1:n
           m=abs(z(i)-y(i));
           if abs(m-max)<tol
               plot(x(i),z(i),'.k');
               plot(x(i),y(i),'.k');
           end;
        end;
        set(t1,'String',strcat(strcat('N=',num2str(index)),strcat(' Error=',num2str(max))));
        poly=0;
        piece=1;
    end;
    
    function push2(src,event); %clear button was hit.
        init(1,1,0);
    end;
    
    function push3(src,event,key); %Cheby or even was hit.
       pts=str2num(get(t2,'String'));
       xdff=x(n)-x(1);
       for j=1:pts
           c=cos( (2*j*pi-pi)/(2*pts))+1;
           c=c*xdff/2+x(1);
           if key
               c=x(1)+(j-1)*xdff/(pts-1);
           end;
           p2=[0 0 0];
           i=1;
           while i~=0
              if x(i)>=c || i==n
                  p2=[x(i) y(i) i];
                  i=0;
              else
                  i=i+1;
              end;
           end;
           duplicate=0;
           for i=1:index
             if points(index,1)==p2(1)
                 duplicate=1;
             end;
           end;
           if ~duplicate
             index=index+1;
             points(index,:)=p2;
             plot(p2(1),p2(2),'.r');
           end;
           set(t1,'String',strcat(strcat('N=',num2str(index)),' Error=NA'));
       end;
           
    end;
        
function [x y]=getfunction(res,x,y);
%
%INPUT ==>(opt:res,opt:x,opt:y) optional resolution (default=0.001)
%OUTPUT ==>[x y] returns a discrete function
%
%Allows user to draw a discrete function with the mouse.
%
figure(2);
ButtonDown=0;
fix=5;
fixd=3;
fixi=3;
shift=0;
tol=.01;
input=0;
try
    res;
    x;
    y;
    input=1;
    startx=y(1);
catch
    res=.001;
end;
if ~input
  x=-1:res:1;
  y=0*x;
end;
n=max(size(x));
oldb=1;
oldy=0;
rank=0;
startx=0;
wait=1;
plot(x,y);
axis([-1.2 1.2 -1.2 1.2]);
hold on;
plot([-1 1 1 -1 -1],[1 1 -1 -1 1],'k');
set(gcf,'WindowButtonDownFcn',@down);
set(gcf,'WindowButtonUpFcn',@up);
set(gcf,'WindowButtonMotionFcn',@motion);
set(gcf,'KeyPressFcn',@press);
set(gcf,'CloseRequestFcn',@close);
t3=uicontrol('Style','Text','String','Draw function with mouse and button; Use Keys c:clear, d:differentiate, i:integrate, s:smooth, x:exit and save','Position',[15,5,500,15]);
if input; replot('Start');end;
while wait
    pause(1);
end;


    function motion(src,event);
        pt=get(gca,'CurrentPoint');
        p=pt(1,1:2);
        if ButtonDown
          %disp(p);
          b=convert(p(1),1);
          
          if b>=oldb
            plot(x(oldb:b+1),y(oldb:b+1),'w');
            d=b-oldb;
            for i=1:d
              y(oldb+i)=oldy+i*(p(2)-oldy)/d;
            end;
            %fprintf('Write %f-->%f=%f-->%f\n',oldb+1,b,oldy,p(2));
            plot(x(oldb:b+1),y(oldb:b+1),'b');
          else
            plot(x(b-1:oldb),y(b-1:oldb),'w');
            d=oldb-b;
            for i=0:d-1
              y(b+i)=p(2)+i*(oldy-p(2))/d;
            end;
            %fprintf('Write %f<--%f=%f<--%f\n',b,oldb-1,p(2),oldy);
            plot(x(b-1:oldb),y(b-1:oldb),'b');
          end;
          
          oldb=b;
          oldy=p(2);
        end;
    end;
    function press(src,event);
       c=event.Character;
       if c=='d'
         ButtonDown=0;
         differentiate();
       elseif c=='i'
         ButtonDown=0;
         shift=0;
         integrate();
       elseif c=='k'
           ButtonDown=0;
           shift=1;
           integrate();
       elseif c=='s'
         ButtonDown=0;
         smooth();
       elseif c=='c'
         ButtonDown=0;
         input=0;
         rank=0;
         for i=1:n
             y(i)=0;
         end;
         cla;
         plot(x,y);
         axis([-1.2 1.2 -1.2 1.2]);
         hold on;
         plot([-1 1 1 -1 -1],[1 1 -1 -1 1],'k');
       else
         close(0,0);
       end;
    end;
    function b=convert(pointx,noend);
    %Converts -1<x<1 into index for x array. 
    %only return 2<=b<=n-1 if noend==1
    %returns 1<=b<=n for noend==0
    %and 0<=b<=n+1 for noend=-1
        i=1;
        b=1;
        while i~=0;
            if pointx<=x(i)
                b=i;
                i=0;
            else
                i=i+1;
                if i>n
                    b=n;
                    i=0;
                end;
            end;
        end;
        if noend==1
          if b==1 b=2; end;
          if b==n b=n-1; end;
        elseif noend==-1
            if pointx>x(n) b=b+1; end;
            if pointx<x(1) b=0; end;
        end;
    end;
    function differentiate();
        if rank==0
            startx=y(fixd);
        end;
        rank=rank+1;
        cleanend(fixd);
        z=zeros(n,1);
        for i=2:n-1
            d1=x(i)-x(i-1);
            d2=x(i+1)-x(i);
            s1=(y(i)-y(i-1))/d1;
            s2=(y(i+1)-y(i))/d2;
            z(i)=0.5*(s1+s2);
            %z(i)=s2;
        end;
        d=x(2)-x(1);
        z(1)=(y(2)-y(1))/d;
        d=x(n)-x(n-1);
        z(n)=(y(n)-y(n-1))/d;
        y=z;
        replot('Dif');
    end;
    function integrate();
        rank=rank-1;
        cleanend(fixi);
        z=zeros(n,1);
        z(1)=0;
        if rank==0
            z(1)=startx;
        end;
        for i=2:n
            d=x(i)-x(i-1);
            z(i)=0.5*(y(i-1)+y(i))*d+z(i-1);
        end;
        if shift
          mx=max(z(10:n-10));
          mn=min(z(10:n-10));
          sh=(mx+mn)/2;
          z=z-sh;
        end;
        y=z;
        replot('Int');
    end;
    function smooth();
      spoints=[0 0];
      sindex=0;
      temp=y(1);
      count=1;
      for i=2:n
        if abs(y(i)-temp)<tol
            count=count+1;
        else
            c=i-ceil(count/2);
            sindex=sindex+1;
            spoints(sindex,:)=[c temp];
            count=0;
            temp=y(i);
         end;
      end;
      firstx=spoints(1,1);
      for i=1:firstx;
         y(i)=spoints(1,2);
      end;
      for i=2:sindex
          newy=spoints(i,2);
          oldy=spoints(i-1,2);
          dy=newy-oldy;
          oldx=spoints(i-1,1);
          dx=spoints(i,1)-oldx;
          for j=1:dx
              y(oldx+j)=oldy+j*dy/dx;
          end;
      end;
      lastx=spoints(sindex,1);
      for i=lastx+1:n
          y(i)=spoints(sindex,2);
      end;
      replot('smooth');
    end;
    
    function down(src,event);
        pt=get(gca,'CurrentPoint');
        p=pt(1,1:2);
        b=convert(p(1),0);
        ButtonDown=1;
        oldb=b;
        oldy=y(b);
        %disp('down');
    end;
    function up(src,event);
        ButtonDown=0;
        %disp('up')
    end;
    function close(src,event);
        ButtonDown=0;
        cleanend(fix);
        wait=0;
        closereq;
    end;
    function cleanend(fix2)
        for i=1:fix2-1
            y(i)=y(fix2);
            y(n+1-i)=y(n+1-fix2);
        end;
    end;
    function replot(s); %s is string to print with max and min
        [mx indmx]=max(y(10:n-10));
        [mn indmn]=min(y(10:n-10));
        %fprintf('%s Max=%f at %f\n',s,mx,indmx);
        %fprintf('%s Min=%f at %f\n\n',s,mn,indmn);
        dy=.1*(mx-mn);
        cla;
        plot(x,y);
        if rank==0 && ~input
            mx=1;mn=-1;dy=.2;
        end;
        axis([-1.2 1.2 mn-dy mx+dy]);
        hold on;
        plot([-1 1 1 -1 -1],[mx mx mn mn mx],'k');
    end;

end
    

end
