function [min_consumption_speed_radps, min_consumption_torque_Nm, min_consumption_power_W] = REVS_calc_engine_min_bsfc(engine, varargin)
%function [min_consumption_radps, min_consumption_torque_Nm, min_consumption_power_W] = REVS_calc_engine_min_bsfc(engine, varargin)
%
%  max_kW      = parse_varargs(varargin, 'max_kw', max_radps * maxtorque_Nm / 1000, 'numeric');
%  z_offset    = parse_varargs(varargin, 'z_offset', 0, 'numeric');
%  increment   = parse_varargs(varargin, 'increment', 10, 'numeric');
%  do_plots    = parse_varargs(varargin, 'do_plots', 0, 'toggle');

    min_radps = min(engine.full_throttle_speed_radps);
    max_radps = max(engine.full_throttle_speed_radps);
    mintorque_Nm = min(engine.full_throttle_torque_Nm);
    maxtorque_Nm = max(engine.full_throttle_torque_Nm);
    
    max_kW      = parse_varargs(varargin, 'max_kw', max_radps * maxtorque_Nm / 1000, 'numeric');
    z_offset    = parse_varargs(varargin, 'z_offset', 0, 'numeric');
    increment   = parse_varargs(varargin, 'increment', 2.5, 'numeric');
    do_plots    = parse_varargs(varargin, 'do_plots', 0, 'toggle');
    

    
    if do_plots
        REVS_plot_engine_bsfc(engine);
        hold on;
    end
    
    min_radps = max(eps, min_radps);
        
    speeds_radps = linspace(min_radps, max_radps, 1000);
	max_torque_Nm = interp1(engine.full_throttle_speed_radps, engine.full_throttle_torque_Nm, speeds_radps); % Preinterpolate for speed
	
	
	min_consumption_power_W = 1000 * max( (0:increment:max_kW), 1);
	min_consumption_flow_gps = zeros( size( min_consumption_power_W ) );
	min_consumption_torque_Nm = zeros( size( min_consumption_power_W ) );
	min_consumption_speed_radps = zeros( size( min_consumption_power_W) );
	
	
%    for power = [5 10:1:max_kW];
    for idx = 1:length(min_consumption_power_W)
	
		torques_Nm = min_consumption_power_W(idx)./speeds_radps;

		fuel_rate_gps =  interp2(engine.fuel_map_speed_radps, engine.fuel_map_torque_Nm, engine.fuel_map_gps, speeds_radps, torques_Nm);
		
	
		% Remove Bad Points
		fuel_rate_gps(torques_Nm > max_torque_Nm) = Inf;		% Above WOT Curve
		if idx > 1
% 			fuel_rate_gps( speeds_radps < min_consumption_speed_radps(idx-1) ) = Inf;	% Force Increasing with speed
		end
		
		
		[min_consumption_flow_gps(idx), min_idx] = min( fuel_rate_gps);
		min_consumption_speed_radps(idx) = speeds_radps( min_idx);
		

	end

	
	% Remove Inf Points
	inf_pts = isinf( min_consumption_flow_gps ) ;
    min_consumption_speed_radps(inf_pts) = [];
	min_consumption_power_W(inf_pts) = [];
	

	for  idx = 1:length(min_consumption_speed_radps)-1	
		min_consumption_speed_radps(idx) = min( min_consumption_speed_radps(idx:end) );	
	end
	
	%Force Unique speeds (Monotonic)
    [ min_consumption_speed_radps, unique_pts] = unique(min_consumption_speed_radps,'last');
	min_consumption_power_W = min_consumption_power_W(unique_pts);
	min_consumption_torque_Nm = min_consumption_power_W ./ min_consumption_speed_radps;

	
	
    if do_plots
        plothg(min_consumption_speed_radps * convert.radps2rpm, min_consumption_torque_Nm, 'm.-', 'linewidth', 3);
        hold off;
    end
    
    
end