function REVS_plot_engine(e, plot_type, varargin)
% % REVS_plot_engine(e, format, varargin)
% no_figure                   = parse_varargs(varargin,'no_figure',false,'toggle');
% contour_levels              = parse_varargs(varargin,'contours',[],'numeric',{'vector'});
% contour_filled              = parse_varargs(varargin,'filled',false,'toggle');
% contours_labels             = parse_varargs(varargin,'no_contour_labels',true,'toggle');
% contour_width               = parse_varargs(varargin,'contour_width',1,'numeric');
% contour_color               = parse_varargs(varargin,'contour_color','auto');
% contour_cmap                = parse_varargs(varargin,'colormap','jet');
% 
% plot_power_lines            = parse_varargs(varargin,'no_power_lines',true,'toggle');
% plot_min_bsfc               = parse_varargs(varargin,'no_min_bsfc',true,'toggle');
% plot_na_curve               = parse_varargs(varargin,'no_na_curve',true,'toggle');
% plot_ctp_curve			  = parse_varargs(varargin,'show_ctp_curve',false,'toggle');
% plot_deac_norm			  = parse_varargs(varargin,'deac_prop',1.0,'numeric',{'scalar'});
% 
% highlight_patch             = parse_varargs(varargin,'highlight_contour',[],'numeric');
% highlight_color             = parse_varargs(varargin,'highlight_color',[0.50,.85,0.50],'numeric');
% highlight_alpha             = parse_varargs(varargin,'highlight_alpha',[],'numeric',{'scalar'});
% 
% line_width				  = parse_varargs(varargin,'line_width',3,'numeric');
% fill_outside				  = parse_varargs(varargin,'show_above_wot',true,'toggle');
% 
% plot_audit                  = parse_varargs(varargin,'plot_audit',false,'');
% plot_audit_points           = parse_varargs(varargin,'plot_audit_points',0,'numeric');
% 
% plot_test_data              = parse_varargs(varargin,'plot_test_data',false,'');
% plot_test_data_points       = parse_varargs(varargin,'plot_test_data_points',0,'numeric',{'initger'});
% 
% show_bmep					  = parse_varargs(varargin,'show_bmep',false,'toggle');
% show_matrix_min_speeds      = parse_varargs(varargin,'show_matrix_min_speeds',false,'toggle');
% idle_plot                   = parse_varargs(varargin,'idle_plot',false,'toggle');
% title_str                   = parse_varargs(varargin,'title','','string');

validate_arg(e, {'class_REVS_engine','class_REVS_engine_legacy','struct'});
validate_arg(plot_type, {'char'}, {'efficiency','fuel flow', 'fuel consumption','fuel per injection','BSFC','brake specific fuel','BS CO2', 'brake specific CO2','CO2'});


no_figure                   = parse_varargs(varargin,'no_figure',(false),'toggle');
contour_levels              = parse_varargs(varargin,'contours',[],'numeric',{'vector'});
contour_filled              = parse_varargs(varargin,'filled',false,'toggle');
contours_labels             = parse_varargs(varargin,'no_contour_labels',true,'toggle');
contour_width               = parse_varargs(varargin,'contour_width',1,'numeric');
contour_color               = parse_varargs(varargin,'contour_color','auto');
contour_cmap                = parse_varargs(varargin,'colormap','jet','colormap');

plot_power_lines            = parse_varargs(varargin,'no_power_lines',true,'toggle');
plot_min_bsfc               = parse_varargs(varargin,'no_min_bsfc',true,'toggle');
plot_na_curve               = parse_varargs(varargin,'no_na_curve',true,'toggle');
plot_ctp_curve				= parse_varargs(varargin,'show_ctp_curve',false,'toggle');
plot_deac_norm				= parse_varargs(varargin,'deac_prop',1.0,'numeric',{'scalar'});

