function result = bisection(fx, y, property, target, lowLimit, highLimit, tol, k_max)
% Description:      Uses the bisection method to solve f(x)=target with use
%                   with tpws.m.  For instance, used to solve for the
%                   density corresponding with a given pressure 
%                   and temperature in Region 3 calculations.
% -------------------------------------------------------------------------
% Syntax:           result = bisection(fx, y, property, target, lowLimit, highLimit, tol, k_max)
%                   fx          Function of x with optional parameter y
%                               Ex: f(x)=y*sin(x) where y is a constant
%                   y           Value for the parameter y in fx
%                   property    Desired thermodynamic property to solve for
%                               (T, p, v, u, h, s, c_p, c_v)
%                   target      Goal seek for f(x)=target
%                   lowLimit    Minimum value for the search interval
%                   highLimit   Maximum value for the search value
%                   tol         Tolerance for convergence
%                   k_max       Maximum number of iterations
% -------------------------------------------------------------------------
% Output:           result      The converged value of x which satisfies 
%                               f(x)=target
% -------------------------------------------------------------------------
% Example:          rho = bisection('region3_calc(x,y)',T,'p',p, 0, 1500, 0.000001, 10000);
%                   Region 3 of the IAPWS IF-97 is given in terms of
%                   f(rho,T). Because it is more reasonable to calculate
%                   the properties of water as a function of (p,T), the
%                   density which corresponds to the (p,T) point is needed.

% Reset the iteration counter
k=0;

% Force the while loop to initialize
guess = target+2*tol;

% Bisection method
while (abs(target - guess) > tol)
    x = (lowLimit + highLimit) / 2;
    tpwsData = eval(fx);
    guess = eval(['tpwsData.',property]);
    
    if (guess < target)
        lowLimit = x;
    else
        highLimit = x;
    end
 
    if k > k_max
        break;
    end 
    
    k=k+1;
    
end

% Return the converged value
result = x;