function v = cordic(alpha,n)
% returns v = [cos(alpha), sin(alpha)] 
% alpha=angle (radians)
% n=number of iterations, a non null positive integer
% Author: Michel Barbeau
% Version: December 23, 2014

% initial unit vector
% in IInd quadrant or IIIrd quadrant?
if (alpha>=pi/2)&&(alpha<=3*pi/2)
	% start with a pi radian rotation
    v = [-1;0];
	alpha=alpha-pi;
else % in Ist or IVrd quadrant
    v = [1;0];
    % in IVth quadrant?
    if alpha>=3*pi/2
        % translate to a negative angle
        alpha=alpha-2*pi;
    end
end

% array of pre-computed arctangents (radians) of negative powers of two
gamma = atan( pow2(0:-1:-(n-1)) );
% Factor K
K=prod( 1./( sqrt( 1 + pow2(-2.*[0:1:n-1]) ) ) );
poweroftwo = 1;
% steer the unit vector to the target angle "alpha"
for i= 1:1:n
    % rotation direction
    sigma = sign(alpha);
    % build rotation matrix
    R = [1, -sigma*poweroftwo; sigma*poweroftwo, 1];
    % apply the rotation
    v = R * v; 
    % for the next iteration,
    % - update angle of unit vetor
    alpha = alpha - sigma * gamma(i); 
    %  - calculate 2^-i
    poweroftwo = poweroftwo / 2;
end
% apply factor K
v = K*v;
return