%
disp('Create some data -- e.g. try to plot a sine curve.')
disp('Suggestion: use about ten data points -- use the') 
disp(' right mopuse button to end creating data-points.')
disp('Strike any key when ready.')
pause
%
% Create a plotting window with some reference marks and guide lines
%
clf
hold on
axis([0,13,-1.2,1.2])
plot([0,13],[0,0],'c-');
plot([0,13],[1 1],'c:');
plot([0,13],[-1 -1],'c:');
for j=1:4 plot(j*[3.1416 3.1416],.04*[-1 1],'c-'),end;
%
% Data creation and input
%
p=drawfig('r*',[0,13,-1.2,1.2])
%
disp('These are your data points.')
%
tt=p(1,:)';
yy=p(2,:)';
[tt,yy]
disp('Here we try to find the best for your data points, using')
disp('functions of the form:')
disp('f(t)=c1+c2*cos(t)+c3*cos(2t)+c4*sin(t)+c5*sin(2t).')
disp('Strike any key when ready.')
pause
%
% Compare to the graph of the correct sine, highlight the errors.
%
ttt=[0:.1:13];
plot(tt,sin(tt),'b*');
plot(ttt,sin(ttt),'b-');
for j=1:length(tt)
    plot([tt(j) tt(j)],[sin(tt(j)),yy(j)],'g-')
    end
err=sqrt(sum((sin(tt)-yy).^2))
text(2,1.1,['Your least squares error from the target-sine is ',num2str(err)])
pause
%
disp(' Next, evaluate the basis functions at the "t-values" of the')
disp(' data-points. Each column contains the corresponding "y-values"')
disp(' for one of the basis functions: 1, cos t, cos 2t, sin t, sin 2t')
disp('Strike any key when ready.')
pause
%
A=[ones(size(tt)),cos(tt),cos(2*tt),sin(tt),sin(2*tt)]
%
disp('We now find solve the overdetermined system A*c=y in the "least')
disp('squares sense, first carrying all steps out by hand.')
disp('First calculate the transpose A'' and the product A''A.') 
disp('Strike any key when ready.')
pause
AT=A'
ATA=A'*A
disp('Calculate A''y, the inverse of A''A.')
disp('Strike any key when ready.')
pause
ATy=A'*yy
ATAinv=inv(A'*A)
disp('Finally calculate the least squares soluntion cc=(A''A)^(-1)A''y')
disp('Strike any key when ready.')
pause
cc=inv(A'*A)*A'*yy
disp('It remains to calculate (plot) the trig polynomial that best fit')
disp('the data. Simply use the coefficients ci computed above and plot')
disp('the linear combination: c1+c2*cos(t)+c3*cos(2t)+c4*sin(t)+c5*sin(2t).')
disp('For curiosity we overlay the individual component functions.')
disp(' ')
disp('Strike any key when ready.')
pause

%
% Recreate most of the original window with data points
%
clf
hold on
axis([0,13,-1.2,1.2])
plot([0,13],[0,0],'g:');
plot([0,13],[1 1],'g:');
plot([0,13],[-1 -1],'g:');

plot(tt,yy,'r*')
y1=cc(1)*ones(size(ttt));
y2=cc(2)*cos(ttt);
y3=cc(3)*cos(2*ttt);
y4=cc(4)*sin(ttt);
y5=cc(5)*sin(2*ttt);
%
% For direct calculation of the new least squares error
yy1=cc(1)*ones(size(tt));
yy2=cc(2)*cos(tt);
yy3=cc(3)*cos(2*tt);
yy4=cc(4)*sin(tt);
yy5=cc(5)*sin(2*tt);
yyappr=yy1+yy2+yy3+yy4+yy5;
%
plot(ttt,y1,'b-');
plot(ttt,y2,'b-');
plot(ttt,y3,'b-');
plot(ttt,y4,'b-');
plot(ttt,y1+y2+y3+y4+y5,'r-');
pause
disp('Strike any key when ready.')
for j=1:length(tt)
    plot([tt(j) tt(j)],[yyappr(j),yy(j)],'g-')
    end
err2=sqrt(sum((yyappr-yy).^2))
text(2,1.1,['Your new least squares error using 5 functions is ',num2str(err2)])
%
%
disp('Now some really neat news: If a system A*c=y is overdetermined,')
disp('MATLAB will automatically calculate the least squares solution')
disp('upon entering A\y. BEWARE: In many cases in school you may need')
disp(' to know that an overdetermined system has NO solution -- don''t')
disp('be fooled and simply type A\y and wait for error messages ....')
disp(' ')
disp('Compare the previously found solution cc and the answer A\y by MATLAB')
disp(' ')
disp('Strike any key when ready.')
pause
%
cc
ccc=A\yy
cc-ccc
disp('Exercises: Compare the previous solution to one using more basis')
disp('           functions, e.g. add in cos 3x and sin 3x. Hint: You ')
disp('           merely need to enlarge the matrix A, and then reexecute')
disp('           the command A\y, as well as graphing the new approximation.')
disp(' ')
disp('           Explain why you expect that the curve exactly passes thru')
disp('           the data points if the number of basis functions equals')
disp('           the number of data points. Try it out!')
disp(' ')
disp('           As another exercise, find the best polynomial fit in a ')
disp('           least squares sense. Use polynomials of various degrees.')