function [engine] = REVS_build_engine( engine_base,  fuel_data, varargin )
%REVS_build_engine Build REVS engine from source data
%
% Inputs:    engine,  map, WOT, CT, varargin
%
% engine - Skeleton engine (structure or class), must define displacement
%
% map - structure containing fuel map test data.  Must contain fields from
%		which speed torque and fuel can be obtained.  Available field names:
%		speed_radps
%		speed_rpm
%		torque_Nm
%		bmep_bar
%		fuel_gps
%		bsfc_gpkWhr
%
%
% WOT - optional structure defining maximum torque curve points or test data.
%		Naming scheme same as map.  If omitted use empty matrix [], and
%		convex hull of map data will be used to approximate curve.
%
% CT -  optional structure defining minimum torque curve points or test data.
%		Naming scheme same as map.  If omitted use empty matrix [], and
%		convex hull of map data will be used to approximate curve. If no fuel
%		map data is avaiable below 0 torque a default curve is scaled based on 
%		engine displacement.
%
% Varargin Options
%
% no_plots - true or false(default): if true, enineMap2Alpha will not make a
% call to makeREVSplots.m.
%

validate_arg( engine_base, {'class_REVS_engine','class_REVS_engine_legacy','struct'} );
validate_arg( fuel_data, {'cell','struct'} );



WOT_data = parse_varargs( varargin, 'WOT', [], {'cell','struct'});
CT_data = parse_varargs( varargin, 'CT', [], {'cell','struct'});



class_type = parse_varargs(varargin,'class','class_REVS_engine','char');

show_plots = parse_varargs(varargin,'no_plots',true,'toggle');
show_point_labels = parse_varargs(varargin,'no_point_labels',true,'toggle');
show_voronoi = parse_varargs(varargin,'show_voronoi',false,'toggle');
show_grid = parse_varargs(varargin,'show_grid',false,'toggle');

method = parse_varargs(varargin,'fit_method','gridfit');


% Scattered interpolant method
interp_method = parse_varargs(varargin,'interp_method','natural','char',{'natural','triangle','linear'});

% Gridfit smoothing
smooth = parse_varargs(varargin,'smoothing',0.0005,'numeric',{'scalar'});
x_smooth = parse_varargs(varargin,'x_smoothing',smooth,'numeric',{'scalar'});
y_smooth = parse_varargs(varargin,'y_smoothing',smooth,'numeric',{'scalar'});

fuel_out_speed_radps = parse_varargs( varargin,'out_speed',[],'numeric',{'vector'});
fuel_out_torque_Nm = parse_varargs( varargin,'out_torque',[],'numeric',{'vector'});

% Naturally Aspirated Curve
na_curve = parse_varargs( varargin,'na_curve',[]);


% Cylinder Deac Data
deac_data = parse_varargs( varargin,'deac',[]);

% Copy provided info into engine class
engine = eval( class_type );

if ~isfield(engine_base,'idle_speed_radps') || isa(engine_base.idle_speed_radps,'class_REVS_dynamic_lookup')
	% No need to modify idle
elseif isnumeric(engine_base.idle_speed_radps) && isscalar( engine_base.idle_speed_radps)
	temp = engine_base.idle_speed_radps;
	engine_base.idle_speed_radps = class_REVS_dynamic_lookup;
	engine_base.idle_speed_radps.axis_1.signal = 'veh_spd_mps';
	engine_base.idle_speed_radps.axis_1.breakpoints = [0,10];
	engine_base.idle_speed_radps.table = temp * ones(2,1);
else
	error('Unknown input for idle_speed_radps, must be a class_REVS_dynamic_lookup or scalar numeric value')
end
fields = fieldnames(engine_base);
for f = 1:length(fields)
	engine.(fields{f}) = engine_base.(fields{f});
end



%% Generate any missing fields in the provided input_data
% data needs  speed in radps, torque in Nm, fuel in g/sec



fuel_data = convert_and_combine( fuel_data, engine.displacement_L, 1);
WOT_data = convert_and_combine( WOT_data, engine.displacement_L, 0);
CT_data = convert_and_combine( CT_data, engine.displacement_L, 0);
deac_data = convert_and_combine( deac_data, engine.displacement_L, 1);


%% Construct Class to fill


if ~isempty(WOT_data)
    tolerance = max(WOT_data.torque_Nm) * 0.01;
else
    tolerance = max( fuel_data.torque_Nm) * 0.01;
