close all
clear all

filenm = 'Mar22-2021_EffMap_AveragedQuadrants';
%filenm = 'Mar22-2021_EffMap_AveragedTorqueDir';

load([filenm '.mat']);



[X,Y] = meshgrid(-Speeds,WindingCurrents);  %Speeds specify driving motor speed, so WEC motor measures velocity with opposite sign.

MechPower = -(2*pi/60)*meanRPM.*meanAuxTorque;
ElecPower = meanVBus.*(meanIBatt+meanLoadDumpCurrent);

Loss = MechPower-ElecPower;


%Compute Loss from paramterized model
LossModel = zeros(size(Loss));  

tau_c = 0.1; k_v = 0.03/1000; k_th = 20;
for i = 1:length(Speeds)
    for j = 1:length(WindingCurrents)
        N = meanRPM(i,j);
FrictionLoss = ((tau_c*tanh(2*pi*N/60/k_th)+k_v*N/10)*N);
ISquaredRLoss = .18*meanIwind(i,j)^2;
SwitchingLoss = 3+0.1*meanVBus(i,j);

LossModel(i,j) = FrictionLoss+ISquaredRLoss+SwitchingLoss;

    end
end


Eff = NaN(size(ElecPower));
for i = 1:size(Speeds)
    for j = 1:size(WindingCurrents)
        if(ElecPower(i,j) >= 0)  %Generating
            Eff(i,j) = ElecPower(i,j)./(MechPower(i,j));
        end
        if(ElecPower(i,j) < 0)  %Motoring
            Eff(i,j) = MechPower(i,j)./ElecPower(i,j);
        end
    
    if((Speeds(i) == 0) || (WindingCurrents(j) == 0))  
        Eff(i,j) = 0;
    end
    
    if(Eff(i,j) <= 0)
        Eff(i,j) = NaN;
    end
    
  end
end


figure('Position',[1000    600   1500     900]);
h1 = plot(X,Y,'ro');
grid on;
hold all
h2 = plot(meanRPM,meanIwind,'kx');
xlabel('RPM');
ylabel('Winding Current (Amps)');
legend([h1(1) h2(1)],'Requested Speed/Current','Measured Speed/Current');
xlim([-6000 6000]);
ylim([min(WindingCurrents) max(WindingCurrents)]);
title('Test Points in Speed-Torque Plane');

figure('Position',[ 1000         300        1000        1200]);
subplot(3,2,1)
surf(meanRPM,meanIwind,MechPower);
xlabel('RPM');  ylabel('Winding Current (Amps)');  zlabel('Power (W)');
title('Mechanical Power');
subplot(3,2,2)
surf(meanRPM,meanIwind,ElecPower);
xlabel('RPM');  ylabel('Winding Current (Amps)');  zlabel('Power (W)');
title('Electrical Power');

subplot(3,2,3)
plot(meanRPM(:,1:2:end),MechPower(:,1:2:end),'-x');
xlabel('RPM');    ylabel('Power (W)');
title('Mechanical Power');
grid on;
j = 0;
for i = 1:length(WindingCurrents)
   j = j+1;
   legstr{j} = [num2str(WindingCurrents(i)) 'A'];
end
%   legend(legstr,'Location','NorthWest');

 subplot(3,2,4)
plot(meanRPM(:,1:end),ElecPower(:,1:end),'-x');
xlabel('RPM');    ylabel('Power (W)');
grid on;
title('Electrical Power');
j = 0;
for i = 1:length(WindingCurrents)
   j = j+1;
   legstr{j} = [num2str(WindingCurrents(i)) 'A'];
end
 %  legend(legstr,'Location','NorthWest');
   
   
subplot(3,2,5)
plot(meanIwind',MechPower','-x');
xlabel('Wind Current (Amps)');    ylabel('Power (W)');
title('Mechanical Power');

grid on;
clear legstr;
for i = 1:length(Speeds)
    legstr{i} = [num2str(Speeds(i)) ' RPM'];
end
 %  legend(legstr);
   
      
subplot(3,2,6)
plot(meanIwind',ElecPower','-x');
xlabel('Wind Current (Amps)');    ylabel('Power (W)');
grid on;
title('Electrical Power');
clear legstr;
for i = 1:length(Speeds)
    legstr{i} = [num2str(Speeds(i)) ' RPM'];
end
 %  legend(legstr);
 
 
 
figure('Position',[1000    600   1500     900]);
subplot(2,2,[1 2])
surf(meanRPM,meanIwind,Loss);
xlabel('RPM');
ylabel('Winding Current (A)');
zlabel('Loss (W)');
title('Total Electrical Loss (W)');

subplot(2,2,3)
plot(meanRPM(:,1:end),Loss(:,1:end),'-x');
xlabel('RPM');    ylabel('Loss (W)');
grid on;

subplot(2,2,4)
plot(meanIwind',Loss','-x');
xlabel('Wind Current (Amps)');    ylabel('Loss (W)');
grid on;

figure('Position',[1000    600   1500     900]);
subplot(2,2,[1 2])
surf(meanRPM,meanIwind,LossModel);
xlabel('RPM');
ylabel('Winding Current (A)');
zlabel('Loss (W)');
title('Total Electrical Loss (W) Parametrized Model');

subplot(2,2,3)
plot(meanRPM(:,1:end),LossModel(:,1:end),'-x');
xlabel('RPM');    ylabel('Loss (W)');
grid on;

subplot(2,2,4)
plot(meanIwind',LossModel','-x');
xlabel('Wind Current (Amps)');    ylabel('Loss (W)');
grid on;


 
figure('Position',[1000    600   1500     900]);
subplot(2,2,[1 2])
surf(meanRPM,meanIwind,Eff);
xlabel('RPM');
ylabel('Winding Current (A)');
zlabel('Eff');
title('Total Efficiency');

subplot(2,2,3)
plot(meanRPM(:,1:end),Eff(:,1:end),'-x');
xlabel('RPM');    ylabel('Eff');
grid on;

subplot(2,2,4)
plot(meanIwind',Eff','-x');
xlabel('Wind Current (Amps)');    ylabel('Eff');
grid on;
