function theta=UV2angle(U, V);
% theta=UV2angle(U, V);
% calculates theta for each pair of U and V vectors
% theta is in radians clockwise from V, or positive Y axis
%
NonzeroV=find(V); zeroV=find(V==0);

% the thetae is measured in radians clockwise from V
theta = zeros(size(U));

if length(zeroV)
    theta(zeroV)= (U(zeroV) ./ abs(U(zeroV))) * pi / 2; 
end
theta(NonzeroV) = atan(U(NonzeroV)./V(NonzeroV));

for ang=1:length(theta)
    if V(ang)<0 % in lower hemi-circle
        theta(ang)=theta(ang)+pi;
    else
        if U(ang) <0 % in upper left quadrant, want positive theta
            theta(ang)=theta(ang)+ 2*pi;
        end
    end
end

%theta=theta * 180/pi;