% This script is used to estimate the efficiency map
% from the 
% motor namaplate parameters
close all; clear all; clc;
% This algorithm is for Permanent magnet Synchronnous Motor (PMSM)

%vozidlo  ENS 12 var.17

U1n = 34;              % nominal voltage V   
U1G=U1n;
I1n = 205;              % nominal current A
cosfin = 0.92;           % nominal power factor 
Pn = 10000;            % nominal power W
R1 = 1.3*2.5e-3;        % stator winding resistance Ohm, 1.3 is the coeffitient to account the increase of stator resistance due to heating
nn = 10000;              % nominal speed min-1
Tn = 9.5;            % nominal torque Nm
nmax=15840;              % maximal revolution min-1
nmin=100;               % maximal revolution min-1
Tmaxrel=3;           % pøetížitelnost motoru



% calculate nominal losses
dPn = sqrt(3)*U1n*I1n*cosfin - Pn;
% calculate nominal Joule losses
%dPjn = 3*R1*I1n^2;
dPjn=0.5*dPn;
% calculate nominal no load losses
%dP0n = dPn - dPjn;
dP0n=0.5*dPn;

disp(['dPn=',num2str(dPn),' dPjn=',num2str(dPjn),' dP0n=',num2str(dP0n)]);

Pmax=Tmaxrel*Pn;
P = linspace (1000,Pmax,15);
n = linspace (nmin,nmax,26);

% prepare vectors for map calculations
Px = P;
for i = 1:length(n)-1
   Px = [Px ;P];
end

nx = n;
for i = 1:length(P)-1
   nx = [nx ;n];
end

nx = transpose(nx);

dP = [zeros(length(n),length(P))];
eta = [zeros(length(n),length(P))];

% calculate efficiency map for all points
for i = 1:length(P)
    for j = 1:length(n)
         Tij = P(i)*60/(2*pi*n(j));
         
        %         losses calculation  
         U1=U1n*n(j)/nn;
         if U1>U1n
            U1=U1n;
         end
         dP(j,i) = dP0n * (U1/U1n)^2 * nn/n(j) + dPjn *((Tij*n(j)*U1n)/(Tn*nn*U1))^2  ;      
        %         efficiency calculation for motor
        eta(j,i) = P(i)/(P(i) + dP(j,i));
        I1=(Tij*2*pi*n(j)/60)/U1;
        
     
        disp(['P(i)=',num2str(P(i)),' Tij=',num2str(Tij),' n(j)=',num2str(n(j)),' U1=',num2str(U1)]);
        disp(['dP(j,i)=',num2str(dP(j,i)),' eta(j,i)=',num2str(eta(j,i))]);
        if P(i)<=Pn&&Tij<=Tn
          disp([' I1=',num2str(I1)]); 
        end
    end
end

% plot as a contour map
[C,h] = contour3(nx,Px/1000,eta*100,10,'k');

clabel(C,h,'LineStyle','--','FontSize',10);
h.LineStyle = '--';
% label only those selected iso countour lines

%h.LevelList = [0.15,0.3,0.45,0.6,0.75,0.8,0.82,0.84,0.86,0.88,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97];
h.LevelList = [15,30,45,60,75,80,82,84,86,88,90,91,92,93,94,95,96,97];


% add surface chart of the efficiency map
hold on;
gf = surf(nx,Px/1000,eta*100);

% make the surface chart transparent so that the iso contour labels are not
% hidden
gf.FaceAlpha = 0.6;

% interpolate the surface chart
shading interp

% add description
title 'Mapa úèinností - Motorický režim'
c = colorbar('southoutside');
c.Label.String = 'Efficiency [%]';
xlabel('{\itOtáèky} [1/min]')
ylabel('{\itVýkon} [kW]')
zlabel('{\itÚèinnost} [-]')
grid on

% pridano
  % sikma primka
  n1=(30*Pn/1000)/(pi*Tn);
  xprimka=[nmin,nn];
  yprimka=[nmin*Pn/(nn*1000),Pn/1000];
  zprimka=[1,1];
    
% vodorovna primka
    xvprimka=[nn,nmax];
    yvprimka=[Pn/1000,Pn/1000];
    zvprimka=[1,1];
   xvch=[xprimka,xvprimka];
   yvch=[yprimka,yvprimka];
   zvch=[zprimka,zvprimka];
        
   plot3(xvch,yvch,zvch,'-k','LineWidth',3)
   view([0,90])
   
   
% Ulozeni
    savefig('Charakteristika_motoru.fig');
    xlswrite('Charakteristika_motoru_p.xls',Px);
    xlswrite('Charakteristika_motoru_n.xls',nx);
    xlswrite('Charakteristika_motoru_T.xls',Tij);

