function Markov(file_header)
% This function ->
%       Markov: Given the packets file, the 2-state Markov model is
%       produced.  Outputting the transitional probabilities with
%       equilibrium probabilitites

filename = [file_header];

fid = fopen(filename,'r');
tline = fgets(fid);

count = 0;   %typical counter
room_back = 0;
mp_back = 0;
line_counter = 0;

%  The Bur vector populates the index of the vector when a Bur of that
%  index has been found

while(1)
   tline = fgets(fid);
   room = sscanf(tline,'%f');  % extract room number 
   mp = Extract(tline, 1);      % extract MP
   
   if((room == room_back) && (mp == mp_back))
       break;
   end
   
   room_back = room;
   mp_back = mp;
   line_counter = line_counter + 1;
end

%dont miss over any lines
frewind(fid)
i = 1;
while(i <= line_counter)
    tline = fgets(fid);
    i = i + 1;
end

i = 1;
default_state = '1';
one = 0; %default state       
zero = 0;
zeronum = 0;
one_one = 0;
one_zero = 0;
zero_one = 0;
zero_zero = 0;

while(~(feof(fid)))
    previous_state = default_state;
    tline = fgets(fid);
    new_tline = tline(9:length(tline));    
     
    % Need to loop for the length subtracted by 2, one for the \0 char and
    % one to eliminate the final comparison of the final state and nothing.
    for (i=1:(length(new_tline) - 2))
        previous_state;
        current_state = new_tline(i);
        if(current_state == '1')
            one = one + 1;
            zeronum = 0;
            if(current_state == previous_state)
              one_one = one_one + 1;
            else
              zero_one = zero_one + 1;
            end
        else
            zero = zero + 1;
            zeronum = zeronum + 1;
            if(current_state == previous_state)
              zero_zero = zero_zero + 1;
            else
              one_zero = one_zero + 1;
            end
        end    
        previous_state = current_state;
    end
end



Markov_ma = zeros(2,2);
Markov_ma(1,1) = zero_zero/zero;
Markov_ma(2,1) = zero_one/zero;
Markov_ma(1,2) = one_zero/one;
Markov_ma(2,2) = one_one/one;

disp('Markov matrix is')
Markov_ma

% find the equilibrium values of the states Ax = b where A = stat_ma
Markov_ma = Markov_ma - eye(2);
stat_ma = ones(3,2);
b = [1;0;0];
stat_ma(2,1) = Markov_ma(1,1);
stat_ma(3,1) = Markov_ma(2,1);
stat_ma(2,2) = Markov_ma(1,2);
stat_ma(3,2) = Markov_ma(2,2);
% stat_ma [zeros, ones]'
equilibrium_states_are = linsolve(stat_ma, b)
fclose(fid);