  function P3Grade()
% Checks out SpecialVec

% Random mu example...
n = 10;
A = randn(n,n);
[U,S,V] = svd(A);
mu = S(n,n) + rand(1)*(S(1,1)-S(n,n));
x = SpecialVec(A,mu);
xCheck = abs(norm(x)-1);
muCheck = abs(norm(A*x)-mu);
spanCheck = norm(V(:,2:n-1)'*x);
fprintf('                   xCheck     muCheck     spanCheck\n')
fprintf('---------------------------------------------------------\n')
fprintf('mu = Random:    %10.2e  %10.2e  %10.2e\n',xCheck,muCheck,spanCheck)

% Smallest possible mu...
mu = S(n,n);
x = SpecialVec(A,mu);
xCheck = abs(norm(x)-1);
muCheck = abs(norm(A*x)-mu);
spanCheck = norm(V(:,2:n-1)'*x);
fprintf('mu = S(n,n):    %10.2e  %10.2e  %10.2e\n',xCheck,muCheck,spanCheck)

% Largest possible mu
mu = S(1,1);
x = SpecialVec(A,mu);
xCheck = abs(norm(x)-1);
muCheck = abs(norm(A*x)-mu);
spanCheck = norm(V(:,2:n-1)'*x);
fprintf('mu = S(1,1):    %10.2e  %10.2e  %10.2e\n',xCheck,muCheck,spanCheck)

  function x = SpecialVecCVL(A,mu)
% A is nxn with SVD U*S*V'.
% mu is a real that satisfies S(n,n) <= mu <= S(1,1)
% x is a column n-vector with the property that  (a) norm(x) = 1, (b) x is a 
% linear combination of V(:,1) and V(:,n), and (c) norm(A*x) = mu.
[n,n] = size(A);
[U,S,V] = svd(A);
sig_1 = S(1,1);
sig_n = S(n,n);
if S(n,n)==S(1,1)
    % In this case any cV(:,1)+sV(:n) works with c^2+s^2=1. Use (c,s)=(1,0)
    x = V(:,1);
    return
end
% if y = A*x and x = c*V(:,1)+s*V(:,n),  then x is a unit vector
% and norm(y)^2 = c^2*sig_1^2 + s^2*sig_n^2. We want
%
%           c^2*sig_1^2 + s^2*sig_n^2 = mu^2
% i.e.,
%           c^2*(sig_1^2 - mu^2) = s^2*(mu^2 - sig_n^2)
%
% so tau^2 = (c/s)^2 = (mu^2 - sig_n^2)/(sig_1^2 - mu^2)
% and
%           s = 1/sqrt(1+tau^2) = sqrt(sig_1^2-mu^2)/sqrt(sig_1^2-sig_n^2)
s = sqrt(sig_1^2-mu^2)/sqrt(sig_1^2-sig_n^2);
c = sqrt(1-s^2);
x = c*V(:,1) + s*V(:,n);

% Points off for using things like fsolve to get c and s
% Roughly 3 points for equations that correctly define c and s and 2 points for 
% solving.
   
   
   