function REVS_plot_engine_diff(e1, e2, compare_mode, varargin)
%function REVS_plot_engine_efficiency(e, varargin)
%
% no_figure                   = parse_varargs(varargin,'no_figure',false,'toggle');
% contour_levels              = parse_varargs(varargin,'contours',[0,10, 15, 20, 23, 25:50],'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,'line_width',1,'numeric');
% contour_color               = parse_varargs(varargin,'line_color','auto');
% power_lines                 = parse_varargs(varargin,'no_power_lines',true,'toggle');
% cmap                        = parse_varargs(varargin,'colormap','jet');
% highlight_patch             = parse_varargs(varargin,'highlight_contour',[],'numeric');
% highlight_color             = parse_varargs(varargin,'highlight_color',[0,.6,0],'numeric');
% highlight_alpha             = parse_varargs(varargin,'highlight_alpha',[],'numeric',{'scalar'});
% full_throttle_width         = parse_varargs(varargin,'line_width',3,'numeric');
% fill_above_full_throttle    = parse_varargs(varargin,'showabovefullthrottle',true,'toggle');
% plot_audit                  = parse_varargs(varargin,'plot_audit',false,'');


if nargin < 3
	compare_mode				= 'efficiency';
else
	validatestring(compare_mode,{'efficiency','bsfc','bsfc_pct','gps','gps_pct','fuel','fuel_pct', 'fuel_vol','fuel_vol_pct','co2','co2_pct' });
end


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,'line_width',1,'numeric');
contour_color               = parse_varargs(varargin,'line_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');
show_bmep					= parse_varargs(varargin,'show_bmep',false,'toggle');

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

full_throttle_width         = parse_varargs(varargin,'line_width',3,'numeric');

fill_above_full_throttle    = parse_varargs(varargin,'showabovefullthrottle',true,'toggle');

correct_fuel				= parse_varargs(varargin,'no_fuel_adjust',true,'toggle');


scaling_mode				= parse_varargs(varargin,'scaling','avg_power','char',{'none','max_power','max_torque','avg_torque','avg_power'});
scaling_window				= parse_varargs( varargin,'scaling_window',linspace(1800, 5500, 25)* convert.rpm2radps, 'numeric' );
scaling_fuel				= parse_varargs( varargin,'scaling_no_fuel_adjust', true, 'toggle');



time_range					= parse_varargs(varargin,'time_range',[]);

heatmap_data				= parse_varargs(varargin,'heatmap',[]);
heatmap_integrate			= parse_varargs(varargin,'heatmap_type','','char');
heatmap_time_range			= parse_varargs(varargin,'heatmap_time_range',time_range);


data_points_data			= parse_varargs(varargin,'data_points',[]);
data_points_decimation		= parse_varargs(varargin,'data_points_decimation',1);
data_points_time_range		= parse_varargs(varargin,'data_points_time_range',time_range);

% 
% 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',{'integer'});
% 
% 
% plot_test_time_range        = 

title_str                   = parse_varargs(varargin,'title','','string');



% Get statistics on engine 1
e1_torque = max(0,interp1( e1.full_throttle_speed_radps, e1.full_throttle_torque_Nm, scaling_window,'linear','extrap' ));
e1_average_torque = sum(e1_torque )./ length(scaling_window);
e1_average_power = sum(e1_torque .* scaling_window )./ length(scaling_window);


if correct_fuel
	e2.fuel_map_gps = e2.fuel_map_gps * e2.fuel.energy_density_MJpkg / e1.fuel.energy_density_MJpkg;
end

opts = {};

if ~scaling_fuel
% 	opts{end+1} = 'no_scale_fuel_adjust'; % For old REVS_modify_engine
	opts{end+1} = 'no_scaling_knock_adjust';
	opts{end+1} = 'no_scaling_friction_adjust';
	opts{end+1} = 'no_scaling_sv_ratio_adjust';
end


% Scale 2nd engine if needed
switch lower(scaling_mode)
	
	case 'max_power'
		e2 = REVS_modify_engine2( e2, 'scale_max_power', e1.max_power_kW, opts{:});
	case 'max_torque'
		e2 = REVS_modify_engine2( e2, 'scale_max_torque', e1.max_torque_Nm, opts{:});
	case 'avg_power'
		e2 = REVS_modify_engine2( e2, 'scale_avg_power', e1_average_power, 'avg_window', scaling_window, opts{:});
	case 'avg_torque'
		e2 = REVS_modify_engine2( e2, 'scale_avg_torque', e1_average_torque, 'avg_window', scaling_window, opts{:});	
		
