function varargout=decompose_orthogonal_matrix(M)

% use [R,S,T] = decompose_orthogonal_matrix(M) to grab the results for M 3x3
% then M=R*S*T
% similar for nxn matrices

[m,n]=size(M);

N=M;
I=eye(n);

a=I(:,1);
b=N(:,1);
if b==a
  varargout{1}=eye(n);
else
  varargout{1}=reflect_median(a,b);
endif

for i=2:n
N=varargout{i-1}*N;
a=I(:,i);
b=N(:,i);
if b==a
  varargout{i}=eye(n);
else
  varargout{i}=reflect_median(a,b);
endif
endfor

endfunction