end


%% Handle WOT Curve
[speed_radps, torque_Nm] = process_full_throttle( WOT_data, fuel_data, tolerance);

% Define Wide-Open Throttle (WOT) Torque Table
engine.full_throttle_speed_radps  = speed_radps;
engine.full_throttle_torque_Nm = torque_Nm;


%% Handle CTP Curve
[speed_radps, torque_Nm] = process_closed_throttle( CT_data, fuel_data, engine.displacement_L, max(engine.full_throttle_speed_radps),  tolerance);

% Define Closed Throttle Torque Table
engine.closed_throttle_speed_radps     = speed_radps;
engine.closed_throttle_torque_Nm       = torque_Nm;


%% Handle Naturally Aspirated Torque Curve
%TODO: add using raw MAP data from fuel map
%TODO: add passing in a curve
if max(engine.full_throttle_torque_Nm) >  15 / ( 4 * pi / (engine.displacement_L /1000)*1e-5 )
	% Assume Turbo if > 15 bar BMEP - set NA torque @ 9 bar BMEP
	engine.naturally_aspirated_speed_radps = engine.full_throttle_speed_radps([1,end]);
	engine.naturally_aspirated_torque_Nm = 9 / ( 4 * pi / (engine.displacement_L /1000)*1e-5 ) * [1,1];
	
	
else
	% Assume NA
	engine.naturally_aspirated_speed_radps = [];
	engine.naturally_aspirated_torque_Nm = [];
end

%% Construct Fuel Map

% Construct Fuel Map Grid
if isempty( fuel_out_speed_radps )
    fuel_out_speed_radps = find_clusters( fuel_data.speed_radps);
    fuel_out_speed_radps(end+1) = min(fuel_out_speed_radps)*0.5;
    fuel_out_speed_radps(end+1) = 0;
    fuel_out_speed_radps(end+1) = max(engine.full_throttle_speed_radps);
    fuel_out_speed_radps = fill_gaps(fuel_out_speed_radps, max(fuel_out_speed_radps)/15);
end

if isempty(fuel_out_torque_Nm)
    fuel_out_torque_Nm = find_clusters( fuel_data.torque_Nm);
    fuel_out_torque_Nm(end+1) = 1.1 * min(engine.closed_throttle_torque_Nm);
    fuel_out_torque_Nm(end+1) = 1.05* max(engine.full_throttle_torque_Nm);
    fuel_out_torque_Nm = fill_gaps(fuel_out_torque_Nm, max(fuel_out_torque_Nm)/20);
end

% Add Closed Throttle Points
CT_interp.speed_radps = fuel_out_speed_radps(:);
CT_interp.torque_Nm = interp1( engine.closed_throttle_speed_radps, engine.closed_throttle_torque_Nm, CT_interp.speed_radps,'linear','extrap');

% Combine Lists of Points
fuel_pts.speed_radps = [fuel_data.speed_radps(:); ];
fuel_pts.torque_Nm = [fuel_data.torque_Nm(:); ];
fuel_pts.fuel_gps = [fuel_data.fuel_gps(:); ];


% interpolate map scatter points and grid to generate a surface matrix
speed_scale = max( fuel_data.speed_radps) * 1.0;
torque_scale = max(fuel_data.torque_Nm) * 1.0;
fuel_per_cyc_gridfit = scatter2surf(fuel_pts.speed_radps, fuel_pts.torque_Nm, fuel_pts.fuel_gps./max(fuel_pts.speed_radps, 0.001) , fuel_out_speed_radps, fuel_out_torque_Nm,'method','gridfit','interp','bicubic','xscale',speed_scale,'yscale',torque_scale,'xsmooth',x_smooth,'ysmooth',y_smooth);
fuel_per_cyc_interp = scatter2surf(fuel_pts.speed_radps, fuel_pts.torque_Nm, fuel_pts.fuel_gps./max(fuel_pts.speed_radps, 0.001), fuel_out_speed_radps, fuel_out_torque_Nm,'method','scatterinterp','interp',interp_method,'xscale',speed_scale,'yscale',torque_scale);

% Meshgrid for output
[fuel_mesh_speed_radps, fuel_mesh_torque_Nm] = meshgrid( fuel_out_speed_radps, fuel_out_torque_Nm );

if strcmpi( method,'gridfit')
    map_out_fuel_gps = fuel_per_cyc_gridfit .* fuel_mesh_speed_radps;
