function drawsmoothcurveAB4(xdata,ydata) %4/28/2022, 4/29/2022  5/10/2022 
%                                      %5/21/2022  6/10/2022
%Three types of connections to a set of data points are graphed.
%1. Straight line connections.
%2. Bezier cubic polynomials between successive points.
%3. Natural spline curve which cubic polynomials between successive points.
%
%
%===>INPUT: xdata and ydata should be rows <===
%
%WARNING:
%The xdata should be ascending values from -15 to 15 because the grid
%is -15 by 15 horizontally and vertically. Any values outside the restraint
%will not appear the plot.
%
%Contibutors to the code:
%Modification of a routine by Tim Sauer.
%
%Modification of code that appeared in the URL
%https://stackoverflow.com/questions/7054272/how-to-draw-smooth-curve-through-n-points-using-javascript-html5-canvas 
%By ==>  https://stackoverflow.com/users/793454/homan
%See:  Homan Chou on Linkedin
%
%In addition contributions from
%DAVID R. HILL, MATH DEPARTMENT, TEMPLE UNIVERSITY, PHILADELPHIA, Pa. 19122
%dhill001@temple.edu
% 
%========================================================================
%

%determine if you want to sort the data for ascending x-coordinates
disp(' ')
disp('Do you want to sort the data so x-coordinates are ascending?')
reply = input('Enter 1 for YES; any other input will be considered NO. ==>  ','s');
       if isempty(reply)
          reply = 'N';
       end
       if reply == '1';
         [xdata,ydata]=sortdata(xdata,ydata)
       end 
       %checking max or min values WARNING:
if max(xdata)>15 | min(xdata)<-15 | max(ydata)>15 | min(ydata)<-15
        disp('Warning: a portion of graphs may not appear on the grid.')
        disp('Some data is not between -15 and 15.')
        disp('PRESS Enter to continue.')
        pause
end
%==========================================       

%checking that the xdata is ascending
% zx=xdata;zy=ydata;
% W=zy;
% [B I]=sort(zx);
% WW=[];
% for jj=1:length(zx)
% WW(jj)=W(I(jj));
% end
% xdata=B;ydata=WW;
% %checking max or min values WARNING:
% if max(xdata)>15 | min(xdata)<-15 | max(ydata)>15 | min(ydata)<-15
%     disp('Warning: a portion of graphs may not appear on the grid.')
%     disp('Some data is not between -15 and 15.')
%     disp('PRESS Enter to continue.')
%     pause
% end
%=====================================    
    
%Setting up the grid to use for graphing
gridutil([-1 1 -1 1]*15,1,1)%function by David Hill; see below
%NOTE: grid size can be changed 
hold on, grid on, axis('square')
addaxes,hold on
title('Press ENTER to graph the connections to the data points.','fontsize',20,'color','red')
pause

%pause  
 %Ipts is a column
 Ipts=[xdata;ydata]'; %making a column

 %=======================================================================
%Modification of code shown at the URL below 
%https://stackoverflow.com/questions/7054272/how-to-draw-smooth-curve-through-n-points-using-javascript-html5-canvas 
%Begin slopes computations for connection;a type of control point

sectionnumber=length(Ipts)-1;
%one less than the number of interpolation points
%NUMBER OF SECTIONS IN THE BEZIER CURVE
tt=1; f=.3; %parameters for slope computations??
m = 0; %initializing for slope
dx1 = 0;
dy1 = 0;
A=[]; %for storing info in the form 
   %[1st interp pt, 1st control type pt,2nd control type pt, 2nd inter pt]
   %the following gets this data
for jj=1:sectionnumber
    preP = [Ipts(jj, 1) Ipts(jj,2)];		
    		    curP = [Ipts(jj+1, 1) Ipts(jj+1,2)];
                    if (jj ==sectionnumber)
                     	dx2 = 0;
    					dy2 = 0;   
                    else   
    				nexP =[Ipts(jj+2, 1) Ipts(jj+2,2)];
    				    m = (preP(2)-nexP(2))/(preP(1)-nexP(1));
    					dx2 = (nexP(1) - curP(1)) * -f;
                        dy2 = dx2 * m * tt;
                    end 
    				
    			    C1=[preP(1)-dx1 preP(2)-dy1];
                    C2=[curP(1)+dx2 curP(2)+dy2];
                    
                    A=[A;[preP C1 C2 curP]];  %SAVING BEZIER DAtA
                    %update
                    dx1 = dx2;
    				dy1 = dy2;
    				preP = curP;
 end
