close all
clear all

%filenm = 'May17-2021_EffMapp_PosToNegTorque_Scale1';
%filenm = 'Mar22-2021_EffMapp_AveragedQuadrants';
%filenm = 'May3-2022_EffMapp_PosToNegTorque';
filenm = 'June10_2022_UWEff';

load([filenm '.mat']);



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

Rll=0.9*2; %motor resistance line-line
Rwl=0.025; %motor lead resistance per leg
Rq=.05; %transistor resistance
WindingPower = meanIwind.^2*Rll/2*1.5;
LinePower = meanIwind.^2*Rwl*1.5;
QPower = meanIwind.^2*Rq*1.5;

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

avgSwitchPower=mean(ElecPowerReal(11,[1:21]));

ElectronicPower = ElecPowerReal.-QPower;

Loss = MechPower+ElecPower;

Eff = NaN(size(ElecPower));

m1 = polyfit(meanRPM([2:11],5),MechPower([2:11],5),1);
m2 = polyfit(meanRPM([11:21],5),MechPower([11:21],5),1);
m3 = polyfit(meanRPM([2:11],17),MechPower([2:11],17),1);
m4 = polyfit(meanRPM([11:21],17),MechPower([11:21],17),1);
m5 = polyfit(meanRPM([12:21],5),meanAuxTorque([12:21],5),1);
m6 = polyfit(meanRPM([2:10],5),meanAuxTorque([2:10],5),1);

SealTorque= m6(2)-m5(2);
SealPower = ((2*pi/60)*abs(meanRPM).*SealTorque);
MechNetPower = MechPower.+SealPower;
MechSealPower=MechPower;
MechSealPower([11:21],:)=MechSealPower([11:21],:).+SealPower([11:21],:);
straightline=m1(1).*meanRPM(:,5).+m1(2);

for i = 1:size(ElecPower,1)
    for j = 1:size(ElecPower,2)
        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([-120 120]);
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:end),MechPower(:,1:end),'-x'); %original
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');
%meanRPM
%ElecPower
 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');
   
   %MechPower
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 Losses (W)');

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

%Loss
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,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;

figure('Position',[ 1000         300        1000        1200]);

subplot(3,2,1)
plot(meanIwind(11,[1:21]),ElecPower(11,[1:21]),'-x'); %total electrical loss at 0rpm
xlabel('Winding Current (Amps)');    ylabel('Power (W)');
title({'Electrical Power Loss at 0RPM';'Plot 1'});

subplot(3,2,2)
plot(meanIwind(11,[1:21]),WindingPower(11,[1:21]),'-x'); %winding i2r loss at 0rpm
xlabel('Winding Current (Amps)');    ylabel('Power (W)');
WindingPowertitle = sprintf("Resistance = %d Ohms L-L",Rll);
title({'Winding i2r Loss at 0RPM';WindingPowertitle;'Plot 2'});

subplot(3,2,3)
plot(meanIwind(11,[1:21]),QPower(11,[1:21]),'-x'); %transistor loss at 0rpm
xlabel('Winding Current (Amps)');    ylabel('Power (W)');
QPowertitle = sprintf("Resistance = %d Ohms",Rq);
title({'Transistor I2R Loss at 0RPM';QPowertitle;'Plot 3'});

subplot(3,2,4)
plot(meanIwind(11,[1:21]),LinePower(11,[1:21]),'-x'); %hookup wire loss at 0rpm
xlabel('Winding Current (Amps)');    ylabel('Power (W)');
LinePowertitle = sprintf("Resistance = %d Ohms",Rwl);
title({'Hookup wire Loss at 0RPM';LinePowertitle;'Plot 4'});

subplot(3,2,5)
plot(meanIwind(11,[1:21]),ElecPowerReal(11,[1:21]),'-x'); %switching loss at 0rpm
xlabel('Winding Current (Amps)');    ylabel('Power (W)');
avgSwitchPowertitle = sprintf("Mean = %d",avgSwitchPower);
title({'Switching Power Loss at 0RPM';avgSwitchPowertitle;'Plot 5'});

subplot(3,2,6)
plot(meanIwind(11,[1:21]),ElectronicPower(11,[1:21]),'-x'); %electronic loss at 0rpm
xlabel('Winding Current (Amps)');    ylabel('Power (W)');
title({'Electronic Power Loss at 0RPM';'Plot 6'});

figure('Position',[ 1000         300        1000        1200]);

subplot(3,2,1)
plot(meanRPM([1:21],11),MechPower([1:21],11),'-x'); %seal loss 0Amps
xlabel('RPM');    ylabel('Power (W)');
title({'Mechanical Power Loss at 0amp';'Plot 1'});

subplot(3,2,2)
plot(meanRPM([1:21],11),MechNetPower([1:21],11),'-x'); %seal loss 0Amps
xlabel('RPM');    ylabel('Power (W)');
SealTorquestr = sprintf("Seal Torque = %d Nm",SealTorque);
title({'Mechanical Power Loss - Seal Drag at 0amp';SealTorquestr;'Plot 2'});

figure('Position',[ 1000         300        1000        1200]);

subplot(3,2,2)
plot(meanRPM([2:21],5),MechPower([2:21],5),meanRPM([2:21],5),MechSealPower([2:21],5),meanRPM([2:21],5),straightline([2:21]));
xlabel('RPM');    ylabel('Power (W)');
%slopes1 = sprintf("Slope -100-0RPM =%d 0-100RPM = %d",m2(1),m1(1));
%slopes2 = sprintf("Slope -100-0RPM =%d 0-100RPM = %d",m3(1),m4(1));
%title({'Mechanical Power';slopes1;slopes2});
title({'Mechanical power';'with Seal Drag included from -100 to 0 RPM';SealTorquestr;'Plot 2'});
grid on;
legend('Measured','adjusted','line');

subplot(3,2,1)
plot(meanRPM([1:21],[5]),meanAuxTorque([1:21],[5]),'-x');
xlabel('RPM');    ylabel('Torque (Nm)');
slopes1 = sprintf("Slope -100-0RPM =%.2d 0-100RPM = %.2d",m5(1),m6(1));
offset1 = sprintf("Offset -100-0RPM =%3d 0-100RPM = %3d",m5(2),m6(2));
title({'Torque Asymmetry @12Amp winding current';slopes1;offset1;'Plot 1'});
grid on;

