function [ max_wheeltorque_Nm, max_wheelpower_kW, max_etorque_Nm, etorques_Nm, wheeltorques_Nm ] = REVS_calc_max_wheeltorque( engine, transmission, drive_axle, varargin )
%function [ max_wheeltorque_Nm, max_wheelpower_kW, max_etorque_Nm, etorques_Nm, wheeltorques_Nm ] = REVS_calc_max_wheeltorque( engine, transmission, drive_axle, varargin )

stall_radps     = parse_varargs(varargin, 'stall_radps', 0, 'numeric');
redline_radps   = parse_varargs(varargin, 'redline_radps', [], 'numeric');
do_plots        = parse_varargs(varargin, 'do_plots', false, 'toggle');
top_speed_mph   = parse_varargs(varargin, 'max_speed_mph', 60, 'numeric');
top_speed_mps   = parse_varargs(varargin, 'max_speed_mps', [], 'numeric');

if isempty(redline_radps)
    redline_radps =  min(engine.full_throttle_speed_radps(end) - 500 * convert.rpm2radps, REVS_calc_engine_power_speed(engine, 0.98, 'max'));
end

if ~isempty(top_speed_mps)
    top_speed_mph = top_speed_mps * convert.mps2mph;
end

v_mph = 1:0.25:top_speed_mph;

v_mps               = v_mph * convert.mph2mps;
wheelspeed_radps    = v_mps / drive_axle.tire.radius_m;
gbo_radps           = wheelspeed_radps * drive_axle.final_drive.gear_ratio;
% gbo_rpm = gbo_radps * convert.radps2rpm;

if transmission.type ~= enum_transmission_type.CVT
    gbi_radps = max(stall_radps, kron(gbo_radps, transmission.gear.ratio(2:end)'));
    gbi_radps(gbi_radps > redline_radps) = 0;   % ignore points above redline
    
    etorques_Nm = interp1(engine.full_throttle_speed_radps, engine.full_throttle_torque_Nm, gbi_radps, 'linear', 0);
    
    wheeltorques_Nm = etorques_Nm .* repmat(transmission.gear.ratio(2:end)', 1, length(gbo_radps)) * drive_axle.final_drive.gear_ratio;
    if stall_radps > 0
        wheeltorques_Nm(1,:) = wheeltorques_Nm(1,:) .* interp1([0 stall_radps 1e6], [1.8, 1, 1], gbo_radps .* transmission.gear.ratio(2));
    end
    
    wheeltorques_Nm_filt = wheeltorques_Nm;
    for g = 1:size(gbi_radps,1)
        pts = find(wheeltorques_Nm(g,:) < 1e-3);
        if ~isempty(pts)
            wheeltorques_Nm_filt(g, pts) = interp1([v_mph(pts(1)-1), v_mph(pts(1)-1) + 5, v_mph(pts(1)-1) + 100], [wheeltorques_Nm(g, pts(1)-1) 0 0], v_mph(pts));
        end
    end
    
    wheeltorques_Nm = wheeltorques_Nm_filt;
    
    max_etorque_Nm = max(etorques_Nm);
    max_wheeltorque_Nm = max(wheeltorques_Nm);
    
    wheelpowers_kW = wheeltorques_Nm .* repmat(wheelspeed_radps, size(wheeltorques_Nm,1), 1) / 1000;
    max_wheelpower_kW = max_wheeltorque_Nm .* wheelspeed_radps / 1000;
else
    gbi_radps = max(stall_radps, kron(gbo_radps, transmission.gear.ratio(3)'));
    gbi_radps(gbi_radps > redline_radps) = redline_radps;   % special case for CVT, able to hold redline
    
    etorques_Nm = interp1(engine.full_throttle_speed_radps, engine.full_throttle_torque_Nm, gbi_radps, 'linear', 0);

    wheeltorques_Nm = etorques_Nm .* repmat(transmission.gear.ratio(3)', 1, length(gbo_radps)) * drive_axle.final_drive.gear_ratio;
    
    if stall_radps > 0
        wheeltorques_Nm(1,:) = wheeltorques_Nm(1,:) .* interp1([0 stall_radps 1e6], [1.8, 1, 1], gbo_radps .* transmission.gear.ratio(3));
    end
    
    max_etorque_Nm = max(etorques_Nm);
    
    wheelpowers_kW = min(engine.max_power_W/1000, wheeltorques_Nm .* wheelspeed_radps / 1000);
    max_wheelpower_kW = wheelpowers_kW; % max_wheeltorque_Nm .* wheelspeed_radps / 1000;

    wheeltorques_Nm = wheelpowers_kW * 1000 ./ wheelspeed_radps;
    
    max_wheeltorque_Nm = wheeltorques_Nm;
end

if do_plots
    fplothg(engine.full_throttle_speed_radps * convert.radps2rpm, engine.full_throttle_torque_Nm); title('Engine Torque');
    fplothg(engine.full_throttle_speed_radps * convert.radps2rpm, engine.full_throttle_torque_Nm .* engine.full_throttle_speed_radps / 1000); title('Engine Power');
    
    fplothg(v_mph, gbi_radps * convert.radps2rpm); title('input RPM');
    
    fplothg(v_mph, wheeltorques_Nm);
    plothg(v_mph, max_wheeltorque_Nm,'m.-');
    
    fplothg(v_mph, wheelpowers_kW);
    plothg(v_mph, max_wheelpower_kW,'m.-');
    
    fplothg(v_mph, gbo_radps);
    plothg(v_mph, wheelspeed_radps,'m.-');
end

end