Xmatrix=[A(:,1) A(:,3) A(:,5) A(:,7)]; %forming x and y coordinates
Ymatrix=[A(:,2) A(:,4) A(:,6) A(:,8)]; %for use in plotting

%PLOT THE INTERPOLATION POINTS
IPTX=Xmatrix(:,1); IPTY=Ymatrix(:,1);

for ii=1:sectionnumber   %<== markersize can be changed
    plot(IPTX(ii),IPTY(ii),'o','markersize',6,'MarkerFaceColor','k')
end
%PLOT LAST INTERPolAtion POINT
plot(Xmatrix(sectionnumber,4),Ymatrix(sectionnumber,4),'o','markersize',6,'MarkerFaceColor','k')
xmess=xlabel('PLOTTING DATA POINTS; ==> PRESS ENTER TO CONTINUE.','fontsize',20);
pause
delete(xmess)

%=====================================
%connect via lines
disp(' ')
disp('Do you want to connect points using lines?')
reply = input('Enter 1 for YES; any other input will be considered NO. ==>  ','s');
       if isempty(reply)
          reply = 'N';
       end
       if reply == '1';
        for jj=1:sectionnumber
            lin1=plot([Ipts(jj,1) Ipts(jj+1,1)],[Ipts(jj,2) Ipts(jj+1,2)],'k:','linewidth',2); 
        end 
        shg  %shg brings the current figure window forward
        xmess=xlabel('CONNECTING with LINES; ==> PRESS ENTER TO CONTINUE.','fontsize',20);
        pause
        delete(xmess)
        %disp('==> PRESS ENTER TO CONTINUE. <===')
        %pause
       end
       if reply~= '1';
       set(gcf,'visible','on');
       disp('==> PRESS ENTER TO CONTINUE. <===')
       pause
       end
%========================================    

%connect via bezier polgons
disp(' ')
disp('Do you want to connect points using cubic Bezier polygons?')
reply = input('Enter 1 for YES; any other input will be considered NO. ==>  ','s');
       if isempty(reply)
          reply = 'N';
       end
       if reply == '1';  %PLOT BEZIER "STYLE" CURVE SECTIONS
       for jj=1:sectionnumber
    x=Xmatrix(jj,:);y=Ymatrix(jj,:); 
    drawsection(x,y)  %calling a function to draw cubic sections
    shg %shg brings the current figure window forward
       end
       xmess=xlabel('CONNECTING with BEZIER; ==> PRESS ENTER TO CONTINUE.','fontsize',20);
       pause
       %going to plot the control points in red
       %
       %need code for unpacking control points
       [rr,col]=size(A);
       %first control point
       cp1x=A(1,5);cp1y=A(1,6);
       %final control point
       cpfx=A(rr,3);cpfy=A(rr,4);
       %use rows between 1 and last
       cpxmid=[];cpymid=[];
       for kk=2:rr-1
           cpxmid=[cpxmid,A(kk,3)];
           cpymid=[cpymid,A(kk,4)];
           cpxmid=[cpxmid,A(kk,5)];
           cpymid=[cpymid,A(kk,6)];
       end
       cpointsx=[cp1x, cpxmid, cpfx];
       cpointsy=[cp1y,cpymid,cpfy];
       
       numcontrolpts=length(cpointsx);%counting control points
       for ii=1:numcontrolpts   %<== markersize changed from 10 to 6
           plot(cpointsx(ii),cpointsy(ii),'o','markersize',6,'MarkerFaceColor','r')
        end
       xmess=xlabel('PLOTTING CONTROL POINTS==> PRESS ENTER TO CONTINUE.','fontsize',20);
       pause
       delete(xmess)
       xmess=xlabel('SHOWING SLOPES at CONTROL POINTS==> PRESS ENTER TO CONTINUE.','fontsize',20);
       %drawing control lines
       %use matrix A
       [rr,col]=size(A);
       %first line
       L1=plot([A(1,5) A(1,7)],[A(1,6) A(1,8)],'r:','linewidth',2);
       %final line
       Lfin=plot([A(rr,1) A(rr,3)],[A(rr,2) A(rr,4)],'r:','linewidth',2);
       %mid lines
       for kk=2:rr-1
          Lmid= plot([A(kk,1) A(kk,3)],[A(kk,2) A(kk,4)],'r:','linewidth',2);
          Lmid1= plot([A(kk,5) A(kk,7)],[A(kk,6) A(kk,8)],'r:','linewidth',2);
       end
       end %end of BEZIER plots
       disp('==> PRESS ENTER TO CONTINUE. <===')
       pause
       
%========================================