highlight_patch             = parse_varargs(varargin,'highlight_contour',[],'numeric');
highlight_color             = parse_varargs(varargin,'highlight_color',[0.50,.85,0.50],'numeric');
highlight_alpha             = parse_varargs(varargin,'highlight_alpha',[],'numeric',{'scalar'});

line_width					= parse_varargs(varargin,'line_width',3,'numeric');
fill_outside				= parse_varargs(varargin,'show_above_wot',true,'toggle');

plot_test_data              = parse_varargs(varargin,'plot_test_data',false);
plot_test_data_points       = parse_varargs(varargin,'plot_test_data_points',0,'numeric',{'integer'});
plot_test_time_range        = parse_varargs(varargin,'time_range',[]);

show_bmep					= parse_varargs(varargin,'show_bmep',false,'toggle');
show_matrix_min_speeds      = parse_varargs(varargin,'show_matrix_min_speeds',false,'toggle');
idle_plot					= parse_varargs(varargin,'idle_plot',false,'toggle');
title_str                   = parse_varargs(varargin,'title','','string');



%% Create Fine Mesh for making smooth images
speeds_rpm = linspace(min(e.fuel_map_speed_radps(:) ), max(e.fuel_map_speed_radps(:)), 200 ) .* convert.radps2rpm;
torques_Nm = linspace(min(e.fuel_map_torque_Nm(:)),max(e.fuel_map_torque_Nm(:)),220);

[mesh_speed_rpm, mesh_torque_Nm] = meshgrid( speeds_rpm, torques_Nm);

% [in_mesh_speed_radps, in_mesh_torque_Nm] = meshgrid(e.fuel_map_speed_radps, e.fuel_map_torque_Nm);
% mesh_fuel_gps = interp2( e.fuel_map_speed_radps * convert.radps2rpm, e.fuel_map_torque_Nm, e.fuel_map_gps ./ ( 5000 + in_mesh_speed_radps .* in_mesh_torque_Nm) , mesh_speed_rpm, mesh_torque_Nm) .* ( 5000 + mesh_speed_rpm .* convert.rpm2radps .*  mesh_torque_Nm);


% Apply clinderdeac proportionately
if ( ~isfield( e,'deac_fuel_map_gps') && ~isprop( e,'deac_fuel_map_gps') )|| all(isnan(e.deac_fuel_map_gps(:))) 
	title_deac_append = '';
	fuel_data = e.fuel_map_gps;
else
	title_deac_append  = sprintf('- %2d%% Cylinder Deac', plot_deac_norm*100);
	fuel_data = plot_deac_norm .* e.deac_fuel_map_gps + ( 1-plot_deac_norm ) .* e.fuel_map_gps;
	fuel_data = min(fuel_data,  e.fuel_map_gps);	% Take Care of Nans
end

mesh_fuel_gps = interp2( e.fuel_map_speed_radps * convert.radps2rpm, e.fuel_map_torque_Nm, fuel_data , mesh_speed_rpm, mesh_torque_Nm);