end

speeds_rpm = linspace(min(min(e1.fuel_map_speed_radps),min(e2.fuel_map_speed_radps)) .* convert.radps2rpm, max(max(e1.fuel_map_speed_radps), max(e2.fuel_map_speed_radps)) .* convert.radps2rpm, 200 );
torques_Nm = linspace(min(min(e1.fuel_map_torque_Nm), min(e2.fuel_map_torque_Nm)),max(max(e1.fuel_map_torque_Nm), max(e2.fuel_map_torque_Nm)),250);

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


%Compute Gridded Maps for comparison
[e1_gps, e1_eff, e1_bsfc] = interp_engine( e1, mesh_speed_rpm, mesh_torque_Nm );
[e2_gps, e2_eff, e2_bsfc] = interp_engine( e2, mesh_speed_rpm, mesh_torque_Nm );


switch lower(compare_mode)
	case 'efficiency'
		plot_data =  e1_eff - e2_eff;
		contour_scale = 1.0;
        title_compare_str = 'Efficiency Difference';
	case 'bsfc'
		plot_data =  e1_bsfc - e2_bsfc;
		contour_scale = 10.0;
        title_compare_str = 'BSFC Difference';
	case {'gps','fuel'}
		plot_data =  e1_gps - e2_gps;
		contour_scale = 0.1;
        title_compare_str = 'Fuel Mass Flow (g/s) Difference';
	case 'bsfc_pct'
		plot_data = 100 * (e1_gps - e2_gps) ./ e2_gps;
		contour_scale = 1.0;
        title_compare_str = 'BSFC % Difference';
	case {'gps_pct','fuel_pct'}
		plot_data = 100 * (e1_gps - e2_gps) ./ e2_gps;
		contour_scale = 1.0;
        title_compare_str = 'Fuel Mass Flow % Difference';
    case 'fuel_vol'
		plot_data = e1_gps./e1.fuel.density_kgpL_15C - e2_gps./e2.fuel.density_kgpL_15C;
		contour_scale = 0.1;
        title_compare_str = 'Fuel Volume Flow (mL/s) Difference';
	case 'fuel_vol_pct'
		plot_data = 100*(e1_gps./e1.fuel.density_kgpL_15C - e2_gps./e2.fuel.density_kgpL_15C) ./ (e2_gps./e2.fuel.density_kgpL_15C);
		contour_scale = 1.0;
        title_compare_str = 'Fuel Volume Flow % Difference';
    case 'co2'
		plot_data = (e1_gps.*e1.fuel.carbon_weight_fraction - e2_gps.*e2.fuel.carbon_weight_fraction).*(44.0095/12.0107);
		contour_scale = 0.1;
        title_compare_str = 'CO2 Emissions (g/s) Difference';
	case 'co2_pct'
		plot_data = 100*(e1_gps.*e1.fuel.carbon_weight_fraction - e2_gps.*e2.fuel.carbon_weight_fraction)./(e2_gps.*e2.fuel.carbon_weight_fraction);
		contour_scale = 1.0;   
        title_compare_str = 'CO2 Emissions % Difference';
    otherwise
		error( 'Unknown comparison mode');
		
end



if isempty(contour_levels)
	contour_levels = [-20,-15, -10:-1, -0.5, -0.2, 0, 0.2, 0.5, 1:10, 15,20] * contour_scale;
end

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));
			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

% Additional Contour Plotting Arguments
addtl_args = { 'linewidth',contour_width, 'color', contour_color};

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


colormap(contour_cmap);

if contours_labels
	clabel(c,h,'labelspacing',144)
end

hold on

if fill_above_full_throttle	
	% 	plot_size = axis;	
	max_x = 2* max(speeds_rpm);
	max_y = 2* max(torques_Nm);
	e1_ft = interp1(e1.full_throttle_speed_radps * convert.radps2rpm,  e1.full_throttle_torque_Nm, speeds_rpm(:) );
	e2_ft = interp1(e2.full_throttle_speed_radps * convert.radps2rpm,  e2.full_throttle_torque_Nm, speeds_rpm(:) );
	
	fill([max_x; max_x; 0; 0; speeds_rpm(:) ],[0; max_y; max_y; 0; max(e1_ft, e2_ft)],[0.8 0.8 1.0])	