elseif strcmpi( method,'interp' )
    map_out_fuel_gps = fuel_per_cyc_interp .* fuel_mesh_speed_radps;
end

% fuel_gridfit = scatter2surf(fuel_pts.speed_radps, fuel_pts.torque_Nm, fuel_pts.fuel_gps , fuel_out_speed_radps, fuel_out_torque_Nm,'method','gridfit','interp','bicubic','xscale',speed_scale,'yscale',torque_scale,'xsmooth',x_smooth,'ysmooth',y_smooth);
% fuel_interp = scatter2surf(fuel_pts.speed_radps, fuel_pts.torque_Nm, fuel_pts.fuel_gps, fuel_out_speed_radps, fuel_out_torque_Nm,'method','scatterinterp','interp',interp_method,'xscale',speed_scale,'yscale',torque_scale);
% 
% % Meshgrid for output
% [fuel_mesh_speed_radps, fuel_mesh_torque_Nm] = meshgrid( fuel_out_speed_radps, fuel_out_torque_Nm );
% 
% if strcmpi( method,'gridfit')
%     map_out_fuel_gps = fuel_gridfit ;
% elseif strcmpi( method,'interp' )
%     map_out_fuel_gps = fuel_interp;
% end

% assign the grid vectors to the surface X and Y definitions.
engine.fuel_map_speed_radps = fuel_out_speed_radps;
engine.fuel_map_torque_Nm = fuel_out_torque_Nm;
engine.fuel_map_gps = max(map_out_fuel_gps, 0); % change negative values to 0



%% Handle Cylinder Deac

if ~isempty(deac_data)
    
    deac_fuel_per_cyc_gridfit = scatter2surf(deac_data.speed_radps, deac_data.torque_Nm, deac_data.fuel_gps./max(deac_data.speed_radps, 0.001), fuel_out_speed_radps, fuel_out_torque_Nm,'method','gridfit','interp','bicubic','xscale',speed_scale,'yscale',torque_scale,'xsmooth',x_smooth,'ysmooth',y_smooth);
    deac_fuel_per_cyc_interp = scatter2surf(deac_data.speed_radps, deac_data.torque_Nm, deac_data.fuel_gps./max(deac_data.speed_radps, 0.001), fuel_out_speed_radps, fuel_out_torque_Nm,'method','scatterinterp','interp',interp_method,'xscale',speed_scale,'yscale',torque_scale);
    
    if strcmpi( method,'gridfit')
        deac_out_fuel_gps = deac_fuel_per_cyc_gridfit .* fuel_mesh_speed_radps;
    elseif strcmpi( method,'interp' )
        deac_out_fuel_gps = deac_fuel_per_cyc_interp .* fuel_mesh_speed_radps;
    end
    
%     deac_fuel_gridfit = scatter2surf(deac_data.speed_radps, deac_data.torque_Nm, deac_data.fuel_gps, fuel_out_speed_radps, fuel_out_torque_Nm,'method','gridfit','interp','bicubic','xscale',speed_scale,'yscale',torque_scale,'xsmooth',x_smooth,'ysmooth',y_smooth);
%     deac_fuel_interp = scatter2surf(deac_data.speed_radps, deac_data.torque_Nm, deac_data.fuel_gps, fuel_out_speed_radps, fuel_out_torque_Nm,'method','scatterinterp','interp',interp_method,'xscale',speed_scale,'yscale',torque_scale);
%     
%     if strcmpi( method,'gridfit')
%         deac_out_fuel_gps = deac_fuel_gridfit ;
%     elseif strcmpi( method,'interp' )
%         deac_out_fuel_gps = deac_fuel_interp;
%     end

    % Fuel map for cylinder deac [g/sec]
    % Always Less than base map
    engine.deac_fuel_map_gps = max(min(deac_out_fuel_gps,engine.fuel_map_gps), 0.0);
    
    % Compute deac region
    k = convhull(deac_data.speed_radps, deac_data.torque_Nm);
    deac_hull_speed_radps = deac_data.speed_radps(k);
    deac_hull_torque_Nm = deac_data.torque_Nm(k);
    next_pt = [2:length(k), 1];
    
    % Compute Centroid
    deac_cent_speed_radps = sum( (deac_hull_speed_radps + deac_hull_speed_radps(next_pt)) .* ( deac_hull_speed_radps  .* deac_hull_torque_Nm(next_pt) - deac_hull_speed_radps(next_pt) .* deac_hull_torque_Nm) )/6  /( sum( deac_hull_speed_radps  .* deac_hull_torque_Nm(next_pt) - deac_hull_speed_radps(next_pt) .* deac_hull_torque_Nm)/2);
    deac_cent_torque_Nm = sum( (deac_hull_torque_Nm + deac_hull_torque_Nm(next_pt)) .* ( deac_hull_speed_radps  .* deac_hull_torque_Nm(next_pt) - deac_hull_speed_radps(next_pt) .* deac_hull_torque_Nm) )/ 6 / ( sum( deac_hull_speed_radps  .* deac_hull_torque_Nm(next_pt) - deac_hull_speed_radps(next_pt) .* deac_hull_torque_Nm)/2);
    
    % Expand by 25%
    deac_expand_factor = 1.25;
    deac_expand_speed_radps = deac_cent_speed_radps + deac_expand_factor * (deac_hull_speed_radps - deac_cent_speed_radps);
    deac_expand_torque_Nm = deac_cent_torque_Nm + deac_expand_factor * (deac_hull_torque_Nm - deac_cent_torque_Nm);
    
    % Expand to next breakpoint in fuel map for interpolation
