%
% PriusEfficiency.m
%
% Prius HSD efficiency charts in Octave / Matlab. It all started with the question "how to minimize the flow over the electric path"
% (See http://www.priusfreunde.de/portal/index.php?option=com_kunena&Itemid=117&func=view&catid=13&id=236212 [German])
%
% The basic concept is that, assuming battery current is zero, all power going over the electic path is either generated by MG1
% consumed by MG2 or vice versa. The mechanical power of MG1 can be calculated by multiplying RPM and torque. MG1 RPM can be calculated
% from ICE RPM and vehicle speed (MG2 RPM is linear speed). MG1 torque can be calculated from ICE torque: MG1 torque is always ~28% of
% the ICE torque. Assuming that the ICE follows its "basic operating line" we know the ICE torque for each ICE RPM. This is all we
% need to calculate for combination of ICE RPM and vehicle speed the ICE power and MG1 power. Comparing those two values gives us an
% idea how the power is distributed between the mechanical and electric link.
%
% Licence CC BY-NC-SA 4.0, author proprius, contact me at http://www.priusfreunde.de or http://priuschat.com
%
% Version 1.0, 2014-02-02, proprius, Initial version for NHW20
% Version 1.1, 2014-04-18, proprius, MG and inverter efficiency now depend on torque and rpm
%                                    Added transmission losses
%                                    Option to use ICE power for Y-axis. That way low torque situations at minimum RPM can be displayed (<= 9.4 kW).
%                                    Added more plots
% 
% The script makes some assumptions and simplifications. So reality will be a bit different. Here is what I can think of:
%
% - Battery current is assumed to be zero
% - Values like rolling and air resistance, fuel and air density, and weight are set to default values. In reality they depend
%   on many factors, e. g. temperature, tires, etc.
% - Constant speed plot is for flat ground and no wind
% - Torque/BSFC values are taken from chart in Toyota paper SAE 2004-01-0064,
%   "Development of New-Generation Hybrid System THS II - Drastic Improvement of Power Performance and Fuel Economy,"
%   (see http://priuschat.com/threads/bsfc-challenge.55947/ )
%   + It is assumed that the basic operating line of followed, which in reality is not always true, see Fig. 22 in the link above.
%   + ICE is assumed to be warmed up
%   + The original chart only goes up to 4500 RPM. I added values for 5000 RPM based on the max output power of 57 kW
%     for the ICE. The BSCF at 5000 RPM is just a good guess based on the 4500 RPM chart.
% - I removed speeds below 10 km/h, because numbers get weird in that range
%
% Resources (Thanks to all the authors):
%
% http://eahart.com/prius/psd/
% http://prius.ecrostech.com/original/Understanding/PowerSplitDevice.htm
% http://priuschat.com/threads/bsfc-challenge.55947/
% http://home.arcor.de/sg_temp/prius/fahrprofilrechner.htm
% http://priuschat.com/threads/mg1-rpm.59608/#post-818017
% http://priuschat.com/threads/ice-to-wheel-transmission-efficiency-value-at-60-mph.115094/
% http://priuschat.com/threads/introduction-to-prius-power-flow.29352/
% http://priuschat.com/threads/two-mg1-questions.117253/#post-1712003
% http://ecee.colorado.edu/~ecen5017/notes/OakRidge_2004Prius.pdf / http://www.osti.gov/scitech/biblio/890029
% http://info.ornl.gov/sites/publications/files/Pub26762.pdf
% http://www.ae.pwr.wroc.pl/filez/20110606092430_HEV_Toyota.pdf
% http://www.hybrid-autos.info/Technik/E-Maschinen/wirkungsgrad-pmsm.html

% Generic constants
global FUEL_DENSITY AIR_DENSITY NM_RPM_PER_KW;

FUEL_DENSITY                         = 0.750;
AIR_DENSITY                          = 1.225;
NM_RPM_PER_KW                        = 9549.3;

function retval = powerFromRpmAndTorque (rpm,torque)
  global NM_RPM_PER_KW;
  retval = (rpm .* torque) ./ NM_RPM_PER_KW;
end

function retval = torqueFromPowerAndRpm (power,rpm)
  global NM_RPM_PER_KW;
  retval = (power ./ rpm) .* NM_RPM_PER_KW;
end

% Vehicle specific constants
global MODEL INVERTER_AVG_EFFICIENCY MG_EFFICIENCY_RPM MG_EFFICIENCY MG2_AVERAGE_EFFICIENTY
global MECHANICAL_EFFICIENCY REDUCTION_GEAR_LOSS SELF_CONSUMPTION ICE_START_ENERGY_LITER;
global SUN_GEAR_TEETH RING_GEAR_TEETH SUN_TO_RING SUN_TO_CARRIER MG1_ICE_TORQUE_FACTOR;
global MAX_TORQUE_MG1 MAX_RPM_MG1 MIN_RPM_MG1 MAX_TORQUE_MG2 MAX_RPM_MG2 MAX_KMH_FOR_MG2;
global MAX_KMH MAX_KMH_SPEEDOMETER MIN_KMH_SPEEDOMETER SPEEDOMETER_FACTOR MIN_RPM_ICE MAX_RPM_ICE;
global ICE_RPM_TABLE ICE_TORQUE_TABLE ICE_BSFC_TABLE;
global ICE_RPM_TABLE_INCLUDING_LOW_TORQUE ICE_TORQUE_TABLE_INCLUDING_LOW_TORQUE;
global ICE_BSFC_TABLE_INCLUDING_LOW_TORQUE ICE_POWER_TABLE_INCLUDING_LOW_TORQUE;
global MIN_WHEELPOWER_ICE_RUNNING;
global CW FRONTAREA ROLLING WEIGHT;
global ICE_DRAG_FUEL_CUT_OFF MAX_KMH_ICE_NOT_TURNING;
global ICE_TORQUE_TABLE_FULL ICE_BSFC_TABLE_FULL;
global ICE_RPM_TABLE_INCLUDING_LOW_TORQUE_FULL;
global ICE_TORQUE_TABLE_INCLUDING_LOW_TORQUE_FULL;
global ICE_BSFC_TABLE_INCLUDING_LOW_TORQUE_FULL;
global INTERPOLATION_FACTOR;

MODEL                                = 'Toyota Prius 2004-2009 (NHW20)';
INVERTER_AVG_EFFICIENCY              = 0.98;
MG2_AVERAGE_EFFICIENTY               = 0.9;
MECHANICAL_EFFICIENCY                = 0.97;
REDUCTION_GEAR_LOSS                  = 1.4;                                                         % kW loss at max speed
SELF_CONSUMPTION                     = 0.5;
ICE_START_ENERGY_LITER               = 0.001 .* 2 .* 240 .* (FUEL_DENSITY ./ 1000);                 % Assume 0.001 kWh electric energy to start ICE, energy storage efficiency of 50% and BSFC of 240 g/kWh for charging
ICE_DRAG_FUEL_CUT_OFF                = 2;                                                           % kW to turn ICE over without fuel
SUN_GEAR_TEETH                       = 30;
RING_GEAR_TEETH                      = 78;
SUN_TO_RING                          = RING_GEAR_TEETH ./ SUN_GEAR_TEETH;                           % 2.6
SUN_TO_CARRIER                       = SUN_TO_RING .+ 1;                                            % 3.6
MG1_ICE_TORQUE_FACTOR                = -(SUN_GEAR_TEETH ./ (SUN_GEAR_TEETH .+ RING_GEAR_TEETH));    % ~ -0.2777
MAX_TORQUE_MG1                       = 45;
MAX_RPM_MG1                          = 10000;
MIN_RPM_MG1                          = -MAX_RPM_MG1;
MAX_TORQUE_MG2                       = 400;
MAX_RPM_MG2                          = 6400;
MAX_KMH_FOR_MG2                      = 181.9;
SPEEDOMETER_FACTOR                   = 1.075;                                                       % speedometer = true speed + 7,5%
MAX_KMH                              = 170;
MAX_KMH_ICE_NOT_TURNING              = 65;
MAX_KMH_SPEEDOMETER                  = MAX_KMH .* SPEEDOMETER_FACTOR;
MIN_KMH_SPEEDOMETER                  = 10;
MIN_RPM_ICE                          = 1200;
MAX_RPM_ICE                          = 5000;
ICE_RPM_TABLE                        = [1200 1290 1700 2000 2220   2800 3230 3500 4500 5000];
ICE_TORQUE_TABLE                     = [  75   77   81   85   86.5   91   94   97   97  109];
ICE_BSFC_TABLE                       = [ 242  240  235  232  230    230  230  230  235  240];
ICE_RPM_TABLE_INCLUDING_LOW_TORQUE   = [1200 1200 1200 1200 1200 1200 1200 1200 1200 ICE_RPM_TABLE];
ICE_TORQUE_TABLE_INCLUDING_LOW_TORQUE = [  14   18   24   35   43   48   55   63   68 ICE_TORQUE_TABLE];
ICE_BSFC_TABLE_INCLUDING_LOW_TORQUE  = [ 500  400  350  300  280  270  260  250  245 ICE_BSFC_TABLE];
ICE_POWER_TABLE_INCLUDING_LOW_TORQUE = powerFromRpmAndTorque(ICE_RPM_TABLE_INCLUDING_LOW_TORQUE, ICE_TORQUE_TABLE_INCLUDING_LOW_TORQUE);
MIN_WHEELPOWER_ICE_RUNNING           = 10;
CW                                   = 0.26;
FRONTAREA                            = 2.23;
ROLLING                              = 0.01;
WEIGHT                               = 1375;
INTERPOLATION_FACTOR                 = 2000;

function retval = reductionGearLoss(kmh)
  global REDUCTION_GEAR_LOSS MAX_KMH;
  retval = (kmh .* REDUCTION_GEAR_LOSS) ./ MAX_KMH; 
end  

function retval = powerPsdToPowerWheel (kmh,powerPsd)
  global MECHANICAL_EFFICIENCY SELF_CONSUMPTION;
  retval = (powerPsd .* MECHANICAL_EFFICIENCY) .- (reductionGearLoss(kmh) .+ SELF_CONSUMPTION);
end

function retval = powerRequiredFromPsdForGivenWheelPower (kmh,powerWheel)
  global MECHANICAL_EFFICIENCY;
  retval = (powerWheel .+ reductionGearLoss(kmh)) ./ MECHANICAL_EFFICIENCY;
end

function retval = efficiencyFactorPmsmRpm (normalizedRpm)
  v = normalizedRpm - 0.2;
  retval = ((6.5 .+ log( normalizedRpm ./ 5 .+ 0.01 )) .* (2 .- v .* v .* 1.4)) .* 0.123;
end 

function retval = efficiencyFactorPmsmTorque (normalizedTorque)
  retval = ((6.5 .+ log( abs(normalizedTorque) ./ 5 .+ 0.01 )) .* (5 .- abs(normalizedTorque))) .* 0.05;
end 

function retval = efficiencyFactorPmsmHighTorqueHighRpm (normalizedRpm,normalizedTorque)
  w = abs(0.75 * normalizedTorque);
  z = (normalizedRpm .* (1.0 + normalizedRpm ./ 5) + 0.11) .* abs(w);
  retval = ((5.0 .+ log( z .* 3 .+ 0.06 )) .* (0.50 .- z.*z.*10)) .* 0.75;
  retval = retval .- normalizedRpm .* 0.2;
end 

function retval = efficiencyFactorPmsmHighTorqueLowRpm (normalizedRpm,normalizedTorque)
  w = 0.8 .* (abs(normalizedTorque) .- 1 );
  retval = 2 .* (normalizedRpm .* normalizedRpm .* 10 .+ w .* w + w .* 0.3) .+ 0.5;
  retval = retval .* (retval < 1) .+ (retval >= 1);
  retval = retval .* (retval > 0);
end 

function retval = efficiencyFactorPmsmGenerator (normalizedRpm,normalizedTorque)
  retval = (normalizedTorque .+ 10) ./ 10;
  retval = retval .* ((normalizedRpm .* 20 .+ normalizedTorque) > 0);
  retval = retval .* (normalizedTorque < 0) .+ (normalizedTorque >= 0);
end 

function retval = efficiencyPmsmNormalized (normalizedRpm,normalizedTorque)
  retval = efficiencyFactorPmsmRpm(normalizedRpm);
  retval = retval .* efficiencyFactorPmsmTorque(normalizedTorque);
  retval = retval .* efficiencyFactorPmsmHighTorqueHighRpm(normalizedRpm,normalizedTorque);
  retval = retval .* efficiencyFactorPmsmHighTorqueLowRpm(normalizedRpm,normalizedTorque);
  retval = retval .* efficiencyFactorPmsmGenerator(normalizedRpm,normalizedTorque);
  retval = (0.6 .+ 0.365 .* retval);
  retval = retval .* (retval > 0.6) .+ (retval <= 0.6) .* 0.6;
end 

function retval = efficiencyMg2 (rpm,nm)
  global MAX_RPM_MG2 MAX_TORQUE_MG2;
  _nm = nm .* ((rpm > 0) - (rpm <= 0));
  retval = efficiencyPmsmNormalized(abs(rpm) ./ MAX_RPM_MG2,_nm ./ MAX_TORQUE_MG2);
end

function retval = efficiencyMg1 (rpm,nm)
  global MAX_RPM_MG1 MAX_TORQUE_MG1;
  _nm = nm .* ((rpm > 0) - (rpm <= 0));
  retval = efficiencyPmsmNormalized(0.8 .* abs(rpm) ./ MAX_RPM_MG1,0.4 .* _nm ./ MAX_TORQUE_MG1);
end

function retval = efficiencyInverterNormalized (x,y)
  v = x .+ abs(y) ./ 40;
  retval = 0.96 .+ (log(v .* 10 .+ 0.02) ./ 50) .* ((3 .- v) ./ 3);
end

function retval = efficiencyInverter (rpm,nm)
  % Use max rpm and torque from MG2 also for MG1 as I assume that efficiency depends on absolute rpm and torque
  global MAX_RPM_MG2 MAX_TORQUE_MG2;
  _nm = nm .* ((rpm > 0) - (rpm <= 0));
  retval = efficiencyInverterNormalized(abs(rpm) ./ MAX_RPM_MG2,_nm ./ MAX_TORQUE_MG2);
end

function retval = gramToLiter (gram)
  global FUEL_DENSITY;
  retval = gram ./ (FUEL_DENSITY .* 1000);
end  

function retval = kmhToMeterPerSec (kmh)
  retval = kmh ./ 3.6;
end

function retval = speedometerToTrueSpeed (kmh)
  global SPEEDOMETER_FACTOR;
  retval = kmh ./ SPEEDOMETER_FACTOR;
end

function retval = powerConstantSpeedOnFlatGround (kmh)
  % Derived from (c) Hjermind 2008-2009 http://home.arcor.de/sg_temp/prius/fahrprofilrechner.htm
  global CW FRONTAREA ROLLING WEIGHT AIR_DENSITY;
  v = kmhToMeterPerSec(kmh);
  air = v .* (AIR_DENSITY ./ 2) .* CW .* FRONTAREA .* v .* v ./ 1000;
  roll = v .* ROLLING .* WEIGHT .* 9.81 ./ 1000;
  retval = air .+ roll; 
end

function retval = kmhToRpmMg2 (kmh)
  global MAX_RPM_MG2 MAX_KMH_FOR_MG2; 
  retval = (kmh .* MAX_RPM_MG2) ./ MAX_KMH_FOR_MG2;
end  

function retval = rpmMg2ToKmh (rpmMg2)
  global MAX_RPM_MG2 MAX_KMH_FOR_MG2;
  retval = MAX_KMH_FOR_MG2 .* (rpmMg2 ./ MAX_RPM_MG2);
end  

function retval = rpmMg1 (rpmMg2,rpmIce)
  global SUN_TO_RING SUN_TO_CARRIER;
  retval = (SUN_TO_CARRIER .* rpmIce) .- (SUN_TO_RING .* rpmMg2);
end

function retval = rpmMg2 (rpmMg1,rpmIce)
  global SUN_TO_RING SUN_TO_CARRIER;
  retval = ((SUN_TO_CARRIER .* rpmIce) .- rpmMg1) ./ SUN_TO_RING;
end

function retval = rpmIce (rpmMg1,rpmMg2)
  global SUN_TO_RING SUN_TO_CARRIER;
  retval = (rpmMg1 .+ (SUN_TO_RING .* rpmMg2)) ./ SUN_TO_CARRIER;
end

function retval = rpmMg1IceToKmh (rpmMg1,rpmIce)
  retval = rpmMg2ToKmh(rpmMg2(rpmMg1,rpmIce));
end

function retval = powerMg1relativeToIce (rpmMg2,rpmIce)
  global MG1_ICE_TORQUE_FACTOR;
  retval = (rpmMg1(rpmMg2,rpmIce) ./ rpmIce) .* MG1_ICE_TORQUE_FACTOR;
end

function retval = powerMg1relativeToIceKmh (kmh, rpmIce)
  retval = powerMg1relativeToIce(kmhToRpmMg2(kmh),rpmIce);
end  

function retval = rpmToTorqueIceInterp1 (rpmIce)
  global ICE_RPM_TABLE ICE_TORQUE_TABLE;
  retval = interp1(ICE_RPM_TABLE,ICE_TORQUE_TABLE,rpmIce) .* (rpmIce > 0);
end

function retval = rpmToBsfcIceInterp1 (rpmIce)
  global ICE_RPM_TABLE ICE_BSFC_TABLE;
  retval = interp1(ICE_RPM_TABLE,ICE_BSFC_TABLE,rpmIce) .* (rpmIce > 0);
end

function retval = powerToRpmIceInterp1 (powerIce)
  global ICE_RPM_TABLE_INCLUDING_LOW_TORQUE ICE_POWER_TABLE_INCLUDING_LOW_TORQUE;
  retval = interp1(ICE_POWER_TABLE_INCLUDING_LOW_TORQUE,ICE_RPM_TABLE_INCLUDING_LOW_TORQUE,powerIce,'extrap');
end

function retval = powerToTorqueIceInterp1 (powerIce)
  global ICE_TORQUE_TABLE_INCLUDING_LOW_TORQUE ICE_POWER_TABLE_INCLUDING_LOW_TORQUE;
  retval = interp1(ICE_POWER_TABLE_INCLUDING_LOW_TORQUE,ICE_TORQUE_TABLE_INCLUDING_LOW_TORQUE,powerIce,'extrap');
end

function retval = powerToBsfcIceInterp1 (powerIce)
  global ICE_BSFC_TABLE_INCLUDING_LOW_TORQUE ICE_POWER_TABLE_INCLUDING_LOW_TORQUE;
  retval = interp1(ICE_POWER_TABLE_INCLUDING_LOW_TORQUE,ICE_BSFC_TABLE_INCLUDING_LOW_TORQUE,powerIce,'extrap');
end

function initialize ()
  global ICE_TORQUE_TABLE_FULL ICE_BSFC_TABLE_FULL;
  global ICE_RPM_TABLE ICE_TORQUE_TABLE ICE_POWER_TABLE_INCLUDING_LOW_TORQUE;
  global ICE_RPM_TABLE_INCLUDING_LOW_TORQUE_FULL;
  global ICE_TORQUE_TABLE_INCLUDING_LOW_TORQUE_FULL;
  global ICE_BSFC_TABLE_INCLUDING_LOW_TORQUE_FULL;
  global INTERPOLATION_FACTOR;
  ICE_RPM_TABLE_FULL = [0:1/INTERPOLATION_FACTOR:max(ICE_RPM_TABLE)];
  ICE_TORQUE_TABLE_FULL = rpmToTorqueIceInterp1(ICE_RPM_TABLE_FULL);
  ICE_BSFC_TABLE_FULL = rpmToBsfcIceInterp1(ICE_RPM_TABLE_FULL);
  ICE_POWER_TABLE_FULL = [0:1/INTERPOLATION_FACTOR:max(ICE_POWER_TABLE_INCLUDING_LOW_TORQUE)+3];
  ICE_RPM_TABLE_INCLUDING_LOW_TORQUE_FULL = powerToRpmIceInterp1(ICE_POWER_TABLE_FULL);
  ICE_TORQUE_TABLE_INCLUDING_LOW_TORQUE_FULL = powerToTorqueIceInterp1(ICE_POWER_TABLE_FULL);
  ICE_BSFC_TABLE_INCLUDING_LOW_TORQUE_FULL = powerToBsfcIceInterp1(ICE_POWER_TABLE_FULL);
end

function retval = rpmToTorqueIce (rpmIce)
  global ICE_TORQUE_TABLE_FULL INTERPOLATION_FACTOR;
  retval = ICE_TORQUE_TABLE_FULL(round(rpmIce .* INTERPOLATION_FACTOR) .+ 1);
end

function retval = rpmToBsfcIce (rpmIce)
  global ICE_BSFC_TABLE_FULL INTERPOLATION_FACTOR;
  retval = ICE_BSFC_TABLE_FULL(round(rpmIce .* INTERPOLATION_FACTOR) .+ 1);
end

function retval = powerToRpmIce (powerIce)
  global ICE_RPM_TABLE_INCLUDING_LOW_TORQUE_FULL INTERPOLATION_FACTOR;
  retval = ICE_RPM_TABLE_INCLUDING_LOW_TORQUE_FULL(round(powerIce .* INTERPOLATION_FACTOR) .+ 1);
end

function retval = powerToTorqueIce (powerIce)
  global ICE_TORQUE_TABLE_INCLUDING_LOW_TORQUE_FULL INTERPOLATION_FACTOR;
  retval = ICE_TORQUE_TABLE_INCLUDING_LOW_TORQUE_FULL(round(powerIce .* INTERPOLATION_FACTOR) .+ 1);
end

function retval = powerToBsfcIce (powerIce)
  global ICE_BSFC_TABLE_INCLUDING_LOW_TORQUE_FULL INTERPOLATION_FACTOR;
  retval = ICE_BSFC_TABLE_INCLUDING_LOW_TORQUE_FULL(round(powerIce .* INTERPOLATION_FACTOR) .+ 1);
end

%function retval = rpmToTorqueIce (rpmIce)
%  global ICE_RPM_TABLE ICE_TORQUE_TABLE;
%  retval = interp1(ICE_RPM_TABLE,ICE_TORQUE_TABLE,rpmIce) .* (rpmIce > 0);
%end
%
%function retval = rpmToBsfcIce (rpmIce)
%  global ICE_RPM_TABLE ICE_BSFC_TABLE;
%  retval = interp1(ICE_RPM_TABLE,ICE_BSFC_TABLE,rpmIce) .* (rpmIce > 0);
%end
%
%function retval = powerToRpmIce (powerIce)
%  global ICE_RPM_TABLE_INCLUDING_LOW_TORQUE ICE_POWER_TABLE_INCLUDING_LOW_TORQUE;
%  retval = interp1(ICE_POWER_TABLE_INCLUDING_LOW_TORQUE,ICE_RPM_TABLE_INCLUDING_LOW_TORQUE,powerIce);
%end
%
%function retval = powerToTorqueIce (powerIce)
%  global ICE_TORQUE_TABLE_INCLUDING_LOW_TORQUE ICE_POWER_TABLE_INCLUDING_LOW_TORQUE;
%  retval = interp1(ICE_POWER_TABLE_INCLUDING_LOW_TORQUE,ICE_TORQUE_TABLE_INCLUDING_LOW_TORQUE,powerIce);
%end
%
%function retval = powerToBsfcIce (powerIce)
%  global ICE_BSFC_TABLE_INCLUDING_LOW_TORQUE ICE_POWER_TABLE_INCLUDING_LOW_TORQUE;
%  retval = interp1(ICE_POWER_TABLE_INCLUDING_LOW_TORQUE,ICE_BSFC_TABLE_INCLUDING_LOW_TORQUE,powerIce);
%end

function retval = powerBattery (kmh,powerWheel)
  global MIN_WHEELPOWER_ICE_RUNNING;
  _rpmMg2 = kmhToRpmMg2(kmh);
  _powerMg2 = powerRequiredFromPsdForGivenWheelPower(kmh,powerWheel);
  _torqueMg2 = torqueFromPowerAndRpm(_powerMg2,_rpmMg2);
  _efficiencyMg2 = efficiencyMg2(_rpmMg2,_torqueMg2);
  _efficiencyInverterMg2 = efficiencyInverter(_rpmMg2,_torqueMg2);
  retval = (_powerMg2 ./ (_efficiencyMg2 .* _efficiencyInverterMg2)) .* (powerWheel < MIN_WHEELPOWER_ICE_RUNNING);
end

function retval = rpmToPowerIce (rpmIce)
  retval = powerFromRpmAndTorque(rpmIce,rpmToTorqueIce(rpmIce));
end

function retval = torqueMg1 (rpmIce)
  global MG1_ICE_TORQUE_FACTOR;
  retval = rpmToTorqueIce(rpmIce) .* MG1_ICE_TORQUE_FACTOR;
end

function retval = torqueMg1ForPowerIce (powerIce)
  global MG1_ICE_TORQUE_FACTOR;
  retval = powerToTorqueIce(powerIce) .* MG1_ICE_TORQUE_FACTOR;
end

function retval = powerMg1 (kmh,rpmIce)
  if min(rpmIce) == 0;
    retval = 0
  else
    retval = rpmToPowerIce(rpmIce) .* powerMg1relativeToIceKmh(kmh,rpmIce);
  end
end  

function retval = powerMg1ForKmhAndPowerIce (kmh,powerIce)
  retval = powerIce .* powerMg1relativeToIceKmh(kmh,powerToRpmIce(powerIce));
end  
 
function retval = electricPathEfficiencyForMg1Mg2 (_rpmMg1,_torqueMg1,_powerMg1,_rpmMg2)
  global INVERTER_AVG_EFFICIENCY MG2_AVERAGE_EFFICIENTY;
  _efficiencyMg1 = efficiencyMg1(_rpmMg1,_torqueMg1);
  _efficiencyInverterMg1 = efficiencyInverter(_rpmMg1,_torqueMg1);
  
  % In order to get the exact efficiency for MG2, we first have to approximate the torque by
  % using an average efficiency for MG2 and MG2 inverter. We also have to distinguish between
  % power flow from MG1 to MG2 and vice versa.
  _mg1IsGenerator = (_torqueMg1 .* _rpmMg1 <= 0);
  _mg1IsMotor = 1 .- _mg1IsGenerator;
  _approximateEfficiencyElectricPath = _efficiencyMg1 .* MG2_AVERAGE_EFFICIENTY .* _efficiencyInverterMg1 .* INVERTER_AVG_EFFICIENCY;
  _approximatePowerMg2WhenMg1IsGenerator = -_powerMg1 .* _approximateEfficiencyElectricPath;
  _approximatePowerMg2WhenMg1IsMotor = -_powerMg1./ _approximateEfficiencyElectricPath;
  _approximatePowerMg2 = _mg1IsGenerator .* _approximatePowerMg2WhenMg1IsGenerator .+ _mg1IsMotor .* _approximatePowerMg2WhenMg1IsMotor;
  _approximateTorqueMg2 = torqueFromPowerAndRpm(_approximatePowerMg2,_rpmMg2);
  
  _efficiencyMg2 = efficiencyMg2(_rpmMg2,_approximateTorqueMg2);
  _efficiencyInverterMg2 = efficiencyInverter(_rpmMg2,_approximateTorqueMg2);
  retval = _efficiencyMg1 .* _efficiencyMg2 .* _efficiencyInverterMg1 .* _efficiencyInverterMg2;
end

function retval = electricPathEfficiency (kmh,rpmIce)
  _rpmMg2 = kmhToRpmMg2(kmh);
  _rpmMg1 = rpmMg1(_rpmMg2,rpmIce);
  _torqueMg1 = torqueMg1(rpmIce);
  _powerMg1 = powerMg1(kmh,rpmIce);
  retval = electricPathEfficiencyForMg1Mg2(_rpmMg1,_torqueMg1,_powerMg1,_rpmMg2) .* (rpmIce > 0) .+ 100 .* (rpmIce == 0);
end

function retval = electricPathEfficiencyForKmhAndPowerIce (kmh,powerIce)
  _rpmMg2 = kmhToRpmMg2(kmh);
  _rpmMg1 = rpmMg1(_rpmMg2,powerToRpmIce(powerIce));
  _torqueMg1 = torqueMg1ForPowerIce(powerIce);
  _powerMg1 = powerMg1ForKmhAndPowerIce(kmh,powerIce);
  retval = electricPathEfficiencyForMg1Mg2(_rpmMg1,_torqueMg1,_powerMg1,_rpmMg2);
end

function [generatorPower motorPower] = generatorAndMotorPowerForPowerMg1AndEfficiency (_powerMg1,_electricPathEfficiency)
  _mg1IsGenerator = (_powerMg1 <= 0);
  _mg1IsMotor = 1 .- _mg1IsGenerator;
  generatorPower = _mg1IsGenerator .* -_powerMg1                            .+ _mg1IsMotor .* (_powerMg1 ./ _electricPathEfficiency);
  motorPower     = _mg1IsGenerator .* -_powerMg1 .* _electricPathEfficiency .+ _mg1IsMotor .* _powerMg1;
end

function [generatorPower motorPower] = generatorAndMotorPower (kmh,rpmIce)
  _powerMg1 = powerMg1(kmh,rpmIce);
  _electricPathEfficiency = electricPathEfficiency(kmh,rpmIce);
  [generatorPower motorPower] = generatorAndMotorPowerForPowerMg1AndEfficiency(_powerMg1,_electricPathEfficiency);
end

function [generatorPower motorPower] = generatorAndMotorPowerForKmhAndPowerIce (kmh,powerIce)
  _powerMg1 = powerMg1ForKmhAndPowerIce(kmh,powerIce);
  _electricPathEfficiency = electricPathEfficiencyForKmhAndPowerIce(kmh,powerIce);
  [generatorPower motorPower] = generatorAndMotorPowerForPowerMg1AndEfficiency(_powerMg1,_electricPathEfficiency);
end

function retval = electricPathLossKw (kmh,rpmIce)
  [generatorPower motorPower] = generatorAndMotorPower(kmh,rpmIce);
  retval = generatorPower .- motorPower;
end

function retval = electricPathLossKwForKmhAndPowerIce (kmh,powerIce)
  [generatorPower motorPower] = generatorAndMotorPowerForKmhAndPowerIce(kmh,powerIce);
  retval = generatorPower .- motorPower;
end

function retval = generatorPower (kmh,rpmIce)
  [generatorPower motorPower] = generatorAndMotorPower(kmh,rpmIce);
  retval = generatorPower;
end

function retval = powerMg2 (kmh,rpmIce)
  _powerMg1 = powerMg1(kmh,rpmIce);
  _mg1IsGenerator = (_powerMg1 <= 0);
  _mg1IsMotor = 1 .- _mg1IsGenerator;
  _electricPathEfficiency = electricPathEfficiency(kmh,rpmIce);
  _powerMg2WhenMg1IsGenerator = -_powerMg1 .* _electricPathEfficiency;
  _powerMg2WhenMg1IsMotor     = -_powerMg1 ./ _electricPathEfficiency;
  retval = _mg1IsGenerator .* _powerMg2WhenMg1IsGenerator .+ _mg1IsMotor .* _powerMg2WhenMg1IsMotor;
end

function retval = torqueMg2 (kmh,rpmIce)
  _powerMg2 = powerMg2(kmh,rpmIce);
  _rpmMg2 = kmhToRpmMg2(kmh);
  retval = torqueFromPowerAndRpm(_powerMg2,_rpmMg2);
end

function retval = powerAfterPsd (kmh,rpmIce)
  retval = rpmToPowerIce(rpmIce) .- electricPathLossKw(kmh,rpmIce);
end 

function retval = powerAfterPsdForKmhAndPowerIce (kmh,powerIce)
  retval = powerIce .- electricPathLossKwForKmhAndPowerIce(kmh,powerIce);
end 

function retval = powerWheelForKmhAndRpmIce (kmh,rpmIce)
  retval = powerPsdToPowerWheel(kmh,powerAfterPsd(kmh,rpmIce));
end

function retval = powerWheelForKmhAndPowerIce (kmh,powerIce)
  retval = powerPsdToPowerWheel(kmh,powerAfterPsdForKmhAndPowerIce(kmh,powerIce));
end

function retval = powerSurplus (kmh,rpmIce)
  retval = powerWheelForKmhAndRpmIce(kmh,rpmIce) .- powerConstantSpeedOnFlatGround(kmh);
end 

function retval = powerSurplusForKmhAndPowerIce (kmh,powerIce)
  retval = powerWheelForKmhAndPowerIce(kmh,powerIce) .- powerConstantSpeedOnFlatGround(kmh);
end 

function literPerHundredKm = consumption (kmh,rpmIce)
  literPerHour = gramToLiter(rpmToBsfcIce(rpmIce) .* rpmToPowerIce(rpmIce));
  literPerHundredKm = (literPerHour .* 100) ./ kmh;
end

function literPerHundredKm = consumptionByPowerIce (kmh,powerIce)
  literPerHour = gramToLiter(powerToBsfcIce(powerIce) .* powerIce);
  literPerHundredKm = (literPerHour .* 100) ./ kmh;
end

function retval = rpmMg1kmh (kmh, rpmIce)
  retval = rpmMg1(kmhToRpmMg2(kmh),rpmIce);
end  

function retval = minKmhSpeedometer (rpmIce)
  global MAX_RPM_MG1 MIN_KMH_SPEEDOMETER SPEEDOMETER_FACTOR;
  retval = rpmMg1IceToKmh(MAX_RPM_MG1,rpmIce) .* SPEEDOMETER_FACTOR;
  retval = retval .* (retval >= MIN_KMH_SPEEDOMETER) + MIN_KMH_SPEEDOMETER .* (retval < MIN_KMH_SPEEDOMETER);
end  

function retval = maxKmhSpeedometer (rpmIce)
  global MIN_RPM_MG1 MAX_KMH_SPEEDOMETER SPEEDOMETER_FACTOR;
  retval = rpmMg1IceToKmh(MIN_RPM_MG1,rpmIce) .* SPEEDOMETER_FACTOR;
  retval = retval .* (retval <= MAX_KMH_SPEEDOMETER) + MAX_KMH_SPEEDOMETER .* (retval > MAX_KMH_SPEEDOMETER);
end  

function retval = efficiency (kmh,rpmIce)
  global ICE_BSFC_TABLE;
  bsfc = rpmToBsfcIce(rpmIce);
  retval = 100 .* (powerWheelForKmhAndRpmIce(kmh,rpmIce) ./ rpmToPowerIce(rpmIce)) ./ (bsfc ./ min(ICE_BSFC_TABLE));
end  

function retval = efficiencyForKmhAndPowerIce (kmh,powerIce)
  global ICE_BSFC_TABLE;
  bsfc = powerToBsfcIce(powerIce);
  retval = 100 .* (powerWheelForKmhAndPowerIce(kmh,powerIce) ./ powerIce) ./ (bsfc ./ min(ICE_BSFC_TABLE));
end  

function retval = accelerationForPowerWheel (kmh,powerWheel)
  global WEIGHT;
  % a = P / (m * v)
  speedMeterPerSecond = kmh ./ 3.6;
  accelerationMeterPerSquareSecond = (powerWheel .* 1000) ./ (WEIGHT .* speedMeterPerSecond);
  retval = accelerationMeterPerSquareSecond .* 3.6;
end

function retval = accelerationForRpmIce (kmh,rpmIce)
  retval = accelerationForPowerWheel(kmh,powerSurplus(kmh,rpmIce));
end

function retval = accelerationForPowerIce (kmh,powerIce)
  if powerIce > 0
    powerWheel = powerSurplusForKmhAndPowerIce(kmh,powerIce);
  else
    powerWheel = powerWheelEngineOff(kmh); 
  end
  retval = accelerationForPowerWheel(kmh,powerWheel);
end

function [acceleraton rpmIce] = accelerationForRpmMg1 (kmh,rpmMg1)
  global MIN_RPM_ICE;
  rpmIce = rpmIce(rpmMg1,kmhToRpmMg2(kmh));
  rpmIce = max(rpmIce,MIN_RPM_ICE);
  acceleraton = accelerationForRpmIce(kmh,rpmIce);
end

function retval = powerWheelEngineOff (kmh)
  global ICE_DRAG_FUEL_CUT_OFF MAX_KMH_ICE_NOT_TURNING;
  % Result will be negative, as power is taken from the wheels
  retval = powerPsdToPowerWheel(kmh,0) .- powerConstantSpeedOnFlatGround(kmh);
  retval = retval .- (kmh > MAX_KMH_ICE_NOT_TURNING) .* ICE_DRAG_FUEL_CUT_OFF;
end

function retval = accelerationEngineOff (kmh)
  retval = accelerationForPowerWheel(kmh,powerWheelEngineOff(kmh));
end  

function [secondsPulseDuration distancePulse litersPulse kmh rpmIceMin rpmIceMax] = pulse (kmh,kmhHigh,rpmMg1,maxStepSize,minStepSize)
  global MIN_RPM_ICE MAX_RPM_ICE;
  secondsPulseDuration = 0;
  distancePulse = 0;
  rpmIceMin = MAX_RPM_ICE;
  rpmIceMax = MIN_RPM_ICE;
  litersPulse = 0;
  stepSize = maxStepSize;
  while kmh <= kmhHigh
    [acceleration rpmIce] = accelerationForRpmMg1(kmh,rpmMg1);
    rpmIceMin = min(rpmIceMin,rpmIce);
    rpmIceMax = max(rpmIceMax,rpmIce);
    if kmh + acceleration * stepSize > kmhHigh
      stepSize = max([(kmhHigh - kmh) / acceleration, minStepSize]);
    end
    litersPerHour = gramToLiter(rpmToBsfcIce(rpmIce) .* rpmToPowerIce(rpmIce));
    litersPulse = litersPulse + (litersPerHour/3600) * stepSize;
    kmh = kmh + acceleration * stepSize;
    distancePulse = distancePulse + (kmh/3.6) * stepSize;
    secondsPulseDuration = secondsPulseDuration + stepSize;
  end
end

function [duration distance liters kmh rpmIceMin rpmIceMax] = pulseOrGlidePowerIce (kmh,kmhLow,kmhHigh,powerIce,maxStepSize,minStepSize)
  duration = 0;
  distance = 0;
  isPulse = accelerationForPowerIce(kmh,powerIce) > 0;
  stepSize = maxStepSize;
  while (isPulse && kmh <= kmhHigh) || (1-isPulse && kmh >= kmhLow)
    acceleration = accelerationForPowerIce(kmh,powerIce);
    if isPulse && (kmh + acceleration * stepSize > kmhHigh)
      stepSize = max([(kmhHigh - kmh) / acceleration, minStepSize]);
    end
    if (1-isPulse) && (kmh + acceleration * stepSize < kmhLow)
      stepSize = max([(kmhLow - kmh) / acceleration, minStepSize]);
    end
    kmh = kmh + acceleration * stepSize;
    distance = distance + (kmh/3.6) * stepSize;
    duration = duration + stepSize;
  end
  if powerIce > 0
    litersPerHour = gramToLiter(powerToBsfcIce(powerIce) * powerIce);
    liters = duration * (litersPerHour/3600);
    rpmIceMin = rpmIceMax = powerToRpmIce(powerIce);
  else
    litersPerHour = 0;
    liters = 0;
    rpmIceMin = rpmIceMax = 0;
  end
end

function [secondsGlideDuration distanceGlide litersGlide kmh rpmIceMinGlide rpmIceMaxGlide] = glide (kmh,kmhLow,maxStepSize,minStepSize)
  [secondsGlideDuration distanceGlide litersGlide kmh rpmIceMinGlide rpmIceMaxGlide] = pulseOrGlidePowerIce(kmh,kmhLow,kmh,0,maxStepSize,minStepSize);
end

function [consumption rpmIce powerIce] = consumptionForKmh (kmh)
  global ICE_POWER_TABLE_INCLUDING_LOW_TORQUE;
  powerIce = fzero(@(powerIce) powerSurplusForKmhAndPowerIce(kmh,powerIce),[min(ICE_POWER_TABLE_INCLUDING_LOW_TORQUE),max(ICE_POWER_TABLE_INCLUDING_LOW_TORQUE)]);
  consumption = consumptionByPowerIce(kmh,powerIce); 
  rpmIce = powerToRpmIce(powerIce);
end

function [savings consumptionPulseAndGlide] = pulseAndGlide (kmhMiddle,kmhDelta,rpmMg1)
  global SPEEDOMETER_FACTOR;
  global ICE_START_ENERGY_LITER;
  %fprintf('kmhMiddle = %.2f kmhDelta = %.2f rpmMg1 = %.2f\n',kmhMiddle,kmhDelta,rpmMg1);
  maxStepSize = 0.01;
  minStepSize = 0.00001;
  kmhHigh = speedometerToTrueSpeed( kmhMiddle + kmhDelta/2 );
  kmhLow = speedometerToTrueSpeed( kmhMiddle - kmhDelta/2 );
  [secondsPulseDuration distancePulse litersPulse kmh rpmIceMinPulse rpmIceMaxPulse] = pulse(kmhLow,kmhHigh,rpmMg1,maxStepSize,minStepSize);
  %[secondsPulseDuration distancePulse litersPulse kmh rpmIceMinPulse rpmIceMaxPulse] = pulseOrGlidePowerIce(kmhLow,kmhLow,kmhHigh,6,maxStepSize,minStepSize);
  [secondsGlideDuration distanceGlide litersGlide kmh rpmIceMinGlide rpmIceMaxGlide] = glide(kmh,kmhLow,maxStepSize,minStepSize);
  %[secondsGlideDuration distanceGlide litersGlide kmh rpmIceMinGlide rpmIceMaxGlide] = pulseOrGlidePowerIce(kmh,kmhLow,kmhHigh,2,maxStepSize,minStepSize);
  %secondsPulseDuration
  %distancePulse
  %rpmIceMinPulse
  %rpmIceMaxPulse
  %secondsGlideDuration
  %distanceGlide
  %rpmIceMinGlide
  %rpmIceMaxGlide
  liters = ICE_START_ENERGY_LITER + litersPulse + litersGlide;
  kmhAverageTrueSpeed = 3.6 * (distancePulse + distanceGlide) / (secondsPulseDuration + secondsGlideDuration);
  %kmhAverage = SPEEDOMETER_FACTOR .* kmhAverageTrueSpeed;
  consumptionPulseAndGlide = liters * (100000 / (distancePulse + distanceGlide));
  [consumptionConstantSpeed rpmIceConstantSpeed powerIceConstantSpeed] = consumptionForKmh(kmhAverageTrueSpeed);
  savings = 100 * (consumptionConstantSpeed - consumptionPulseAndGlide) / consumptionConstantSpeed;
end

function plotConstantSpeed (speed,rpmIce,lineColor,textColor)
  _powerSurplus = powerSurplus(speedometerToTrueSpeed(speed),rpmIce);
  [C,h] = contour(speed,rpmIce,_powerSurplus,[0,0]);
  clabel(C,h,'rotation',45);
  line = findobj(h,'String','0');
  set(line,'Color',textColor); 
  set(line,'String','Constant Speed');
  set(h,'LineWidth',3);
  set(h,'LineColor',lineColor);
end

function plotConstantSpeedByPowerIce (speed,powerIce,lineColor,textColor)
  _powerSurplus = powerSurplusForKmhAndPowerIce(speedometerToTrueSpeed(speed),powerIce);
  [C,h] = contour(speed,powerIce,_powerSurplus,[0,0]);
  clabel(C,h,'rotation',45);
  line = findobj(h,'String','0');
  set(line,'Color',textColor); 
  set(line,'String','Constant Speed');
  set(h,'LineWidth',3);
  set(h,'LineColor',lineColor);
end

function plotRpmMg1 (speed,rpmIce,rpmMg1,lineColor,textColor)
  _rpmMg1kmh = rpmMg1kmh(speedometerToTrueSpeed(speed),rpmIce);
  [C,h] = contour(speed,rpmIce,_rpmMg1kmh,[rpmMg1,rpmMg1]);
  clabel(C,h,'rotation',45,'LabelSpacing',170);
  line = findobj(h,'String',sprintf('%d',rpmMg1));
  set(line,'Color',textColor);
  set(line,'String',sprintf('MG1 %d rpm',rpmMg1));
  set(h,'LineWidth',3);
  set(h,'LineColor',lineColor);
end  

function plotRpmMg1ByPowerIce (speed,powerIce,rpmMg1,lineColor,textColor)
  _rpmMg1kmh = rpmMg1kmh(speedometerToTrueSpeed(speed),powerToRpmIce(powerIce));
  [C,h] = contour(speed,powerIce,_rpmMg1kmh,[rpmMg1,rpmMg1]);
  clabel(C,h,'rotation',45,'LabelSpacing',170);
  line = findobj(h,'String',sprintf('%d',rpmMg1));
  set(line,'Color',textColor);
  set(line,'String',sprintf('MG1 %d rpm',rpmMg1));
  set(h,'LineWidth',3);
  set(h,'LineColor',lineColor);
end  

function axisSpeedometerAndRpm ()
  global MAX_KMH_SPEEDOMETER MIN_RPM_ICE MAX_RPM_ICE;
  xlabel('Speed as displayed on dashboard (km/h)');
  set(gca,'XTick',0:10:MAX_KMH_SPEEDOMETER);
  ylabel('ICE RPM');
  set(gca,'YTick',[MIN_RPM_ICE MAX_RPM_ICE:-500:MIN_RPM_ICE]);
end

function axisSpeedometerAndPowerIce ()
  global MAX_KMH_SPEEDOMETER;
  xlabel('Speed as displayed on dashboard (km/h)');
  set(gca,'XTick',0:10:MAX_KMH_SPEEDOMETER);
  ylabel('ICE Power (kW)');
  set(gca,'YTick',[0:2:57]);
end

function axisSpeedometerAndTorqueMg2 ()
  global MAX_KMH_SPEEDOMETER MAX_TORQUE_MG2;
  xlabel('Speed as displayed on dashboard (km/h)');
  set(gca,'XTick',0:10:MAX_KMH_SPEEDOMETER);
  ylabel('MG2 Torque (Nm)');
  set(gca,'YTick',[0:50:MAX_TORQUE_MG2]);
end

function axisSpeedometerAndMg2Power ()
  global MAX_KMH_SPEEDOMETER MAX_TORQUE_MG2;
  xlabel('Speed as displayed on dashboard (km/h)');
  set(gca,'XTick',0:10:MAX_KMH_SPEEDOMETER);
  ylabel('MG2 Power (kW)');
  set(gca,'YTick',[0:1:50]);
end

function axisSpeedometerAndMg1Rpm ()
  global MAX_KMH_SPEEDOMETER;
  xlabel('Average speed as displayed on dashboard (km/h)');
  set(gca,'XTick',20:5:130);
  ylabel('MG1 RPM');
  set(gca,'YTick',[-1200:200:6000]);
end

function axisRpmAndTorqueInverter ()
  global MAX_RPM_MG1 MAX_TORQUE_MG2;
  xlabel('RPM');
  set(gca,'XTick',[0:500:MAX_RPM_MG1 MAX_RPM_MG1]);
  ylabel('Torque (Nm)');
  set(gca,'YTick',[0:50:MAX_TORQUE_MG2]);
end

function axisRpmAndTorqueMg1 ()
  global MAX_RPM_MG1 MAX_TORQUE_MG1;
  xlabel('MG1 RPM');
  set(gca,'XTick',[-MAX_RPM_MG1:2000:MAX_RPM_MG1]);
  ylabel('MG1 Torque (Nm)');
  set(gca,'YTick',[-MAX_TORQUE_MG1:5:MAX_TORQUE_MG1]);
end

function axisRpmAndTorqueMg2 ()
  global MAX_RPM_MG2 MAX_TORQUE_MG2;
  xlabel('MG2 RPM');
  set(gca,'XTick',[0:500:MAX_RPM_MG2 MAX_RPM_MG2]);
  ylabel('MG2 Torque (Nm)');
  set(gca,'YTick',[-MAX_TORQUE_MG2:50:MAX_TORQUE_MG2]);
end

function licence (x,y)
  hLicence = text(x,y,'Created by proprius (CC BY-NC-SA 4.0)','clipping','off');
  set(hLicence,'FontSize',6);
end

function setPaperUnits ()
  set(gcf,'PaperUnits','inches','PaperPosition',[0 0 9 7]);
end  

function plotEfficiency (speed,rpmIce)
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndRpm();
  title({MODEL, 'Efficiency map for ICE electric path MG1<->Inverter<->DC<->Inverter<->MG2, and transmission', '100 = no electric/mechanical losses and ICE most efficient'});
  _efficiency = efficiency(speedometerToTrueSpeed(speed),rpmIce);
  [C,h] = contourf(speed,rpmIce,_efficiency,[[20:2:75] [76:93]]);
  ch = clabel(C,h,'rotation',0);
  set(ch,'BackgroundColor',[1 1 .6])
  plotConstantSpeed(speed,rpmIce,'b','w');
  plotRpmMg1(speed,rpmIce,0,'g','w');
  plotRpmMg1(speed,rpmIce,800,[1.0 0.9 0.0],'w');
  licence(10,1020);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusEfficiency.png');
end

function plotEfficiencyByPowerIce (speed,powerIce)
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndPowerIce();
  title({MODEL, 'Efficiency map for ICE, electric path MG1<->Inverter<->DC<->Inverter<->MG2, and transmission', '100 = no electric/mechanical losses and ICE most efficient'});
  _efficiency = efficiencyForKmhAndPowerIce(speedometerToTrueSpeed(speed),powerIce);
  [C,h] = contourf(speed,powerIce,_efficiency,[[0:4:60] [60:2:94]]);
  ch = clabel(C,h,'rotation',0);
  set(ch,'BackgroundColor',[1 1 .6])
  plotConstantSpeedByPowerIce(speed,powerIce,'b','w');
  plotRpmMg1ByPowerIce(speed,powerIce,0,'g','w');
  %plotRpmMg1ByPowerIce(speed,powerIce,800,[1.0 0.9 0.0],'w');
  licence(10,-1);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusEfficiencyByPowerIce.png');
end

function plotElecticPathEfficiency (speed,rpmIce)
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndRpm();
  title({MODEL, 'Electric path efficiency'});
  _electricPathEfficiency = electricPathEfficiency(speedometerToTrueSpeed(speed),rpmIce) .* 100;
  [C,h] = contourf(speed,rpmIce,_electricPathEfficiency,[[38:4:74][74:2:82]]);
  ch = clabel(C,h,'rotation',0);
  set(ch,'BackgroundColor',[1 1 .6])
  plotConstantSpeed(speed,rpmIce,'b','w');
  plotRpmMg1(speed,rpmIce,0,'g','w');
  plotRpmMg1(speed,rpmIce,800,[1.0 0.9 0.0],'w');
  licence(10,1020);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusElectricPathEfficiency.png');
end

function plotElecticPathEfficiencyByPower (speed,powerIce)
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndPowerIce();
  title({MODEL, 'Electric path efficiency'});
  _electricPathEfficiency = electricPathEfficiencyForKmhAndPowerIce(speedometerToTrueSpeed(speed),powerIce) .* 100;
  [C,h] = contourf(speed,powerIce,_electricPathEfficiency,[[38:4:74][74:2:82]]);
  ch = clabel(C,h,'rotation',0);
  set(ch,'BackgroundColor',[1 1 .6])
  plotConstantSpeedByPowerIce(speed,powerIce,'b','w');
  plotRpmMg1ByPowerIce(speed,powerIce,0,'g','w');
  licence(10,-1);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusElectricPathEfficiencyByPower.png');
end

function plotConsumption (speed,rpmIce)
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndRpm();
  title({MODEL, 'Fuel consumption (L/100 km)'});
  _consumption = consumption(speedometerToTrueSpeed(speed),rpmIce);
  [C,h] = contourf(speed,rpmIce,_consumption,[0:30]);
  clabel(C,h,'rotation',0,'color','w');
  plotConstantSpeed(speed,rpmIce,'r','w');
  plotRpmMg1(speed,rpmIce,0,'g','w');
  plotRpmMg1(speed,rpmIce,800,[1.0 0.9 0.0],'w');
  licence(10,1020);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusConsumption.png');
end

function plotConsumptionByPowerIce (speed,powerIce)
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndPowerIce();
  title({MODEL, 'Fuel consumption (L/100 km)'});
  _consumption = consumptionByPowerIce(speedometerToTrueSpeed(speed),powerIce);
  [C,h] = contourf(speed,powerIce,_consumption,[0:30]);
  clabel(C,h,'rotation',0,'color','w');
  plotConstantSpeedByPowerIce(speed,powerIce,'r','w');
  plotRpmMg1ByPowerIce(speed,powerIce,0,'g','w');
  licence(10,-1);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusConsumptionByPowerIce.png');
end


function plotPowerAfterPsd (speed,rpmIce)
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndRpm();
  title({MODEL, 'Power after PSD (kW)'});
  _powerAfterPsd = powerAfterPsd(speedometerToTrueSpeed(speed),rpmIce);
  [C,h] = contourf(speed,rpmIce,_powerAfterPsd,[0:2:57]);
  clabel(C,h,'rotation',0);
  plotConstantSpeed(speed,rpmIce,'r','w');
  plotRpmMg1(speed,rpmIce,0,'g','w');
  plotRpmMg1(speed,rpmIce,800,[1.0 0.9 0.0],'w');
  licence(10,1020);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusPowerAfterPsd.png');
end

function plotElectricPathLoss (speed,rpmIce)
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndRpm();
  title({MODEL, 'Electric path loss (kW)'});
  _electricPathLossKw = electricPathLossKw(speedometerToTrueSpeed(speed),rpmIce);
  [C,h] = contourf(speed,rpmIce,_electricPathLossKw,[0:1:13]);
  clabel(C,h,'rotation',0,'color','w');
  plotConstantSpeed(speed,rpmIce,'r','w');
  plotRpmMg1(speed,rpmIce,0,'g','w');
  plotRpmMg1(speed,rpmIce,800,[1.0 0.9 0.0],'w');
  licence(10,1020);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusElectricPathLoss.png');
end

function plotGeneratorPower (speed,rpmIce)
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndRpm();
  title({MODEL, 'Generator power(kW)'});
  _generatorPower = generatorPower(speedometerToTrueSpeed(speed),rpmIce);
  [C,h] = contourf(speed,rpmIce,_generatorPower,[0:1:35]);
  clabel(C,h,'rotation',0,'color','w');
  plotConstantSpeed(speed,rpmIce,'r','w');
  plotRpmMg1(speed,rpmIce,0,'g','w');
  plotRpmMg1(speed,rpmIce,800,[1.0 0.9 0.0],'w');
  licence(10,1020);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusGeneratorPower.png');
end

function plotPowerSurplus (speed,rpmIce)
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndRpm();
  title({MODEL, 'Power surplus: Power after PSD minus power required for constant speed on flat ground (kW).'});
  _powerSurplus = powerSurplus(speedometerToTrueSpeed(speed),rpmIce);
  [C,h] = contourf(speed,rpmIce,_powerSurplus,[-40:2:40]);
  clabel(C,h,'rotation',0);
  plotConstantSpeed(speed,rpmIce,'b','k');
  plotRpmMg1(speed,rpmIce,0,[0.0 0.5 0.0],'k');
  plotRpmMg1(speed,rpmIce,800,'k','k');
  licence(10,1020);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusPowerSurplus.png');
end

function plotTorqueMg1 (speed,rpmIce)
  global MODEL MAX_TORQUE_MG1;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndRpm();
  title({MODEL, 'MG1 torque (Nm)'});
  _torqueMg1 = torqueMg1(rpmIce);
  [C,h] = contourf(speed,rpmIce,_torqueMg1,[-MAX_TORQUE_MG1:1:MAX_TORQUE_MG1]);
  clabel(C,h,'rotation',0);
  plotConstantSpeed(speed,rpmIce,'b','k');
  plotRpmMg1(speed,rpmIce,0,[0.0 0.5 0.0],'k');
  plotRpmMg1(speed,rpmIce,800,'k','k');
  licence(10,1020);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusTorqueMg1.png');
end

function plotTorqueMg2 (speed,rpmIce)
  global MODEL MAX_TORQUE_MG2;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndRpm();
  title({MODEL, 'MG2 torque (Nm)'});
  _torqueMg2 = torqueMg2(speedometerToTrueSpeed(speed),rpmIce);
  [C,h] = contourf(speed,rpmIce,_torqueMg2,[[-60:5:50] [60:10:200] 225 [250:50:MAX_TORQUE_MG2]]);
  clabel(C,h,'rotation',0,'Color','w');
  plotConstantSpeed(speed,rpmIce,'r','w');
  plotRpmMg1(speed,rpmIce,0,'g','w');
  plotRpmMg1(speed,rpmIce,800,[1.0 0.9 0.0],'w');
  licence(10,1020);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusTorqueMg2.png');
end

function plotEfficiencyInverter ()
  global MODEL MAX_TORQUE_MG2 MAX_RPM_MG1;
  mg2TorqueRange = 0:MAX_TORQUE_MG2;
  mg1RpmRange = 0:MAX_RPM_MG1;
  [mg1Rpm,mg2Torque] = meshgrid(mg1RpmRange,mg2TorqueRange);
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisRpmAndTorqueInverter();
  title({MODEL, 'Inverter efficiency'});
  _efficiencyInverter = efficiencyInverter(mg1Rpm,mg2Torque) .* 100;
  [C,h] = contourf(mg1Rpm,mg2Torque,_efficiencyInverter,[60:100]);
  clabel(C,h,'rotation',0,'Color','w');
  licence(10,-25);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusEfficiencyInverter.png');
end

function plotEfficiencyMg1 ()
  global MODEL MAX_TORQUE_MG1 MAX_RPM_MG1;
  mg1TorqueRange = -MAX_TORQUE_MG1:MAX_TORQUE_MG1;
  mg1RpmRange = -MAX_RPM_MG1:MAX_RPM_MG1;
  [mg1Rpm,mg1Torque] = meshgrid(mg1RpmRange,mg1TorqueRange);
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisRpmAndTorqueMg1();
  title({MODEL, 'MG1 efficiency'});
  _efficiencyMg1 = efficiencyMg1(mg1Rpm,mg1Torque) .* 100;
  [C,h] = contourf(mg1Rpm,mg1Torque,_efficiencyMg1,[60:2:100]);
  clabel(C,h,'rotation',0,'Color','w');
  licence(-9990,-49);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusEfficiencyMg1.png');
end

function plotEfficiencyMg2 ()
  global MODEL MAX_TORQUE_MG2 MAX_RPM_MG2;
  mg2TorqueRange = -MAX_TORQUE_MG2:MAX_TORQUE_MG2;
  mg2RpmRange = 0:MAX_RPM_MG2;
  [mg2Rpm,mg2Torque] = meshgrid(mg2RpmRange,mg2TorqueRange);
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisRpmAndTorqueMg2();
  title({MODEL, 'MG2 efficiency'});
  _efficiencyMg2 = efficiencyMg2(mg2Rpm,mg2Torque) .* 100;
  [C,h] = contourf(mg2Rpm,mg2Torque,_efficiencyMg2,[60:2:100]);
  clabel(C,h,'rotation',0,'Color','w');
  licence(0,-435);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusEfficiencyMg2.png');
end

function plotEfficiencyMg2AndInverter (kmhSamples)
  global MODEL MAX_TORQUE_MG2 MAX_KMH_SPEEDOMETER;
  mg2TorqueRange = 0:MAX_TORQUE_MG2;
  [speed,mg2Torque] = meshgrid(linspace(0,MAX_KMH_SPEEDOMETER,kmhSamples), mg2TorqueRange);
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndTorqueMg2();
  title({MODEL, 'MG2 efficiency'});
  mg2RpmRange = kmhToRpmMg2(speedometerToTrueSpeed(speed));
  _efficiencyMg2 = efficiencyMg2(mg2RpmRange,mg2Torque) .* efficiencyInverter(mg2RpmRange,mg2Torque) .* 100;
  [C,h] = contourf(speed,mg2Torque,_efficiencyMg2,[60:100]);
  clabel(C,h,'rotation',0,'Color','w');
  licence(10,-10);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusEfficiencyMg2AndInverter.png');
end

function plotEfficiencyMg2AndInverterPowerAxis (kmhSamples)
  global MODEL MAX_TORQUE_MG2 MAX_KMH_SPEEDOMETER;
  mg2PowerRange = 0:0.05:50;
  [speed,mg2Power] = meshgrid(linspace(0,MAX_KMH_SPEEDOMETER,kmhSamples), mg2PowerRange);
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndMg2Power();
  title({MODEL, 'MG2 efficiency'});
  mg2RpmRange = kmhToRpmMg2(speedometerToTrueSpeed(speed));
  _efficiencyMg2 = efficiencyMg2(mg2RpmRange,torqueFromPowerAndRpm(mg2Power,mg2RpmRange)) .* efficiencyInverter(mg2RpmRange,torqueFromPowerAndRpm(mg2Power,mg2RpmRange)) .* 100;
  [C,h] = contourf(speed,mg2Power,_efficiencyMg2,[60:100]);
  clabel(C,h,'rotation',0,'Color','w');
  licence(10,-1);
  hold off;
  setPaperUnits();
  saveas(hFigure,'PriusEfficiencyMg2AndInverterPowerAxis.png');
end

function plotTorqueAndPower (rpmIceRange)
  global MODEL MIN_RPM_ICE MAX_RPM_ICE ICE_TORQUE_TABLE;
  hFigure = figure;
  hold on;
  grid on;
  title(MODEL);
  [hAxes,hLine1,hLine2] = plotyy(rpmIceRange,rpmToTorqueIce(rpmIceRange),rpmIceRange,rpmToPowerIce(rpmIceRange),'plot');
  set(hLine1,'LineWidth',2);
  set(hLine2,'LineWidth',2);
  set(hAxes(1),'XTick',[MIN_RPM_ICE MAX_RPM_ICE:-500:MIN_RPM_ICE]);
  set(hAxes(2),'XTick',0);
  set(hAxes(1),'YLim', [60 120]);
  set(hAxes(2),'YLim', [0 60]);
  set(hAxes(1),'YTick',[60:5:120]);
  set(hAxes(2),'YTick',[0:5:60]);
  set(get(hAxes(1),'Xlabel'),'String','ICE RPM');
  set(get(hAxes(1),'Ylabel'),'String','ICE Torque (Nm)');
  set(get(hAxes(2),'Ylabel'),'String','ICE Power (kW)');
  
  colorTorque = 'b';
  colorPower = 'r';
  set(hAxes(1),'ycolor',colorTorque); % not working
  set(get(hAxes(1),'YLabel'),'color', colorTorque); 
  set(hLine1,'Color',colorTorque);
  set(hAxes(2),'ycolor',colorPower); % not working
  set(get(hAxes(2),'YLabel'),'color', colorPower); 
  set(hLine2,'Color',colorPower);
  
  hold off;
  set(gcf,'PaperUnits','inches','PaperPosition',[0 0 5.5 4]);
  saveas(hFigure,'PriusTorqueAndPower.png');
end

function plotPulseAndGlideSavings (kmhMiddle,kmhDelta,mg1Rpm,savings)  
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndMg1Rpm();
  title({MODEL, sprintf('Pulse and Glide savings (%%) compared to constant speed. Delta min/max speed = %d km/h', kmhDelta)});
  [C,h] = contourf(kmhMiddle,mg1Rpm,savings,[-10:2:50]);
  clabel(C,h,'rotation',-90,'Color',[0.85 0.0 0.0]);
  licence(20,-1600);
  hold off;
  setPaperUnits();
  saveas(hFigure,sprintf('PriusEfficiencyPulseAndGlideSavings%dkmh.png',kmhDelta));
end  

function plotPulseAndGlideConsumption (kmhMiddle,kmhDelta,mg1Rpm,consumption)  
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndMg1Rpm();
  title({MODEL, sprintf('Pulse and Glide consumption (L/100 km). Delta min/max speed = %d km/h', kmhDelta)});
  [C,h] = contourf(kmhMiddle,mg1Rpm,consumption,[2:0.1:6.5]);
  clabel(C,h,'rotation',-75,'Color',[0.85 0.0 0.0]);
  licence(20,-1600);
  hold off;
  setPaperUnits();
  saveas(hFigure,sprintf('PriusEfficiencyPulseAndGlideConsumption%dkmh.png',kmhDelta));
end  

function plotPulseAndGlide (kmhDelta)
  global MAX_KMH_SPEEDOMETER;
  %mg1RpmRange = [15:45]; kmhMiddleRange = [65.4:0.2:66.2];
  %mg1RpmRange = [120:10:150]; kmhMiddleRange = [88:93];
  %mg1RpmRange = [-1000:1:1000]; kmhMiddleRange = [20:0.2:120];
  %mg1RpmRange = [-1000:10:1000]; kmhMiddleRange = [20:0.2:120];
  %mg1RpmRange = [-1000:10:1000]; kmhMiddleRange = [20:120];
  mg1RpmRange = [-1200:100:6000]; kmhMiddleRange = [20:130];
  %mg1RpmRange = [-1000:1000:6000]; kmhMiddleRange = [20:10:130];
  [kmhMiddle,mg1Rpm] = meshgrid(kmhMiddleRange,mg1RpmRange);
  %[savings consumption] = arrayfun(@(kmhMiddle,mg1Rpm) pulseAndGlide(kmhMiddle,kmhDelta,mg1Rpm),kmhMiddle,mg1Rpm);
  [savings consumption] = pararrayfun(6,@(kmhMiddle,mg1Rpm) pulseAndGlide(kmhMiddle,kmhDelta,mg1Rpm),kmhMiddle,mg1Rpm);
  %plot(kmhMiddleRange, savings);
  %figure;
  %plot(kmhMiddleRange, consumption);
  plotPulseAndGlideSavings(kmhMiddle,kmhDelta,mg1Rpm,savings);
  plotPulseAndGlideConsumption(kmhMiddle,kmhDelta,mg1Rpm,consumption);
end

initialize();

powerIceRange = 2:0.02:57; kmhSamples = 1000; % high quality plot
%powerIceRange = 2:0.1:57; kmhSamples = 100; % mid quality plot
%powerIceRange = 2:57; kmhSamples = 10; % low quality plot
powerIce = repmat((powerIceRange)',1,kmhSamples);
rpmIceRange = powerToRpmIce( powerIceRange );

%rpmIceRange = MIN_RPM_ICE:.2:MAX_RPM_ICE; kmhSamples = 1000;  % very high quality plot
%rpmIceRange = MIN_RPM_ICE:20:MAX_RPM_ICE; kmhSamples = 100;   % mid quality plot
%rpmIceRange = MIN_RPM_ICE:1000:MAX_RPM_ICE; kmhSamples = 5;    % low quality plot
speed = linspace(minKmhSpeedometer(rpmIceRange),maxKmhSpeedometer(rpmIceRange),kmhSamples);
rpmIce = repmat((rpmIceRange)',1,kmhSamples);

%Z = powerMg1(speedometerToTrueSpeed(speed),rpmIce);                  [C,h] = contourf(speed,rpmIce,Z,[-30:2:30]);                     clabel(C,h,'rotation',0);
%Z = powerMg1relativeToIceKmh(speedometerToTrueSpeed(speed),rpmIce);  [C,h] = contourf(speed,rpmIce,Z,[0:0.1:2.5]);                    clabel(C,h,'rotation',0);
%Z = rpmMg1kmh(speedometerToTrueSpeed(speed),rpmIce);                 [C,h] = contour (speed,rpmIce,Z,[-MIN_RPM_MG1:400:MAX_RPM_MG1]); clabel(C,h,'rotation',0);

plotEfficiencyInverter();
plotEfficiencyMg1();
plotEfficiencyMg2();

plotTorqueAndPower(rpmIceRange);
%plotEfficiency(speed,rpmIce);
plotEfficiencyByPowerIce(speed,powerIce);
%plotElecticPathEfficiency(speed,rpmIce);
plotElecticPathEfficiencyByPower(speed,powerIce);
%plotElectricPathLoss(speed,rpmIce);
%plotConsumption(speed,rpmIce);
plotConsumptionByPowerIce(speed,powerIce);
%plotPowerAfterPsd(speed,rpmIce);
%plotPowerSurplus(speed,rpmIce);
%plotTorqueMg1(speed,rpmIce);
%plotTorqueMg2(speed,rpmIce);
%plotGeneratorPower(speed,rpmIce);

%plotEfficiencyMg2AndInverter(kmhSamples);
%plotEfficiencyMg2AndInverterPowerAxis(kmhSamples);

plotPulseAndGlide(10);
plotPulseAndGlide(20);
%plotPulseAndGlide(30); % Crashes!
 
pause