end

if plot_na_curve
	plot(ax, e1.naturally_aspirated_speed_radps * convert.radps2rpm, e1.naturally_aspirated_torque_Nm,'r--','linewidth',full_throttle_width)
	plot(ax, e2.naturally_aspirated_speed_radps * convert.radps2rpm, e2.naturally_aspirated_torque_Nm,'b--','linewidth',full_throttle_width)
end

if plot_min_bsfc
	[min_consumption_radps, min_consumption_torque_Nm, ~] = REVS_calc_engine_min_bsfc(e1);
	plot(min_consumption_radps * convert.radps2rpm, min_consumption_torque_Nm, 'm-', 'linewidth', 3);
	[min_consumption_radps, min_consumption_torque_Nm, ~] = REVS_calc_engine_min_bsfc(e2);
	plot(min_consumption_radps * convert.radps2rpm, min_consumption_torque_Nm, 'c-', 'linewidth', 3);
end

% Full throttle Curves
plot(ax, e1.full_throttle_speed_radps * convert.radps2rpm, e1.full_throttle_torque_Nm,'r-','linewidth',full_throttle_width)
plot(ax, e2.full_throttle_speed_radps * convert.radps2rpm, e2.full_throttle_torque_Nm,'b-','linewidth',full_throttle_width)

% Closed throttle Curves
plot(ax, e1.closed_throttle_speed_radps * convert.radps2rpm, e1.closed_throttle_torque_Nm,'r-','linewidth',full_throttle_width)
plot(ax, e2.closed_throttle_speed_radps * convert.radps2rpm, e2.closed_throttle_torque_Nm,'b-','linewidth',full_throttle_width)


xlim([0.8*min([e1.idle_speed_radps.table(:);e2.idle_speed_radps.table(:)]), max( e1.full_throttle_speed_radps(end), e2.full_throttle_speed_radps(end) ) ] .* convert.radps2rpm);
ylim([1.05* min([e1.closed_throttle_torque_Nm(:);e2.closed_throttle_torque_Nm(:)]), 1.05* max([e1.full_throttle_torque_Nm(:);e2.full_throttle_torque_Nm(:)])]);

% Plot heat map overlay
if ~isempty(heatmap_data)
	% Nothing to do
	if isa(heatmap_data,'class_test_data')
		
		% select variables from class_test_data
		heatmap_torque = heatmap_data.engine.crankshaft_torque_Nm;
		heatmap_speed = heatmap_data.engine.speed_rpm;
		heatmap_time = heatmap_data.time;
		
		if isempty( heatmap_integrate )
			heatmap_value = heatmap_data.engine.crankshaft_power_kW;
		else
			heatmap_value = eval(['heatmap_data.engine.', heatmap_integrate]);
		end
						
	elseif isstruct(heatmap_data) && isfield(heatmap_data,'mech') && isfield( heatmap_data.mech, 'engine')
		
		% select variables form audit structure
		heatmap_torque = heatmap_data.mech.engine.crankshaft_torque_Nm;
		heatmap_speed = heatmap_data.mech.engine.speed_rpm;
		heatmap_time = heatmap_data.time;
		
		if isempty( heatmap_integrate )
			heatmap_value = heatmap_data.mech.engine.crankshaft_gross_kW;
		else
			heatmap_value = eval(['heatmap_data.mech.engine.', heatmap_integrate]);
		end

		
		
		% 	h = newhist2d(plot_audit.mech.engine.crankshaft_rpm,linspace(500,6500,50),plot_audit.mech.engine.crankshaft_gross_Nm, linspace(0,e1.max_torque_Nm,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));
		% 	alpha(0.5);
		% 	colormap([1 1 1; jet])
		% 	speed_data =plot_audit.mech.engine.crankshaft_rpm(data_pts);
		% 	torque_data = plot_audit.mech.engine.crankshaft_gross_Nm(data_pts);
		%
		% 	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
		%
		%
		% 	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
	
	% Remove energy from excluded region
	if ~isempty( heatmap_time_range)
		heatmap_value( heatmap_time < heatmap_time_range(1) ) = 0;
		heatmap_value( heatmap_time > heatmap_time_range(end) ) = 0;
	end
	
	% Reduce weighting of idle
	heatmap_value = heatmap_value .* (0.5+0.5*(heatmap_speed > 1000));	
			
	caxis_min = interp1([2,length(jet)+1], [min(contour_levels), max(contour_levels)], 1,'linear','extrap');
	caxis_max = max(contour_levels);
	
	[~,ch] = newhist2d(heatmap_speed, linspace(500,6500,40), heatmap_torque, linspace(e1.min_torque_Nm,e1.max_torque_Nm,30), heatmap_value, heatmap_time,caxis_min, caxis_max);
	alpha(ch, 0.75);
	colormap([1,1,1; jet])
	caxis([caxis_min, caxis_max]);
	
