%
% LEAF motor imitation 
% P= 80 kW, T = 280 Nm Nmax = 10500 rpm
% Foula for IPM
%
figure;
close all;
clear all
rpm = 10500.; %motor max speed in rpm
w_max = rpm*2*pi/60 % motor speed in rad/s
T_b = 280 % motor torque
x=linspace(1,w_max,w_max); % rpm from 0 to max
y=linspace(0,T_b,T_b); % torque from 0 to T_b max
% Creating mesh
[w,T] = meshgrid(x,y);
% Torque saturation i.e. field weakening
for c=1:w_max
  for r=1:T_b
    if (T(r,c)*w(r,c)>80000);
      T(r,c) = T_b/w(r,c);
    end
  end
end
% output power from the motor
Output_power=(w.*T); % power = torque*rotation speed
% Calculating loss components. Coefficient fitted to meet Leaf measured map
base_loss=3700.;
L_1 = 0.0033*(base_loss); % -0.0033
L_w1_T0 =0.52*(w./w_max)*(base_loss); % 0.239
L_w0_T1 =0.02*(T/T_b)*base_loss; % 0.047
L_w1_T1 = -1.022*(T./T_b).*(w./w_max)*(base_loss-4000); % -1.022
L_w1_T2 = 1.3*(T./T_b).^2.*(w./w_max)*base_loss; % 1.071
L_w0_T2 = 0.103*(T./T_b).^2 *(base_loss); % 0.103
L_w2_T0 = -0.334*(w./w_max).^2*(base_loss); %-0.334
L_w0_T3 = 0.45*(T./T_b).^3*(base_loss); % 0.339
L_w2_T1 = 0.5*(T./T_b).*(w./w_max).^2*(base_loss); % 0.534
L_w3_T0 = 0.12*(w./w_max).^3*(base_loss+4000); %0.171
%
% Combining all loss components
%
P_loss=L_1+L_w1_T0+L_w0_T1+L_w1_T1+L_w1_T2+L_w0_T2...
+L_w2_T0+L_w0_T3+L_w2_T1+L_w3_T0;
Input_power = Output_power+P_loss;
%
%Creating efficiency map
%
map_mot=Output_power./Input_power;
%
% Further modifiy map to meet target
%
map_mot = map_mot-0.01;
map_mot(map_mot>0.9642)=0.97;
map_mot(map_mot<=0.5)=0.5;
hold on;
%
% Set contour lines to be drawnow
%
V1=(0:0.05:.9);
V2=(0.9:0.02:1);
V=[V1,V2];
box off
grid off
%
% Plot using contour plot and color map
%
contourf (w,T,map_mot,V);
V=[V1,V2,80000];
axis([0 w_max 0 T_b])
xlabel('speed [rad/s]');
ylabel('Torque [N.m]');
title('LEAF motor efficiency map')
hold off
%