%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Copyright (c) 2004.  Takayuki Osogami.  All Rights Reserved.   %
%                                                                   %
% Permission to use, copy, modify, and distribute this source code  %
% and its documentation for any purpose, without fee, and without   %
% a written agreement, is hereby granted, provided that the above   %
% copyright notice, this paragraph and the following two paragraphs %
% appear in all copies, modifications, and distributions.           %
% Created by Takayuki Osogami, Department of Computer Science,      %
% Carnegie Mellon University.                                       %
%                                                                   %
% IN NO EVENT SHALL THE AUTHOR BE LIABLE TO ANY PARTY FOR ANY       %
% DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS   %
% SOURCE CODE AND ITS DOCUMENTATION, EVEN IF THE AUTHOR HAS BEEN    %
% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                        %
%                                                                   %
% THE AUTHOR SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING,      %
% BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY     %
% AND FITNESS FOR A PARTICULAR PURPOSE.  THE SOURCE CODE AND        %
% ACCOMPANYING DOCUMENTATION, IF ANY, PROVIDED HEREUNDER IS         %
% PROVIDED "AS IS."  THE AUTHOR HAS NO OBLIGATION TO PROVIDE        %
% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Usage: [H] = homoNeuts(A0,A1,A2,G,epsilon)
%
% Given a homogeneous QBD process the function returns
% the first three moments of the first passage time
% to the one level below.
%
% INPUT: A0,A1,A2: The repeating part of the generator matrix,
%        G: G matrix of the repeating part
%        epsilon: accuracy (1e-6 recommended)
%
% author: Takayuki Osogami
%         Department of Computer Science
%         Carnegie Mellon University
%         osogami@cs.cmu.edu
% date: July 19, 2004

function [H] = homoNeuts(A0,A1,A2,G,epsilon)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Getting n-th moments (B0,B1,B2),                          %%
%% where [Bm^n]_{ij} = E[time to leave i | transition i->j]. %%
%% When n = 0, [Bm^0] = Pr(i->j | transition i->*).          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[F,P0,P1,P2] = FPmatrix(A0,A1,A2);

S = size(P1);
dim = S(1);

for n = 1:3

 for i = 1:dim
  rate = F(i,i);
  moment(i,i) = factorial(n) / (rate^n);
 end

 B0{n} = moment * P0;
 B1{n} = moment * P1;
 B2{n} = moment * P2;

end


%%%%%%%%%%%%%%%%%%
%% first moment %%
%%%%%%%%%%%%%%%%%%

C1 = B2{1} + B1{1}*G + B0{1}*G*G;
H{1} = Mmatrix(P1,P0,G,C1,epsilon);


%%%%%%%%%%%%%%%%%%%
%% second moment %%
%%%%%%%%%%%%%%%%%%%

M11 = MRImatrix(H,1,1,G);
M12 = MRImatrix(H,1,2,G);
M13 = MRImatrix(H,1,3,G);

C2 = B2{2} + B1{2}*G + B0{2}*G*G + 2*(B1{1}*M11 + B0{1}*M12 + P0*M11*M11);
H{2} = Mmatrix(P1,P0,G,C2,epsilon);


%%%%%%%%%%%%%%%%%%
%% third moment %%
%%%%%%%%%%%%%%%%%%

M21 = MRImatrix(H,2,1,G);
M22 = MRImatrix(H,2,2,G);
M23 = MRImatrix(H,2,3,G);

C3 = B2{3} + B1{3}*G + B0{3}*G*G + 3*(B1{2}*M11 + B0{2}*M12 + B1{1}*M21 + B0{1}*M22 + P0*M11*M21 + P0*M21*M11);
H{3} = Mmatrix(P1,P0,G,C3,epsilon);