switch upper(plot_type)
	case {'EFF','EFFICIENCY'}
		plot_data = 100 * ((mesh_speed_rpm .* mesh_torque_Nm) ./ 9539 ) ./ ( mesh_fuel_gps * e.fuel.energy_density_MJpkg );
		plot_data(mesh_torque_Nm < 0) = nan;
		contour_levels =set_if_empty(contour_levels,[5,10, 15, 20, 23, 26:2:30, 32:50]);
		title_plot_type_append = 'Brake Thermal Efficiency ( % )';
	case {'FUEL FLOW', 'FUEL CONSUMPTION'}
		plot_data =  mesh_fuel_gps;
		increment = max(0.5,round(max(plot_data(:))/30));
		contour_levels =set_if_empty(contour_levels,[0:increment:max(plot_data(:))]);
		title_plot_type_append = 'Fuel Flow ( g/s )';
    case {'FUEL PER INJECTION'}
		plot_data =  mesh_fuel_gps *60 ./ mesh_speed_rpm / e.num_cylinders / 2*1000;
		max_contour = ceil(max(plot_data(2:end,end)));
		increment = round(max(0.5,max_contour/30));
		contour_levels =set_if_empty(contour_levels,0:increment:max_contour);
		title_plot_type_append = 'Fuel Flow ( mg / inj )';
	case {'BSFC','BRAKE SPECIFIC FUEL'}
		plot_data = (mesh_fuel_gps .* 3600 .* 9549) ./ max(eps,mesh_speed_rpm .* mesh_torque_Nm);
		plot_data(mesh_torque_Nm < 0) = nan;
		title_plot_type_append = 'Brake Specific Fuel Consumption ( g/kW*hr )';
		contour_levels =set_if_empty(contour_levels,[200:10:250, 275:25:375, 400,450,500,600,750,1000]);
    case {'CO2'}
        plot_data =  mesh_fuel_gps .* e.fuel.carbon_weight_fraction * (44.0095/12.0107);
        increment = max(0.5,round(max(plot_data(:))/30));
		contour_levels =set_if_empty(contour_levels,[0:increment:max(plot_data(:))]);
		title_plot_type_append = 'CO_2 Emissions ( g/s )';      
    case {'BS CO2', 'BRAKE SPECIFIC CO2'}     
        plot_data = (mesh_fuel_gps .* e.fuel.carbon_weight_fraction * (44.0095/12.0107) *  3600 .* 9549) ./ max(eps,mesh_speed_rpm .* mesh_torque_Nm);
        plot_data(mesh_torque_Nm < 0) = nan;
        title_plot_type_append = 'Brake Specific CO_2 ( g/kW*hr )';
		contour_levels =set_if_empty(contour_levels,[600:25:750, 800:50:900, 1000, 1250, 1500, 2000, 3000,5000]); 
    otherwise
		error('Unknown plot format %s',plot_type)
end

% = max(max(mesh_bsfc));


if no_figure
	fig = gcf;
	ax = gca;
else
	fig = figure;
	ax = axes;
end

if ~isempty(highlight_patch)
	c = contourc(speeds_rpm, torques_Nm, plot_data , contour_levels);
	idx = 1;
	contour_size = size(c,2);
	while idx < contour_size
		if c(1,idx) == highlight_patch
			num_pts = c(2,idx);
			x_pts = c(1,idx+(1:num_pts));
			y_pts = c(2,idx+(1:num_pts));
			
% 			ft = interp1(e.full_throttle_speed_radps*convert.radps2rpm, e.full_throttle_torque_Nm, x_pts) + 0.05*max(e.full_throttle_torque_Nm);
% 			y_pts = min( y_pts, ft);
			
			h = patch(x_pts,y_pts,highlight_color);
			if ~isempty( highlight_alpha )
				alpha(h,highlight_alpha);
			end
			
		end
		idx = idx + c(2,idx) +1;
	end
	hold on
end


for s = 1:length(speeds_rpm)
	ft = interp1(e.full_throttle_speed_radps*convert.radps2rpm, e.full_throttle_torque_Nm, speeds_rpm(s)) + 0.05*max(e.full_throttle_torque_Nm);
	for t = 1:length(torques_Nm)
		if mesh_torque_Nm(t,s) > ft
			plot_data(t,s) = NaN;
		else
			plot_data(t,s) = plot_data(t, s);
		end
	end
end
% Additional Contour Plotting Arguments
addtl_args = { 'linewidth',contour_width, 'color', contour_color};


