# copyright Barry L. Nelson and Michael R. Taaffe 1999

with(linalg):
with(plots):

PhPhInfSojourn  := proc(A, lambda, B, mu, arrivalT)
#
# this procedure returns the mean, variance and a plot structure of
# the cdf of the sojourn time of a virtual arrival at time arrivalT
# to the Ph_t/Ph_t/Infinity queue as given in Nelson and 
# Taaffe (1999)
#
# last update 6/30/99
# 
# Inputs:
# A = Markov chain transition matrix for the arrival process
# B = Markov chain transition matrix for the service process
# lambda = arrival rate vector
# mu = service rate vector
# arrivalT = arrival time of virtual arrival
#
#
# Output:
# [1] mean of sojourn time
# [2] variance of sojourn time
# [3] plot structure for the cdf of sojourn time
#
           local  mA, mB, N, k, i, j, g1, a, c, d, b, x,
                  InitP, InitN, InitNA, InitNN, mdeq, cdf, first,
                  second, Tinf;
#
# create index mapping function
           g1 := proc(i,j)
              RETURN(N.i.N.j);
           end;
 #
# determine number of phases in phase representations
           mA := rowdim(A) -1;
           mB := rowdim(B) -1;
# extract beta
          b := transpose(submatrix(B, mB+1..mB+1, 1..mB));
# create the initial conditions
          InitP := vector(mA, 0);
          InitN := vector(mB);
          for i from 1 to mB do
             InitN[i] := b[i,1];
          od;
          InitNA := matrix(mA, mB, 0);
          InitNN := vector(mA);
          for i from 1 to mA do
             InitNN[i] := matrix(mB, mB, 0);
          od;
# create the numerical procedure
          mdeq := PhPhInfNumericInit(A, lambda, B, mu, arrivalT, InitP, InitN, InitNA, InitNN); 
#
# create the N vector for plotting
          N := matrix(mB, 1, [N.(1..mB)]);
          N := N(t);
#
# compute first and second moment
         cdf := subs(mdeq, 1 - sum(N[a,1], a=1..mB)):
         first := subs(mdeq, M1(t));
         second := subs(mdeq, M2(t));
         Tinf := arrivalT;
         while cdf(Tinf) < 0.99999 do
          Tinf := Tinf + 0.1
         od;
         RETURN(first(Tinf), second(Tinf)- first(Tinf)^2, 
         odeplot(mdeq, [t, 1 - sum(N[a,1], a=1..mB)], arrivalT..arrivalT + Tinf));
         end;