end

% Plot operational data point cloud
if ~isempty( data_points_data)

	
	if isa(data_points_data,'class_test_data')
		
		% select variables from class_test_data
		data_points_torque = data_points_data.engine.crankshaft_torque_Nm;
		data_points_speed = data_points_data.engine.speed_rpm;

		if ~isempty( data_points_label )
			data_points_value = data_points_data.engine.crankshaft_power_kW;
		end
		
				
	elseif isstruct(data_points_data) && isfield(data_points_data,'mech') && isfield( data_points_data.mech, 'engine')
		
		% select variables form audit structure
		data_points_torque = data_points_data.mech.engine.crankshaft_torque_Nm;
		data_points_speed = data_points_data.mech.engine.speed_rpm;

		if ~isempty( data_points_label )
			data_points_value = data_points_data.mech.engine.crankshaft_power_kW;
		end
	
	end

	
	pts_idx = 1:data_points_decimation:length( data_points_speed);
	
	plothg( data_points_speed(pts_idx), data_points_torque(pts_idx), 'r.','MarkerSize',1 )
	if ~isempty( data_points_value )
		
	end
	
	
end


if isempty(e1.name)
	e1.name = inputname(1);
end

if isempty(e2.name)
	e2.name = inputname(2);
end

if isempty( title_str)

if strcmpi( e1.fuel.id, e2.fuel.id)
	fuel_adjust_str = '';
elseif correct_fuel
	fuel_adjust_str = sprintf('\nAdjusted For Fuel Heating Value');
else
	fuel_adjust_str = sprintf('\nNot Adjusted For Fuel Heating Value');	
end
	
% if strcmpi(scaling_mode,'none')
% 	scale_str = '';
% else
% 	scale_str = sprintf( ' scaled to %0.3f L', e2.displacement_L);
% end	

% REVS_modify_engine2 appends scaling information to engine name
scale_str = '';


title_str = sprintf('%s\n {\\color{red}%s} - {\\color{blue}%s}%s%s',title_compare_str, strrep(e1.name,'_',' '), strrep(e2.name,'_',' '), scale_str, fuel_adjust_str);

end


title(title_str);

if plot_power_lines
	plot_size = axis;
	powerlines
end

xlabel('Speed ( RPM )');


	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, ylim(ax) *  4 * pi / (e1.displacement_L *100) );
		set(ax2,'XTickLabel','');
				
		addlistener(ax,'YLim','PostSet', @(src,evt)bmep_axis_callback(src,evt, ax,ax2, 4 * pi / (e1.displacement_L *100) ) );
		set(ax2,'HandleVisibility','off');
		
	else
		ylabel(ax,'Torque ( Nm )');
	end

end


function [gps, eff, bsfc] = interp_engine( e, mesh_speed_rpm, mesh_torque_Nm )

gps = interp2( e.fuel_map_speed_radps * convert.radps2rpm, e.fuel_map_torque_Nm, e.fuel_map_gps , mesh_speed_rpm, mesh_torque_Nm);

% Nan out points above WOT
ft = interp1(e.full_throttle_speed_radps*convert.radps2rpm, e.full_throttle_torque_Nm, mesh_speed_rpm) + 0.05*max(e.full_throttle_torque_Nm);
gps( mesh_torque_Nm > ft) = nan;

bsfc = gps * 3600 ./ max(eps, (mesh_speed_rpm .* mesh_torque_Nm) ./ 9539 );


eff = 100 * ((mesh_speed_rpm .* mesh_torque_Nm) ./ 9539 ) ./ ( gps * e.fuel.energy_density_MJpkg );


end

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