% if isstruct(plot_audit)
% 	newhist2d(plot_audit.mech.engine.crankshaft_rpm,linspace(500,6500,50),plot_audit.mech.engine.crankshaft_gross_Nm, linspace(0,250,50),plot_audit.mech.engine.crankshaft_gross_kW.*data_pts.*(0.5+0.5*(plot_audit.mech.engine.crankshaft_rpm > 1000))*6, plot_audit.time, min(contour_levels), max(contour_levels));
% 	hold on
%     	
% 	if plot_audit_points
% 		% plot(ax,speed_data(1:plot_audit_points:end), torque_data(1:plot_audit_points:end),'r.','MarkerSize',1);
%         
%         speed_data =plot_audit.mech.engine.crankshaft_rpm(data_pts);
%         torque_data = plot_audit.mech.engine.crankshaft_gross_Nm(data_pts);
%         plothg(speed_data(1:plot_audit_points:end), torque_data(1:plot_audit_points:end),'r.','MarkerSize',1);
% 
% 	end
% 	
% else
    
if isa(plot_test_data,'class_test_data')
    
    

    if ~isempty( plot_test_time_range)
        data_pts = (plot_test_data.time >= plot_test_time_range(1)) & (plot_test_data.time <= plot_test_time_range(end));
    else
        data_pts = true(size(plot_test_data.time));
    end
        
	[~,ch] = newhist2d(plot_test_data.engine.speed_rpm, linspace(500,6500,50), plot_test_data.engine.crankshaft_torque_Nm, linspace(0,e.max_torque_Nm,50),plot_test_data.engine.crankshaft_power_kW.*data_pts.*(0.5+0.5*(plot_test_data.engine.speed_rpm > 1000))*6, plot_test_data.time, min(contour_levels), max(contour_levels));
	hold on

	if plot_test_data_points
%		plot(plot_test_data.engine.speed_rpm(1:plot_test_data_points:end), plot_test_data.engine.crankshaft_torque_Nm(1:plot_test_data_points:end),'r.','MarkerSize',1);
        
        speed_data = plot_test_data.engine.speed_rpm(data_pts);
        torque_data = plot_test_data.engine.crankshaft_torque_Nm(data_pts);
        plothg(speed_data(1:plot_test_data_points:end), torque_data(1:plot_test_data_points:end),'r.','MarkerSize',1);
	end
end





if contour_filled
	[c,h] = contourf(ax, mesh_speed_rpm, mesh_torque_Nm, plot_data , contour_levels ,addtl_args{:});
else
	[c,h] = contour(ax, mesh_speed_rpm, mesh_torque_Nm, plot_data , contour_levels ,addtl_args{:});
end

%assignin('base','h',h)
cmap = contour_cmap;
if ischar( cmap)
	cmap = eval(cmap);
end


% Linearize colormap based on contour levels
contour_levels_norm = (contour_levels - min(contour_levels)) ./ ( max(contour_levels) - min(contour_levels));
cmap_convert = interp1( contour_levels_norm, linspace(0,1,length(contour_levels_norm)), linspace(0,1,size(cmap,1)));
cmap = interp1( linspace(0,1,size(cmap,1)),cmap, cmap_convert);

% cmap_convert = linspace(0,1,length(contour_levels));
% cmap_pts = interp1((1:size(cmap,1))/size(cmap,1),cmap,cmap_convert);

colormap(cmap);

if contours_labels
	clabel(c,h,'labelspacing',144,'color',[0.5,0.5,0.5],'fontweight','bold')
end

hold on

% Draw Max NA Torque
if plot_na_curve
	plot(ax, e.naturally_aspirated_speed_radps * convert.radps2rpm, e.naturally_aspirated_torque_Nm,'k--','linewidth',line_width)
end

% Draw Minimum BSFC Line
if plot_min_bsfc
	[min_consumption_radps, min_consumption_torque_Nm, ~] = REVS_calc_engine_min_bsfc(e);
	h = plot(ax,min_consumption_radps * convert.radps2rpm, min_consumption_torque_Nm, '-','color',[1.0, 0.5, 0.8], 'linewidth', 3);
end





if fill_outside
	%     fill([plot_size(2); plot_size(2); 0; 0; e.full_throttle_speed_radps(:) * convert.radps2rpm],[0; plot_size(4)*2; plot_size(4)*2; 0; e.full_throttle_torque_Nm(:)],[0.8 0.8 1.0])