% %connect via natural spline polygons
disp(' ')
disp('Do you want to connect points using cubic Natural Spline polygons?')
reply = input('Enter 1 for YES; any other input will be considered NO. ==>  ','s');
       if isempty(reply)
          reply = 'N';
       end
       if reply == '1';  %PLOT BEZIER "STYLE" CURVE SECTIONS
           [DATA coeff] = nsplinefuncshort(xdata,ydata)
           xmess=xlabel('CONNECTING with SPLINES; ==> PRESS ENTER TO CONTINUE.','fontsize',20);
           shg %shg brings the current figure window forward
           pause
       delete(xmess)
       end
       if reply ~= '1';
           disp('==> PRESS ENTER TO CONTINUE. <===')
       
       set(gcf,'visible','on');
       pause
       end
%========================================

disp('The data set for points used.') 
[xdata;ydata]
xlabel('Lines in black, Bezier polygons in blue, Spline polygons in green dots.','fontsize',20,'color','magenta')
A
end

%end of main
%========================================================================

function [xdata,ydata]=sortdata(xdata,ydata)
%checking that the xdata is ascending
zx=xdata;zy=ydata;
W=zy;
[B I]=sort(zx);
WW=[];
for jj=1:length(zx)
WW(jj)=W(I(jj));
end
xdata=B;ydata=WW;
%checking max or min values WARNING:
if max(xdata)>15 | min(xdata)<-15 | max(ydata)>15 | min(ydata)<-15
    disp('Warning: a portion of graphs may not appear on the grid.')
    disp('Some data is not between -15 and 15.')
    disp('PRESS Enter to continue.')
    pause
end
end %function
%=====================================    



%========================================================================
function drawsection(x,y)
%DRAW the parameterized cubic; code from Tim Sauer
t=0:.01:1;
    bx=3*(x(2)-x(1));
    by=3*(y(2)-y(1));
    cx=3*(x(3)-x(2))-bx;
    cy=3*(y(3)-y(2))-by;
    dx=x(4)-x(1)-bx-cx;
    dy=y(4)-y(1)-by-cy;
    xpoly=x(1)+t.*(bx+t.*(cx+t*dx));
    ypoly=y(1)+t.*(by+t.*(cy+t*dy));
    seg=plot(xpoly,ypoly,'b-','linewidth',1); %graph the segment 
    
end

%========================================================================
function gridutil(ax,xinc,yinc)  %6/16/98, 2/13/2005, 5/29/07, 6/20/2015
%GRIDUTIL   A utility to create a two dimensional plotting grid.
%           
%           
%           Use in the form  ==>  gridutil(ax,xinc,yinc)  <==
%
%               where    
%                      ax is a row with four entries
%                      ax = [xmin  xmax  ymin  ymax]
%
%                      xinc is the increment for tick marks on the x-axis
%
%                      yinc is the increment for tick marks on the x-axis
%  
%  The color and intensity of the grid lines can be changed; see last
%  command line. (gridalpha value is in [0, 1])
%
%  By: David R. Hill, Math Dept, Temple University,
%      Philadelphia, Pa. 19122   Email: dhill001@temple.edu

if nargin<3
   disp('ERROR in gridutil: not enough input arguments.')
   return
end  
%hfig=figure('units','normal','position',[0 0 1 1]);
hfig=figure('units','normal','position',[0 0 1 .9],'color','white');
axis('off')
grbox=subplot('position',[.15 .15 .7 .7]);
axis(ax);
set(grbox,'xtick',[ax(1):xinc:ax(2)],'ytick',[ax(3):yinc:ax(4)],...
    'xgrid','on','ygrid','on');
%set(grbox,'gridlinestyle','-') %sets the line style to solid
set(grbox,'gridlinestyle','-','gridColor','black','gridalpha',1)
%added 6/20/2015 from Technical Support Case #01405773 at the MathWorks
end

function dolines 
disp(' ')
disp('Do you want to connect points using lines?')
%disp('Enter 1 for YES; any other input will be considered NO.')
reply = input('Enter 1 for YES; any other input will be considered NO. ==>  ','s');
       if isempty(reply)
          reply = 'N';
       end
       if reply == '1';
        for jj=1:sectionnumber
            lin1=plot([Ipts(jj,1) Ipts(jj+1,1)],[Ipts(jj,2) Ipts(jj+1,2)],'k:','linewidth',2); 
        end   
       end
end      

function stlines
%Connect interpolation points with straight lines
 for jj=1:sectionnumber
   lin1=plot([Ipts(jj,1) Ipts(jj+1,1)],[Ipts(jj,2) Ipts(jj+1,2)],'k:','linewidth',2); 
 end 
end

