function [ engine_out ] = REVS_modify_engine2( engine_in, varargin )
%REVS_modify_engine takes in an engine structure and its displacement(L), apply engine modifiers,
%   and returns the result in a new engine structure.
%
%   [ engine_new ] = REVS_modify_engine( engine, varargin )
%
%   Avaiable engine modifiers can be applied by varargin as follows:
%   Scale engine: 'scale_factor', n
%               where n = NewSize/OrigSize


%% Parse varargs
% Check if enough varargin was passed
narginchk( 1, inf );


% Validate required arguments
validate_arg( engine_in, {'class_REVS_engine','class_REVS_engine_legacy','struct'});

% Engine Scaling Options
num_cylinders = parse_varargs( varargin, 'num_cylinders', [], 'numeric', {'integer', 'scalar'} );

% Scale engine size: scale_factor = NewSize/OrigSize (scale_factor > 0)
scale_factor = parse_varargs( varargin, 'scale_factor', [], 'numeric' );
scale_displacement = parse_varargs( varargin, 'scale_displacement', [], 'numeric' );
scale_max_power = parse_varargs( varargin, 'scale_max_power',[], 'numeric' );
scale_avg_power = parse_varargs( varargin, 'scale_avg_power',[], 'numeric' );
scale_max_torque = parse_varargs( varargin, 'scale_max_torque',[], 'numeric' );
scale_avg_torque = parse_varargs( varargin, 'scale_avg_torque',[], 'numeric' );

% window to use for averaging power or torque
averaging_window = parse_varargs( varargin, 'avg_window',[1800,6000] * convert.rpm2radps, 'numeric' );

scaling_knock_adjust = parse_varargs( varargin,'no_scaling_knock_adjust',true,'toggle');
scaling_friction_adjust = parse_varargs( varargin,'no_scaling_friction_adjust',true,'toggle');
scaling_sv_ratio_adjust = parse_varargs( varargin,'no_scaling_sv_ratio_adjust',true,'toggle');

% Fricion Adjustment
input_friction_adjust_Nm = parse_varargs( varargin, 'friction_adjust_Nm', 0.0, 'numeric');
input_friction_adjust_bar = parse_varargs( varargin, 'friction_adjust_bar', 0.0, 'numeric');


% Fill in if only window end points provided
if numel(averaging_window) ==2
	averaging_window = linspace(averaging_window(1), averaging_window(end), 25);
end




%% Copy and create a new engine structure
engine_out = engine_in;


%% Compute Deired Scaling

engine_torque = interp1( engine_in.full_throttle_speed_radps, engine_in.full_throttle_torque_Nm, averaging_window,'linear','extrap' );
engine_power = engine_torque .* averaging_window;

engine_average_torque = sum( engine_torque)./ length(averaging_window);
engine_average_power = sum( engine_power )./ length(averaging_window);


if ~isempty( scale_displacement )
	scale_factor = scale_displacement / engine_in.displacement_L;
elseif ~isempty( scale_max_power )
	scale_factor = scale_max_power / engine_in.max_power_W;
elseif ~isempty( scale_max_torque )
	scale_factor = scale_max_torque / engine_in.max_torque_Nm;
elseif ~isempty( scale_avg_power )
	scale_factor = scale_avg_power / engine_average_power;
elseif ~isempty( scale_avg_torque )
	scale_factor = scale_avg_torque / engine_average_torque;
elseif isempty(scale_factor) && ~isempty( num_cylinders)
	scale_factor = num_cylinders ./ engine_in.num_cylinders;
elseif isempty(scale_factor)
	scale_factor = 1.0;
end


%% Adjust Physical Characteristics
engine_out.displacement_L					= engine_in.displacement_L * (scale_factor);
engine_out.inertia_kgm2						= engine_in.inertia_kgm2 .* (scale_factor);

%% Set number of cylinders

disp_limit_cyl =     [	3,	  4,	  6,   	  8,	16];
disp_limit_max = 1.0*[1.5,	2.7,	4.5,	7.5,	15];
%disp_limit_max = [1.5,	2.7,	3.55,	7.5,	15.0];
disp_limit_min = 1.0*[0.6,	1.4,	2.5,	3.8,	8.0];
%disp_limit_min = [0.6,	1.0,	2.5,	3.56,	8.0];


if ~isempty( num_cylinders)
	% Use provided Cylinder Count
	engine_out.num_cylinders = num_cylinders;
