% Computer problem 3.5 on p.153 Heath
% A least squares problem.
% K. Ming Leung, 04/26/03, 04/20/05 (modified orbit plotting function)
% A skeleton program for the problem.

clear all; hold off; bigscreen;
% the original data
Xvec = [1.02 0.95 0.87 0.77 0.67 0.56 0.44 0.30 0.16 0.01]';
Yvec = [0.39 0.32 0.27 0.22 0.18 0.15 0.13 0.12 0.13 0.15]';
% low and upper limits along x and y in contour plot
lowLimx = -1; highLimx = 1.5;
lowLimy = -0.5; highLimy = 2;

% dimensions of the linear least squares problem
m = length(Xvec); n = 5;
% matrix A for the linear least squares problem
....
% parameter obtained from the original data
Poriginal = ...
% plot the orbit
conicSectionPlot_(Poriginal,lowLimx,highLimx,lowLimy,highLimy,'b-');
hold on;
plot(Xvec,Yvec,'ok'); % plot original data


% Introduce random noise (perturbaton)
% set the lower and upper limits of the random noice
low = -0.005; high = -low;
rand('state',17);
Xpvec = Xvec + low + rand(m,1)*(high-low); 
Ypvec = Yvec + low + rand(m,1)*(high-low);
plot(Xpvec,Ypvec,'xr');    % plot the randomized data
% set up matrix and vector of the linear least squares problem
Pperturbed = 
conicSectionPlot_(Pperturbed,lowLimx,highLimx,lowLimy,highLimy,'r-');

% SVD for the original problem
PlotStyle = ['b--'; 'r--'; 'g--'; 'c--'; 'm--'];
% find the SVD and compute the parameters Psvd 
% see p.138 Heath
for k=1:nn, 
    conicSectionPlot_(Psvd,lowLimx,highLimx,lowLimy,highLimy,PlotStyle(k,:)); 
end 
%
% SVD for the perturbed problem
PlotStyleP = ['b:'; 'r:'; 'g:'; 'c:'; 'k:'];
% repeat the above using the perturbed data

% total least squares treatment
% compute the parameters Pt
% plot the resulting orbit
conicSectionPlot_(Pt,lowLimx,highLimx,lowLimy,highLimy,'y-'); 