%	fill([ realmax; realmax; 0; 0; e.full_throttle_speed_radps(:) * convert.radps2rpm; e.full_throttle_speed_radps(end) * convert.radps2rpm; ],[0; realmax; realmax; 0; e.full_throttle_torque_Nm(:); 0],[0.8 0.8 1.0])
%	fill([ realmax; 0; 0; e.closed_throttle_speed_radps(:) * convert.radps2rpm; realmax ],[ -realmax; -realmax; e.closed_throttle_torque_Nm(1); e.closed_throttle_torque_Nm(:); e.closed_throttle_torque_Nm(end)],[0.8 0.8 1.0])
	fill([ 10000; 10000; 0; 0; e.full_throttle_speed_radps(:) * convert.radps2rpm; e.full_throttle_speed_radps(end) * convert.radps2rpm; ],[0; 10000; 10000; 0; e.full_throttle_torque_Nm(:); 0],[0.8 0.8 1.0])
	fill([ 10000; 0; 0; e.closed_throttle_speed_radps(:) * convert.radps2rpm; 10000 ],[ -10000; -10000; e.closed_throttle_torque_Nm(1); e.closed_throttle_torque_Nm(:); e.closed_throttle_torque_Nm(end)],[0.8 0.8 1.0])
end

% Draw line @ 0 - just because
plot([-10000, 10000],[0,0],'k:','linewidth',line_width);


% Draw CTP Torque
plot(ax, e.closed_throttle_speed_radps * convert.radps2rpm, e.closed_throttle_torque_Nm,'k-','linewidth',line_width)

% Draw WOT Torque
plot(ax, e.full_throttle_speed_radps * convert.radps2rpm, e.full_throttle_torque_Nm,'k-','linewidth',line_width)

% Set Extents
if ~isfield( e,'idle_speed_radps') && ~isprop( e ,'idle_speed_radps' )
	x_min = 0;
elseif isstruct(e.idle_speed_radps) && isfield( e.idle_speed_radps,'table')
	x_min =min(e.idle_speed_radps.table(:))*0.8  *convert.radps2rpm;
elseif isnumeric( e.idle_speed_radps ) && isscalar( e.idle_speed_radps )
	x_min = e.idle_speed_radps*0.8*convert.radps2rpm;
else
	x_min = 0;
end

x_max = e.full_throttle_speed_radps(end) * convert.radps2rpm;
y_min = plot_ctp_curve * 1.1 * min(e.closed_throttle_torque_Nm);
y_max = 1.05 *max(e.full_throttle_torque_Nm);


xlim( [x_min,x_max]);
ylim( [y_min,y_max]);

if plot_power_lines
	max_pow = max( e.full_throttle_speed_radps .* e.full_throttle_torque_Nm );
	power_line_increment = max(5,5* round( max_pow/(1000*8) / 5));
% 	powerlines(x_min, x_max, y_min,y_max,'increment',power_line_increment)
	powerlines;
end

	if show_bmep
		ax2 = axes('position',[0.08,0.11, 0.001, 0.815]);
		set(ax2,'YColor','b');
		set(ax2,'Color',get(fig,'color'));
		ylabel('{\color{black}Torque ( Nm )}       {\color{blue}BMEP ( Bar )}');
		ylim(ax2, [y_min;y_max] *  4 * pi / (e.displacement_L *100) );
		set(ax2,'XTickLabel','');
				
		addlistener(ax,'YLim','PostSet', @(src,evt)bmep_axis_callback(src,evt, ax,ax2, 4 * pi / (e.displacement_L *100) ) );
		set(ax2,'HandleVisibility','off');
		
	else
		ylabel(ax,'Torque ( Nm )');
	end


% Apply Labels
xlabel(ax,'Speed ( RPM )');

if (~isprop(e,'name') && ~isfield( e, 'name' )) || isempty(e.name)
	e.name = 'Engine';
