function M = moment(P,n)

% MOMENT - compute moments of a matrix polynomial
%
%        M = moment(P,n)
%
% The nth moment of P is defined as
%              _
%             \    n
%        M  = /_  k  P . 
%         n    k      k      
%
% n must be an integer >= 0 (default is n=0).

% Copyright (c) 2004 by Fritz Keinert (keinert@iastate.edu),
% Dept. of Mathematics, Iowa State University, Ames, IA 50011.
% This software may be freely used and distributed for non-commercial
% purposes, provided this copyright statement is preserved, and
% appropriate credit for its use is given.
%
% Last update: Feb 20, 2004

if (nargin < 2)
    n = 0;
end

M = zeros(size(P));
for k = 1:size(P.coef,3)
    true_k = k + P.min - 1;
    M = M + true_k^n * P.coef(:,:,k);
end