%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    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: [tau,T] = matching3EC(mean,m2,m3)
% 
% Given the mean and the normalized moments, (mean,m2,m3),
% the function returns the parameters (tau,T)
% of the EC that matches the normalized moments, where
%
%                                 p1
% --lambda--...--lambda--lambda1--lambda2-->
%                                 \ 
%                                  ---------->
%
% tau = (p,0,...,0)
%
%
%     [-lambda lambda   0               ..                   0     ]
%     [   0    -lambda lambda                                :     ]
% T = [   :                    ..                            :     ]
%     [   :                      -lambda lambda   0          0     ]
%     [   :             ..          0    -lambda lambda      0     ]
%     [   :                               0    -lambda1 p1*lambda1 ]
%     [   0             ..                        0      -lambda2  ]
%
%
% author: Takayuki Osogami
%         Department of Computer Science
%         Carnegie Mellon University
%         osogami@cs.cmu.edu
% date: June 6, 2004

function [tau,T] = matching3EC(mean,m2,m3)

if( isPH2(m2,m3) )

  [tau,T] = matching3PH2(mean,m2,m3);
  return;

end

%% determined the mass probability at zero.

EPSILON = 0.001;

if( (m3 > 2*m2-1) && (1/(m2-1)-floor(1/(m2-1)) < EPSILON^2) )

  p = (m2*m2+2*m2-1) / (2*m2*m2);

elseif( m3 < 2*m2-1 )

  p = 1 / (2*m2-m3);

else

  p = 1;

end

mean = mean / p;
m2 = m2 * p;
m3 = m3 * p;


%% determine the number of phases, n

 if( (m3 == 2*m2 - 1) && (m2 <= 2) )

  n = floor(m2/(m2-1) + EPSILON^2);

 else

  n = floor(m2/(m2-1) + 1 - EPSILON^2);

 end


%% calculate the mean and the normalized moments, (meanX,m2X,m3X), of 2-phase PH Coxian, X.

 m2X = ((n-3)*m2-(n-2)) / ((n-2)*m2-(n-1));
 meanX = mean / ((n-2)*m2X-(n-3));
 alpha = (n-2)*(m2X-1)*(n*(n-1)*m2X^2-n*(2*n-5)*m2X+(n-1)*(n-3));
 beta = ((n-1)*m2X-(n-2)) * ((n-2)*m2X-(n-3))^2;
 m3X = (beta*m3-alpha) / m2X;

%% calculate the rate of Erlang

 lambda = 1 / ((m2X-1)*meanX);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% When X is exponential, EC is just an Erlang %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if( (m2X == 2) && (m3X == 3) )

  mu = 1 / meanX;

  T = zeros(n-1,n-1);
  for i = 1:n-2
   T(i,i) = -mu;
   T(i,i+1) = mu;
  end
  T(n-1,n-1) = -mu;

  tau = zeros(1,n-1);
  tau(1) = p;

  return;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% When X is not exponential (2-phase PH) %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[tauX,TX] = matching3PH2(meanX,m2X,m3X);

T = [zeros(n-2,n-2) zeros(n-2,2); zeros(2,n-2) TX];
for i = 1:n-2
  T(i,i) = -lambda;
  T(i,i+1) = lambda;
end

tau = zeros(1,n);
tau(1) = p;