moivre( n, z )

function zroots = moivre( n, z )

%MOIVRE: uses Moivre's formula to calculate the roots of a complex number. 
% 
% Usage:   zroots = moivre( n, z )
%
%          Example: the following statement calculates the cubic roots of 1+i:
%          >> moivre(3,1+i)
%          The function takes the integer part of n, so the statement 
%          >> moivre(3.9,1+i)
%          produces the same result as the first one.
%          If the user inputs a multimensional array, where a particular dimension 
%          is equal to 1 the function stores the n roots along that dimension. 
%          Otherwise a new dimension is added to the output array. 
%          For instance the statement 
%          >> moivre(6,[1+i 1-i -1]) 
%          produces a 6x3 matrix, while the output of 
%          >> moivre(6,[1+i 1-i -1].') 
%          produces a 3x6 matrix. However, the output of  
%          >> moivre(6,rand(10)+sqrt(-1)*rand(10))
%          produces a 10x10x6 array.  
%          

%***************************************************************************************
%***************************************************************************************
% First version: 28/07/2008
% 
% Contact: [e-mail]
% 
% Any suggestions to improve the performance of this 
% code will be greatly appreciated. 
%
% Reference: Milton Abramowitz and Irene A. Stegun, Handbook of Mathematical Functions, 
%            (1964) Dover Publications, New York. ISBN 0-486-61272-4. (p. 74).
% 
%***************************************************************************************

zroots = []; 

imunit = sqrt( -1 ); 

n = fix( n ); 

r = abs( z ); x = angle( z ); s = size( z );

[mins index] = min(s); 

if mins == 1 
s(index) = n; 
news = s; 
else 
news = [s n];
end 

if n <  0, disp( 'This functions does not accept negative roots...'               ), return, end 
if n == 0, disp( 'Forgive my ignorance, I have no idea what a zero-th root is...' ), return, end 

for k = 0:n-1
    
    argumento = ( x + 2*k*pi )/n;
    
    zroots = [zroots; ( r.^(1/n) ).*( cos(argumento)+imunit*sin(argumento) )]; 

end 

zroots = reshape(zroots,news);