%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    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] = matching3PH2(mean,m2,m3)
%
% Given mean and normalized moments, (mean,m2,m3),
% the function returns a '2-phase' PH distribution (tau,T) that has
% the first three moments (mean and normalized moments).
%
% If there is no 2-phase PH distribution with the input moments,
% the function returns the "closest" 2-phase PH distribution.
%
% Note: If the exponential distribution (1-phase PH) matches the input,
% the function still returns 2-phase PH distribution.
%
% When the output PH is (0,0), there was an error in the function
% or the input was illegal.
%
% author: Takayuki Osogami
%         Department of Computer Science
%         Carnegie Mellon University
%         osogami@cs.cmu.edu
% date: June 6, 2004


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

%% If the input moment is zero, return error.

if( mean == 0 )

  Illegal_mean = mean
  tau = 0;
  T = 0;
  return;

end

%% If there is no 2-phase PH distribution with the input moments,
%% calculate the normalized moments of the "closest" 2-phase PH distribution.

EPSILON = 0.001;

if( m2 < 3/2 )

  m2 = 3/2;
  m3 = 2;

elseif( m2 >= 2 - EPSILON^2 )

  if( abs(m2-2) < EPSILON^2 )
    m2 = (1+EPSILON) * 2;
  end

  if( m3 <= 3/2 * m2 + EPSILON^2 )
    m3 = (1+EPSILON) * 3/2 * m2;
  end

else

  if( m3 > 6*(m2-1)/m2 )
   m3 = 6*(m2-1)/m2;
   m3 = 2*m2 - 1;
  elseif( m3 < (9*m2-12+3*(2-m2)*sqrt(2*(2-m2))) / m2 )
   m3 = (9*m2-12+3*(2-m2)*sqrt(2*(2-m2))) / m2;
   m3 = 2*m2 - 1;
  end

end

u = (6-2*m3) / (3*m2-2*m3);
v = (12-6*m2) / (m2*(3*m2-2*m3));

sqrtD = sqrt(u^2-4*v);

mu1 = (u+sqrtD)/2;
mu2 = (u-sqrtD)/2;
p = mu2*(mu1-1)/mu1;

MU1 = mu1 / mean;
MU2 = mu2 / mean;

tau = [1 0];
T = [-MU1 p*MU1;
     0    -MU2];