function P4Grade()
% Checks out BestPlane
m = 100;
% Generate and example...
qx = randn(m,1); x0 = sum(qx)/m;
qy = randn(m,1); y0 = sum(qy)/m;
v = 3*(qx-x0) + 2*(qy-y0);
qz = -(v - sum(v)/m)+ .001*randn(m,1);
Q = [qx qy qz];
% Solve two ways..
[a_YOU,Q0_YOU,phiMin_YOU] = BestPlane(Q);
a_YOU = a_YOU/norm(a_YOU);
[a_CVL,Q0_CVL,phiMin_CVL] = BestPlaneCVL(Q);
clc
fprintf('a_YOU  = [%10.6f, %10.6f, %10.6f]\n',a_YOU)
fprintf('a_CVL  = [%10.6f, %10.6f, %10.6f]\n\n',a_CVL)
fprintf('Q0_YOU = [%10.6f, %10.6f, %10.6f]\n',Q0_YOU)
fprintf('Q0_CVL = [%10.6f, %10.6f, %10.6f]\n\n',Q0_CVL)
fprintf('phiMin_YOU = %10.6f\n',phiMin_YOU)
fprintf('phiMin_CVL = %10.6f\n\n',phiMin_CVL)

function [a,Q0,phiMin] = BestPlaneCVL(Q)
% Q is an m-by-3 matrix. 
% determines a row-3 vector a with the property that
%   phi(a) = norm( (Q-ones(m,1)*Q0)*a' )/norm(a) 
% is minimized where Q0 = sum(Q). phiMin is the minimum value of phi 
[m,n] = size(Q);
Q0 = sum(Q)/m;
[~,S,V] = svd(Q - ones(m,1)*Q0);
a = V(:,3);
phiMin = S(3,3);

% -1 for [U,S,V] = svd(Matrix) because that is O(m^2n) we do not need U
% OK to use thin SVD:  svd(matrix,0)