else
	engine_out.num_cylinders = engine_in.num_cylinders;
	% Add cylinders is displacement above limit
	while engine_out.displacement_L	> interp1( disp_limit_cyl, disp_limit_max, engine_out.num_cylinders)
		engine_out.num_cylinders = engine_out.num_cylinders + 1 + (engine_out.num_cylinders >= 4);
	end
	% Remove cylinders if displacement below limit
	while engine_out.displacement_L	< interp1( disp_limit_cyl, disp_limit_min, engine_out.num_cylinders)
		engine_out.num_cylinders = engine_out.num_cylinders - 1 - (engine_out.num_cylinders >= 6);
	end
	
end



%% Generate cylinder scale factors
cyl_volume_scale_factor = scale_factor * engine_in.num_cylinders / engine_out.num_cylinders;
cyl_linear_scale_factor = cyl_volume_scale_factor ^ (1/3);

%% Adjust Cylinder Physical Characteristics
engine_out.bore_mm							= engine_in.bore_mm .* cyl_linear_scale_factor;
engine_out.stroke_mm						= engine_in.stroke_mm .* cyl_linear_scale_factor;	


if ~isempty( num_cylinders) || engine_out.num_cylinders ~= engine_in.num_cylinders
	engine_out.name = sprintf( '%s converted to %0.2fL %d Cylinder', engine_in.name, engine_out.displacement_L, engine_out.num_cylinders);
elseif scale_factor ~= 1
	engine_out.name = sprintf( '%s converted to %0.2fL', engine_in.name, engine_out.displacement_L);
end	 	
	
% Make meshgrid of speed and torque for later use
[speed_mesh_radps, brake_torque_mesh_Nm] = meshgrid( engine_in.fuel_map_speed_radps, engine_in.fuel_map_torque_Nm);
ind_torque_mesh_Nm = brake_torque_mesh_Nm - interp1( engine_in.closed_throttle_speed_radps, engine_in.closed_throttle_torque_Nm, speed_mesh_radps,'linear','extrap');


%% Friction Scaling

% SWRI: Friction Adjustment 0.1 bar FMEP for 10% Stroke increase - converted to torque
% Seemed big - try half
% 	scaling_friction_adjust_Nm =  scaling_friction_adjust * (cyl_stroke_scale_factor-1) * 100 * engine_out.displacement_L ./ (4 * pi );


% New A. Moskalik Method - Heywood Paper
FMEP_factors = [nan, nan, 41.60,41.02,nan, 41.38,nan, 40.97;  nan,nan,2.82,2.55,nan,2.34,nan, 2.21];
input_engine_FMEP_kPa =		[1,engine_in.num_cylinders/ engine_in.displacement_L]	* FMEP_factors( :,engine_in.num_cylinders) ; 
output_engine_FMEP_kPa =	[1,engine_out.num_cylinders/ engine_out.displacement_L] * FMEP_factors( :,engine_in.num_cylinders) ; 

scaling_friction_adjust_Nm =  scaling_friction_adjust * (output_engine_FMEP_kPa - input_engine_FMEP_kPa) * engine_out.displacement_L ./ (4 * pi );


%% Other Friciton Adjustment

total_friction_adjust_Nm = scaling_friction_adjust_Nm + input_friction_adjust_Nm + input_friction_adjust_bar * 100 * engine_out.displacement_L ./ (4 * pi );



%% Calculate BSFC adjustment for heat loss due to surface to volume effect

if scaling_sv_ratio_adjust 
	
	% ricardo_bsfc_sc_coef = [-0.1058 0.3458 -0.4564 (1.2168-0.000400000000000011)];
	% if cyl_disp_scale_factor <= 1
	%	bsfc_mult = polyval( [-0.1058 0.3458 -0.4564 (1.2168-0.000400000000000011)], cyl_disp_scale_factor);
	% else
	%	bsfc_mult = 1/polyval( [-0.1058 0.3458 -0.4564 (1.2168-0.000400000000000011)], 1/cyl_disp_scale_factor);
	% end
	
	% A. Moskalik first take from Novak & Blumberg SAE # 780943
	% bsfc_mult =  1 + scaling_sv_ratio_adjust *  -0.27* log(cyl_disp_scale_factor.^(1/3) ) ;

	% A Moskalik - makes it reversible
	% Original 
  	heat_loss_adjust =  (cyl_volume_scale_factor ^-(1/11.5));
	
	% Add Speed Rolloff
	heat_loss_adjust =  1+ (heat_loss_adjust -1 ) .* ( 1250 .* convert.rpm2radps ./ max(speed_mesh_radps,100) );

	
else 
	heat_loss_adjust = 1.0;
	
end


%% Adjust for increased knock when upscaling