%     deac_mask_base = inpolygon( fuel_mesh_speed_radps, fuel_mesh_torque_Nm, deac_expand_speed_radps, deac_expand_torque_Nm);
%     deac_mask = zeros( size(deac_mask_base) + 2);
%     base_pts_1 = 2:(size(deac_mask,1)-1);
%     base_pts_2 = 2:(size(deac_mask,2)-1);
%     shift = [0,0;0,1;1,1;1,0;1,-1;0,-1;-1,-1;-1,0;-1,1];
%     for i = 1:size(shift,1)
%         deac_mask( base_pts_1 + shift(i,1), base_pts_2 + shift(i,2)) = deac_mask( base_pts_1 + shift(i,1), base_pts_2 + shift(i,2)) | deac_mask_base;
%     end
%     deac_mask = deac_mask( base_pts_1, base_pts_2 );
    
    % Mask Outside excess extrapolation
	deac_mask = inpolygon( fuel_mesh_speed_radps, fuel_mesh_torque_Nm, deac_expand_speed_radps, deac_expand_torque_Nm);
    engine.deac_fuel_map_gps(~deac_mask) = nan;
    
end




%% Plot Results for Examination

if show_plots
    
    if show_voronoi
        figure;
        [vx, vy] = voronoi(fuel_data.speed_radps/speed_scale,fuel_data.torque_Nm/torque_scale);
        superplot(fuel_data.speed_radps * convert.radps2rpm,fuel_data.torque_Nm,'r+7');
        ax = axis;
        hold on
        plot(vx*speed_scale* convert.radps2rpm,vy*torque_scale,'b-');
        axis(ax);
        title( 'Voronoi Diagram');
        xlabel('Speed (RPM)');
        ylabel('Torque (Nm)');
    end
    
    if show_grid
        
        % Show Grid
        REVS_plot_engine(engine,'fuel flow', 'show_bmep','show_CTP_curve','no_min_bsfc','no_power_lines');
        hold on
        superplot( fuel_data.speed_radps * convert.radps2rpm, fuel_data.torque_Nm,'rro4');
        
        for i = 1:length(engine.fuel_map_speed_radps)
            superplot(convert.radps2rpm * 	engine.fuel_map_speed_radps(i),[-intmax,intmax],'gy-1');
        end
        
        for i = 1:length(engine.fuel_map_torque_Nm)
            superplot([-intmax,intmax],engine.fuel_map_torque_Nm(i),'gy-1');
        end
              
        for i = 1:length(fuel_data.fuel_gps)
            text( fuel_data.speed_radps(i) * convert.radps2rpm+ 35, fuel_data.torque_Nm(i), num2str(fuel_data.fuel_gps(i),'%03.2f' ));
        end
                
    end
    
    % Plot BSFC
    if length(unique(fuel_data.bsfc_gpkWhr)) < 20
        REVS_plot_engine_bsfc(engine,'show_CTP_curve','no_min_bsfc','show_bmep','deac_prop',0.0,'contours',unique(fuel_data.bsfc_gpkWhr));
    else
        REVS_plot_engine_bsfc(engine,'show_CTP_curve','no_min_bsfc','show_bmep','deac_prop',0.0);
    end
    
    
    if ~isempty(WOT_data)
        superplot( WOT_data.speed_radps * convert.radps2rpm, WOT_data.torque_Nm,'or.');
    end
    
    superplot( engine.full_throttle_speed_radps * convert.radps2rpm, engine.full_throttle_torque_Nm ,'bo10-1');
    
    if ~isempty(CT_data)
        superplot( CT_data.speed_radps * convert.radps2rpm, CT_data.torque_Nm,'or.');
    end
    
    superplot( engine.closed_throttle_speed_radps * convert.radps2rpm, engine.closed_throttle_torque_Nm ,'bo10-1');
    
    datasets = unique(fuel_data.dataset)';
    data_colors = {'k','lg','m','db','dg','c'};
    line_handles = [];
    for d = datasets
        pts = fuel_data.dataset == d;
        line_handles(end+1) = superplot( fuel_data.speed_radps(pts) * convert.radps2rpm, fuel_data.torque_Nm(pts), [data_colors{d},'o3']);
        if show_point_labels
            text( fuel_data.speed_radps(pts) * convert.radps2rpm, fuel_data.torque_Nm(pts), num2str(fuel_data.bsfc_gpkWhr(pts),'%03.0f' ));
        end
    end
    
    legend( line_handles, strcat( {'Dataset '},int2str(datasets')));
    
	%% PLot fuel per injection
	REVS_plot_engine(engine,'fuel per injection','show_CTP_curve','no_min_bsfc','show_bmep','deac_prop',0.0);

	
	
	
	
    %% Plotting for cylinder deac
    if ~isempty(deac_data)
        % Plot BSFC
        REVS_plot_engine_bsfc(engine,'show_CTP_curve','no_min_bsfc','show_bmep');
        
        % 	superplot(WOT.pts_speed_radps  * convert.radps2rpm,	WOT.pts_torque_Nm,'or.');
        superplot( engine.full_throttle_speed_radps * convert.radps2rpm, engine.full_throttle_torque_Nm ,'bo10-1');
        
        
        % 	superplot(CT.pts_speed_radps  * convert.radps2rpm,	CT.pts_torque_Nm,'or.');
        superplot( engine.closed_throttle_speed_radps * convert.radps2rpm, engine.closed_throttle_torque_Nm ,'bo10-1');
        
        line_handles = [];
        datasets = unique(deac_data.dataset)';
        for d = datasets
            pts = deac_data.dataset == d;
            line_handles(end+1) = superplot( deac_data.speed_radps(pts) * convert.radps2rpm, deac_data.torque_Nm(pts), [data_colors{d},'o3']);
            if show_point_labels
                text( deac_data.speed_radps(pts) * convert.radps2rpm, deac_data.torque_Nm(pts), num2str(deac_data.bsfc_gpkWhr(pts),'%03.0f' ));
            end
        end
        
        legend( line_handles, strcat( {'Dataset '},int2str(datasets')));
        
        
        superplot( deac_hull_speed_radps * convert.radps2rpm , deac_hull_torque_Nm, 'gy--');
        superplot( deac_expand_speed_radps * convert.radps2rpm , deac_expand_torque_Nm, 'gy--');
        
    end
    
    
    
    
end



end


function [speed_radps, torque_Nm] = process_full_throttle( WOT, map, tol)

if ~isempty(WOT) && length(WOT.speed_radps ) > 50
    % > 50 points => Test Data need to reduce
    
    % Add point at origin
    speed_radps = [WOT.speed_radps(:); 0];
    torque_Nm = [WOT.torque_Nm(:); 0];
    
    % Filter provided data by averaging into bins
    [speed_radps,torque_Nm] =  gen_avg_curve(speed_radps, torque_Nm, 50 * convert.rpm2radps);
    
    % Begin constructing curve by selecting points with largest error
    [speed_radps,torque_Nm] = build_curve(speed_radps,torque_Nm, tol);
    
    
elseif ~isempty(WOT) && length(WOT.speed_radps ) > 5
    % Use what was provided
    [speed_radps, sort_idx] = unique(WOT.speed_radps);
    torque_Nm = WOT.torque_Nm(sort_idx);
    
else
    %Combine map data with any points provided
    
    if isempty( WOT)
        xscat = [max(map.speed_radps)*1.1;map.speed_radps;0];
        yscat = [-1;map.torque_Nm;0];
    else
        xscat = [max(map.speed_radps)*1.1;map.speed_radps;0;WOT.speed_radps];
        yscat = [-1;map.torque_Nm;0;WOT.torque_Nm];
    end
    
    % Find WOT from speed and load convex hull
    
    
    k = convhull(xscat, yscat);
    
    speed_radps = xscat(k);
    torque_Nm = yscat(k);
    
    speed_radps = speed_radps( torque_Nm >= 0 );
    torque_Nm = torque_Nm( torque_Nm >= 0 );
    
    [speed_radps, sort_idx] = unique(speed_radps);
    torque_Nm = torque_Nm(sort_idx);
    
end

% Cleanup zero speed
if torque_Nm(1) ~= 0
    speed_radps = [0;speed_radps(:)];
    torque_Nm = [0;torque_Nm(:)];
elseif speed_radps(1) ~= 0
    speed_radps(1) = 0;
end

% Roll off on top end
if torque_Nm(end) > 0
    speed_radps = [speed_radps(:); speed_radps(end) .* [1.05; 1.15] ];
    torque_Nm = interp1( speed_radps(1:end-2), torque_Nm(:), speed_radps, 'linear','extrap' );
    torque_Nm(end-1) = min( torque_Nm(end-1), 0.90*torque_Nm(end-2 ) * speed_radps(end-2) / speed_radps(end-1) );		% Force downward slope / No Increase in power
    torque_Nm(end) = min(torque_Nm(end), 0.0 );								% Force end to Zero If not already
end


end



function [speed_radps, torque_Nm] = process_closed_throttle( CT, map, displacement_L, max_speed_radps, tol)

if ~isempty(CT)
    
    % Filter provided data by averaging into bins
    [speed_radps,torque_Nm] =  gen_avg_curve(CT.speed_radps, CT.torque_Nm, 50 * convert.rpm2radps);
    
    % Begin constructing curve by selecting points with largest error
    [speed_radps,torque_Nm] = build_curve(speed_radps,torque_Nm, tol);
    
elseif sum( map.torque_Nm < 0 ) > 5
    % Negative Torques were mapped - Use Convex Hull
    xscat = [max(map.speed_radps);map.speed_radps;0];
    yscat = [min(map.torque_Nm);map.torque_Nm;1];
    
    k = convhull(xscat, yscat);
    speed_radps = xscat(k);
    torque_Nm = yscat(k);
    
    speed_radps = speed_radps( torque_Nm < 0 );
    torque_Nm = torque_Nm( torque_Nm < 0 );
    
    % Filter provided data by averaging into bins
    [speed_radps,torque_Nm] =  gen_avg_curve(speed_radps, torque_Nm, 50 * convert.rpm2radps);
    
    % Begin constructing curve by selecting points with largest error
    [speed_radps,torque_Nm] = build_curve(speed_radps,torque_Nm, tol);
    
else
    
    % Scale "default" closed throttle curve by engine engine.displacement_L
    speed_radps = convert.rpm2radps*[657.2 812.1 1050 1991 3008 4064 4555];
    torque_Nm = [-29.57 -34.01 -36.12 -40.44 -47.98 -57.72 -62.77] * displacement_L / 3.3;
    
end


% Extend from 0 to max WOT speed
speed_extend_radps = unique([0; speed_radps(:); max_speed_radps]);
torque_Nm = interp1( speed_radps, torque_Nm, speed_extend_radps,'linear','extrap');
speed_radps = speed_extend_radps;


end






function out = convert_and_combine(in, displacement_L, has_fuel)

if isempty(in)
    out = [];
    return;
end

if ~iscell(in)
    in = {in};
end


out.speed_radps = [];
out.torque_Nm = [];
out.dataset = [];

if has_fuel
    out.fuel_gps = [];
    out.bsfc_gpkWhr = [];
end



for i = 1:length(in)
    
    
    % Generate speed in radps if missing
    if ~isfield( in{i},'speed_radps') && isfield(in{i}, 'speed_rpm' )
        in{i}.speed_radps = in{i}.speed_rpm * convert.rpm2radps;
    end
    
    % Generate torque in Nm if missing
    if ~isfield( in{i},'torque_Nm') && isfield( in{i},'bmep_bar')
        in{i}.torque_Nm = in{i}.bmep_bar / ( 4 * pi / (displacement_L /1000)*1e-5 );
	end
    
	
	if has_fuel
	
		% Serialize fuel if table provided
		if isfield( in{i},'fuel_gps') && ~isvector( in{i}.fuel_gps)
			[in{i}.fuel_gps, in{i}.torque_Nm, in{i}.speed_radps] = table2columns(  in{i}.fuel_gps, in{i}.torque_Nm, in{i}.speed_radps );
		end

		% Serialize BSFC if table provided
		if isfield( in{i},'bsfc_gpkWhr') && ~isvector( in{i}.bsfc_gpkWhr)
			[in{i}.bsfc_gpkWhr, in{i}.torque_Nm, in{i}.speed_radps] = table2columns(  in{i}.bsfc_gpkWhr, in{i}.torque_Nm, in{i}.speed_radps );
		end

		% Serialize BTE if table provided
		if isfield( in{i},'efficiency_pct') && ~isvector( in{i}.efficiency_pct)
			[in{i}.efficiency_pct, in{i}.torque_Nm, in{i}.speed_radps] = table2columns(  in{i}.efficiency_pct, in{i}.torque_Nm, in{i}.speed_radps );
		end



		% GEnerate fuel in g/sec if missing
		if isfield( in{i},'fuel_gps') 
			% Have fuel data directly
		elseif isfield( in{i},'bsfc_gpkWhr')
			in{i}.fuel_gps = in{i}.bsfc_gpkWhr .* in{i}.speed_radps  .* in{i}.torque_Nm /1000 / 3600;
		elseif isfield( in{i},'efficiency_pct')
			in{i}.fuel_gps = (in{i}.speed_radps  .* in{i}.torque_Nm) ./( in{i}.efficiency_pct / 100 .* in{i}.fuel_prop.energy_density_MJpkg *1000 );
		end

		% GEnerate BSFC if missing
		if ~isfield( in{i},'bsfc_gpkWhr') 
			in{i}.bsfc_gpkWhr = in{i}.fuel_gps ./ ( in{i}.speed_radps  .* in{i}.torque_Nm /1000 / 3600 );
		end
    
	end
    
    % Remove Nans
    nan_idx = isnan(in{i}.torque_Nm) | isnan(in{i}.speed_radps) ;
    
    if has_fuel && isfield( in{i},'fuel_gps')
        nan_idx = nan_idx | isnan(in{i}.fuel_gps);
        in{i}.fuel_gps(nan_idx) = [];
		in{i}.bsfc_gpkWhr(nan_idx) = [];
    end
    
    in{i}.speed_radps(nan_idx) = [];
    in{i}.torque_Nm(nan_idx) = [];
    
    
    out.speed_radps = [out.speed_radps; in{i}.speed_radps(:)];
    out.torque_Nm = [out.torque_Nm; in{i}.torque_Nm(:)];
    out.dataset = [out.dataset; i * ones(size(in{i}.speed_radps(:)))];
    
    if has_fuel
        out.fuel_gps = [out.fuel_gps; in{i}.fuel_gps(:)];
        out.bsfc_gpkWhr = [out.bsfc_gpkWhr; in{i}.bsfc_gpkWhr(:)];
    end
    
end


end


function [spd, trq] = gen_avg_curve( spd_in, trq_in, step)

spd = 0:step:max(spd_in);
trq = zeros(size(spd));

i = length( spd);
while i > 0
    pts = abs(spd_in - spd(i)) < step ;
    
    if sum(pts) <= 0
        trq(i) = [];
        spd(i) = [];
    else
        trq(i) = mean( trq_in(pts ) );
        spd(i) = mean( spd_in(pts ) );
        i = i-1;
    end
end
end

function out = fill_gaps( data, max_gap )

	data = unique(data(:));

	big_gap_left = diff(data) > max_gap;
	big_gap_left(end+1) = false; % Tack on end to match size
	big_gap_right = circshift(big_gap_left, 1);
	
	add_points = (data(big_gap_left) + data(big_gap_right))/2;
	
	out = unique([data; add_points(:)]);
	
end


function [spd, trq] = build_curve(spd_in, trq_in, tol )

% Begin constructing curve by selecting points with largest error
spd = [spd_in(1);spd_in(end);];
trq = [trq_in(1);trq_in(end);];

max_err = inf;
while max_err > tol
    
    err = trq_in - interp1(spd, trq, spd_in,'linear');
    
    [max_err, idx] = max( abs(err));
    
    spd(end+1) = spd_in(idx);
    trq(end+1) = trq_in(idx);
    
    [spd, sort_idx] = unique(spd);
    trq = trq(sort_idx);
    
end

end