end

if isempty(title_str)
    title_str = sprintf('%s %s\n%s',strrep(e.name,'_',' '), title_deac_append, title_plot_type_append);
end
 title(title_str);

if idle_plot
    idle_fuel_gps = interp2( e.fuel_map_speed_radps, e.fuel_map_torque_Nm, e.fuel_map_gps , e.idle_speed_radps.table(1) * ones(size(torques_Nm)), torques_Nm);
    
    fplothg(torques_Nm, idle_fuel_gps);
        
    xyt([e.name ' Torque at Idle, ' num2str(e.idle_speed_radps.table(1) * convert.radps2rpm) ' RPM'],'Fuel Rate (g/s)');
    text(-40,1.9,['Disp [L]: ' num2str(e.displacement_L)]);
    text(-40,1.8125,['Inertia [kg-m^2]: ' num2str(e.inertia_kgm2)]);
    text(-40,1.7,['Speed [RPM]: ' num2str(e.idle_speed_radps.table(1) * convert.radps2rpm)]);
    text(-40,1.6,['Motoring Torque [Nm]: ' num2str(interp1(e.closed_throttle_speed_radps, e.closed_throttle_torque_Nm,e.idle_speed_radps.table(1),'linear','extrap'))]);
    
%    legend off;
    ylim([0 2]);
    xlim([-50 250]);
    
    eff = (torques_Nm * e.idle_speed_radps.table(1) / 10)./(idle_fuel_gps * e.fuel.energy_density_MJpkg);
    plothg(torques_Nm, eff/25,'g.-');
    plot(torques_Nm(eff == max(eff)), max(eff)/25, 'kx');
    text(torques_Nm(eff == max(eff)), max(eff)/25 + .075, num2str(max(eff)));
    
%     xyt([e.name ' Torque at Idle, ' num2str(e.idle_speed_radps.table(1) * convert.radps2rpm) ' RPM'],'Efficiency (%)');
%     text(-40,47,['Disp [L]: ' num2str(e.displacement_L)]);
%     text(-40,45,['Inertia [kg-m^2]: ' num2str(e.inertia_kgm2)]);
%     text(-40,43,['Speed [RPM]: ' num2str(e.idle_speed_radps.table(1) * convert.radps2rpm)]);
%     text(-40,41,['Motoring Torque [Nm]: ' num2str(interp1(e.closed_throttle_speed_radps, e.closed_throttle_torque_Nm,e.idle_speed_radps.table(1),'linear','extrap'))]);
%     vlineat(interp1(e.closed_throttle_speed_radps, e.closed_throttle_torque_Nm,e.idle_speed_radps.table(1),'linear','extrap'),'r--');
%     vlineat(interp1(e.full_throttle_speed_radps, e.full_throttle_torque_Nm,e.idle_speed_radps.table(1),'linear','extrap'),'r--');
%     legend off;
%     ylim([0 50]);
%     xlim([-50 250]);

    vlineat(interp1(e.closed_throttle_speed_radps, e.closed_throttle_torque_Nm,e.idle_speed_radps.table(1),'linear','extrap'),'r--');
    vlineat(interp1(e.full_throttle_speed_radps, e.full_throttle_torque_Nm,e.idle_speed_radps.table(1),'linear','extrap'),'r--');

    legend('Fuel Rate','Eff. div 25','Torque Limits');

end

if show_matrix_min_speeds
    vlineat(REVS_calc_engine_power_speed(e, 0.12, 'min') * convert.radps2rpm,'r');
    vlineat(REVS_calc_engine_power_speed(e, 0.18, 'min') * convert.radps2rpm,'r');
end

end

function bmep_axis_callback(src,evt, ax, ax2, scale)
	
	ylim(ax2, ylim(ax) *  scale);
end

function out = set_if_empty( in, empty_val)

if isempty(in)
	out = empty_val;
else
	out = in;
end

end