  function Q_new = Extend(Q,v)
% Extend an orthonormal basis so that the new expanded subspace includes
%   a given a new vector.
% Q is nxk with orthonormal columns and v is not in the range of Q.
% Q_new = [Q q_kp1] with the property that v is in the rand of Q_new

% Since v is in the range of Q_new we must have  v = Q*alfa + tau*q_kp1
%    for some alfa (kx1) and scalar tau.
% Apply Q' to both sides says Q'v = alfa.
% This says q_kp1 is a unit vector in the direction of v - Q*(Q'v)

z = v - Q*(Q'*v);  % Note importance of parentheses!
q_kp1 = z/norm(z); % Normalize z 
Q_new = [Q q_kp1]; % Append