function [DATA coeff] = nsplinefuncshort(x,y) %<< 11/8/98 & 11/22/07 & 8/14/09 & 11/18/09
                               %<< 4/19/2022   5/9/2022
%NSPLINEFUNCSHORT.M  
%A "cut down of a standalone" natrual spline program nspline.m by 
%the aouthor. This function is a part of another program for connecting
%points of a data set of points on a grid.
%
%         The natural cubic spline expects data in vectors
%         x and y.  The x-values must be in increasing
%         order. The output is a table of coefficients for the
%         piecewise cubics. A list of the points on the spline
%         is generated. Graphing takes place within a program to
%         connect data points. See function drawsmoothcurve(xdata,ydata).
%
%   Use in the form   ===>  [DATA coeff]=nsplinefuncshort(x,y)  <===
%
%OUTPUT: DATA, a list of the coordinates of the points on the spline curve
%        and coeff a table of the coordinates of the cubic polynomial
%        pieces that make up the natural spline.
%   
%   By: David R. Hill, Math Dept, Temple University,
%       Philadelphia, Pa. 19122   Email: dhill001@temple.edu

n=length(x);m=length(y);
if n~=m
   disp('Error: input vectors of different lengths.')
end
if any(sort(x)~=x)
    disp('Error: independent variable must be in increasing order.')
    return
end
[xm,xn]=size(x);
if xm~=n | xn~=1,x=x';
end
%make sure we have columns
[ym,yn]=size(y);
if ym~=m | yn~=1,y=y';
end
%Calculate the Natural Spline
   nm1=n-1;
   for j=1:nm1,h(j)=x(j+1)-x(j);end
   d=[1];d1=[0];d2=[];
   for j=1:n-2,d=[d 2*(h(j)+h(j+1))];d1=[d1 h(j+1)];d2=[d2 h(j)];end
   d=[d 1];d2=[d2 0];
   A=diag(d)+diag(d1,1)+diag(d2,-1);b=zeros(n,1);
   for j=2:n-1
      b(j)=(3/h(j))*(y(j+1)-y(j)) -(3/h(j-1))*(y(j)-y(j-1));
   end
   C=A\b;   %==== KEY STEP
   for j=nm1:-1:1
      B(j)=(y(j+1)-y(j))/h(j) - h(j)*(C(j+1)+2*C(j))/3;
      D(j)=(C(j+1)-C(j))/(3*h(j)); 
   end
   coeff=[D(1:nm1)' C(1:nm1) B(1:nm1)' y(1:nm1)];
  disp('The coefficients of the cubic pieces ')
  disp('D*(x-x(i))^3 + C*(x-x(i))^2 + B*(x-x(i)) + A are:')
  z='           i           x           D         ';
  z=[z '  C           B           A'];
  disp(z)  %display string containing column labels
  format short g
  z=1:nm1;disp([z' x(1:nm1) coeff])
  disp('spline pieces')
  %pause
  %=================
   a=x(1);b=x(n);  %interval over which the spline is evaluated
   t=a:(b-a)/100:b;t=[t x'];t=sort(t);
   %set 100 points in addition to the original data; old value
  for k=1:length(t)
     for ks=1:nm1       
        if t(k)>=x(ks) & t(k)<x(ks+1)       
           ft=polyval(coeff(ks,:),t(k)-x(ks));       
        end       
     end       
     if t(k)<x(1)       
        ft=y(1)+B(1)*(t(k)-x(1));       
     end       
     if t(k)>=x(n)       
        ft=(B(n-1)+2*C(n-1)*h(n-1)+3*D(n-1)*h(n-1)^2);
        ft=ft*(t(k)-x(n))+y(n); 
     end       
     xp(k)=t(k);yp(k)=ft;   %<=== ordered pairs for spline figure
  end  
  %ADDING the DATA accumulation
  DATA=[xp;yp];       %<==============================
plot(x,y,'o','markersize',10,'MarkerFaceColor','k') %plot data points
%plot(xp,yp,'or','linewidth',1) %plot spline curve
plot(xp,yp,'o','markersize',4,'MarkerFaceColor','g') %plot spline curve
%=====================            

end

function addaxes
%addaxes   Draw X and Y axes on the current plot.  2/9/99 
%last updated 8/16/05  last updated  3/2/2022
%
%      Add x and y axes to a 2D plot
%
%      Use in the form    ==>  addaxes  <==
%
%  By: David R. Hill, Mathematics Dept., Temple Univ.
%      Philadelphia, PA. 19122    Email: dhill001@temple.edu

figure(gcf)
hold on,gca;
plot(xlim,[0 0],'-k',[0 0],ylim,'-k','linewidth',2);
hold off
end