if scaling_knock_adjust && cyl_volume_scale_factor > 1

	knock_compare_scale = cyl_volume_scale_factor ^(1/3);
	
	equiv_ind_torque_mesh_Nm = ind_torque_mesh_Nm .* knock_compare_scale;
	equiv_brake_torque_mesh_Nm = equiv_ind_torque_mesh_Nm + interp1( engine_in.closed_throttle_speed_radps, engine_in.closed_throttle_torque_Nm, speed_mesh_radps,'linear','extrap');
	equiv_fuel_gps = interp2( engine_in.fuel_map_speed_radps, engine_in.fuel_map_torque_Nm, engine_in.fuel_map_gps, speed_mesh_radps, equiv_brake_torque_mesh_Nm, 'spline' );
	
	knk_adjust_orig = equiv_fuel_gps ./ max(eps, ( engine_in.fuel_map_gps .* knock_compare_scale));
	knk_adjust_orig(1,:) = 1.0;
	knk_adjust_orig = max( knk_adjust_orig, 1.0);
	
	[min_consumption_speed_radps, min_consumption_torque_Nm] = REVS_calc_engine_min_bsfc(engine_in);
	approx_mbt_trq = interp1( min_consumption_speed_radps, min_consumption_torque_Nm, speed_mesh_radps, 'linear','extrap');
	
	wot_trq_Nm = interp1( engine_in.full_throttle_speed_radps, engine_in.full_throttle_torque_Nm,  speed_mesh_radps,  'linear','extrap' );
	
 	select_pts = knk_adjust_orig > 1.0 &  speed_mesh_radps < 3500 * convert.rpm2radps & equiv_brake_torque_mesh_Nm > approx_mbt_trq * 0.95 & equiv_brake_torque_mesh_Nm < wot_trq_Nm ;
% 	select_pts = knk_adjust_orig > 1.0 &  speed_mesh_radps < 3000 * convert.rpm2radps & equiv_brake_torque_mesh_Nm > approx_mbt_trq * 0.95 & equiv_brake_torque_mesh_Nm < wot_trq_Nm ;

   b = regress( knk_adjust_orig(select_pts), [ ones(sum(select_pts(:)),1),  equiv_ind_torque_mesh_Nm(select_pts), speed_mesh_radps(select_pts), equiv_ind_torque_mesh_Nm(select_pts) .* speed_mesh_radps(select_pts)]);
%   	b = regress( knk_adjust_orig(select_pts), [ ones(sum(select_pts(:)),1),  equiv_ind_torque_mesh_Nm(select_pts), speed_mesh_radps(select_pts), 1./max(speed_mesh_radps(select_pts),50), equiv_ind_torque_mesh_Nm(select_pts) .* speed_mesh_radps(select_pts)]);
	
	
	knk_adjust = [ ones(size(speed_mesh_radps(:))),  equiv_ind_torque_mesh_Nm(:),   speed_mesh_radps(:), equiv_ind_torque_mesh_Nm(:) .* speed_mesh_radps(:) ] * b;	
% 	knk_adjust = [ ones(size(speed_mesh_radps(:))),  equiv_ind_torque_mesh_Nm(:),  speed_mesh_radps(:), 1./max(speed_mesh_radps(:),50), equiv_ind_torque_mesh_Nm(:) .* speed_mesh_radps(:) ] * b;	

	
	knk_adjust = max(1.0, reshape( knk_adjust, size(speed_mesh_radps)));
	
	
%  	knk_adjust = knk_adjust_orig;
	
	
	% Fade out High Speed Effect (Extrapolation)
  	knk_adjust = 1 + (knk_adjust - 1) .* interp1([0,3000, 4500, 10000]*convert.rpm2radps,[1,1,0,0], speed_mesh_radps,'linear','extrap');

else 
	knk_adjust = 1.0;	
	
end


%% Scale Torque and adjust for fricion
engine_out.full_throttle_torque_Nm			= engine_in.full_throttle_torque_Nm  .* (scale_factor) - total_friction_adjust_Nm; %for torque curve
engine_out.closed_throttle_torque_Nm		= engine_in.closed_throttle_torque_Nm  .* (scale_factor) - total_friction_adjust_Nm; %for torque curve
engine_out.naturally_aspirated_torque_Nm	= engine_in.naturally_aspirated_torque_Nm  .* (scale_factor) - total_friction_adjust_Nm; %for torque curve	
   


engine_out.fuel_map_torque_Nm				= engine_in.fuel_map_torque_Nm .* (scale_factor) - total_friction_adjust_Nm;   %for BSFC map
engine_out.fuel_map_gps						= engine_in.fuel_map_gps .*(scale_factor).* heat_loss_adjust .* knk_adjust;			





end
