clear all;
close all;

U= 5:48;            %motor terminal voltage
M = 0:100:5000;     %Load torque
M(1) = 10.0;            %Start at 10 mNm

pumpVolumePerRev = .865 /1000     %Liters / rev
pumpEff = .65   
revsPerLiter = 1 / (pumpVolumePerRev * pumpEff);

n = zeros(4, length(U), length(M));
I = zeros(4, length(U), length(M));
eta = zeros(4, length(U), length(M));
joulesPerLiter = zeros(4, length(U), length(M));
litersPerJoule = zeros(4, length(U), length(M));

%Maxon EC90 244879 48 Volt
kn = 44;    %Speed constant rpm/Volt
R = 2.28;   %Terminal resistance ohms
kM = 217;   %Torque constant mNm/A
Io = 135;   %No load current mA
Mr = kM * (Io / 1000); %Friction torque
Rexternal = .2  %ohms

%calculate speed-torque curves
%calcualte efficiency curves
for i=1:length(U)
    for j=1:length(M)
        I(1,i,j) = M(j) / kM;
        Umotor = U(i) - (I(1,i,j)*Rexternal);
        n(1,i,j) = (kn*Umotor)-((30000/pi())*(R/kM^2)*M(j));
        if(n(1,i,j) < 0)
            n(1,i,j) = 0;
        end
        eta(1,i,j) = (pi()/30000) * (n(1,i,j) * (M(j) - Mr)) / ((U(i) * I(1,i,j)));
        if(eta(1,i,j) < 0)
            eta(1,i,j) = 0;
        end
        joulesPerLiter(1,i,j) = I(1,i,j) * U(i) * (revsPerLiter *60 / n(1,i,j));
        litersPerJoule(1,i,j) = (revsPerLiter*60.0)/(I(1,i,j) * U(i));
    end
end

% surf(M, U, squeeze(n(1,:,:)));
% figure
% surf(M, U, squeeze(eta(1,:,:)));

%Maxon EC90 429271 36 Volt
kn = 88;    %Speed constant rpm/Volt
R = .522;   %Terminal resistance ohms
kM = 109;   %Torque constant mNm/A
Io = 348;   %No load current mA
Mr = kM * (Io / 1000); %Friction torque

%calculate speed-torque curves
%calcualte efficiency curves
for i=1:length(U)
    for j=1:length(M)
        I(2,i,j) = M(j) / kM;
        Umotor = U(i) - (I(2,i,j)*Rexternal);
        n(2,i,j) = (kn*Umotor)-((30000/pi())*(R/kM^2)*M(j));
        if(n(2,i,j) < 0)
            n(2,i,j) = 0;
        end
        eta(2,i,j) = (pi()/30000) * (n(2,i,j) * (M(j) - Mr)) / ((U(i) * I(2,i,j)));
        if(eta(2,i,j) < 0)
            eta(2,i,j) = 0;
        end
        joulesPerLiter(2,i,j) = I(2,i,j) * U(i) * (revsPerLiter *60 / n(2,i,j));
    end
end

% figure
% surf(M, U, squeeze(n(1,:,:)));
% figure
% surf(M, U, squeeze(eta(1,:,:)));

%Maxon EC90 323772 24 Volt
kn = 135;    %Speed constant rpm/Volt
R = .343;   %Terminal resistance ohms
kM = 70.5;   %Torque constant mNm/A
Io = 544;   %No load current mA
Mr = kM * (Io / 1000); %Friction torque

%calculate speed-torque curves
%calcualte efficiency curves
for i=1:length(U)
    for j=1:length(M)
        I(3,i,j) = M(j) / kM;
       Umotor = U(i) - (I(3,i,j)*Rexternal);
        n(3,i,j) = (kn*Umotor)-((30000/pi())*(R/kM^2)*M(j));
        if(n(3,i,j) < 0)
            n(3,i,j) = 0;
        end
        eta(3,i,j) = (pi()/30000) * (n(3,i,j) * (M(j) - Mr)) / ((U(i) * I(3,i,j)));
        if(eta(3,i,j) < 0)
            eta(3,i,j) = 0;
        end
        joulesPerLiter(3,i,j) = I(3,i,j) * U(i) * (revsPerLiter *60 / n(3,i,j));
    end
end

%Maxon RE65 70V 353300
kn = 44.7;    %Speed constant rpm/Volt
R = 1.06;   %Terminal resistance ohms
kM = 214;   %Torque constant mNm/A
Io = 153;   %No load current mA
Mr = kM * (Io / 1000); %Friction torque

%calculate speed-torque curves
%calcualte efficiency curves
for i=1:length(U)
    for j=1:length(M)
        n(4,i,j) = (kn*U(i))-((30000/pi())*(R/kM^2)*M(j));
        if(n(4,i,j) < 0)
            n(4,i,j) = 0;
        end
        I(4,i,j) = M(j) / kM;
        eta(4,i,j) = (pi()/30000) * (n(4,i,j) * (M(j) - Mr)) / ((U(i) * I(4,i,j)));
        if(eta(4,i,j) < 0)
            eta(4,i,j) = 0;
        end
        joulesPerLiter(4,i,j) = I(4,i,j) * U(i) * (revsPerLiter *60 / n(4,i,j));
    end
end

%Plot Efficiency
surf(M,U,squeeze(eta(1,:,:)))
title('Efficiency 48 Volt Winding');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('Efficiency');

figure
surf(M,U,squeeze(eta(2,:,:)))
title('Efficiency 36 Volt Winding');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('Efficiency');

figure
surf(M,U,squeeze(eta(3,:,:)))
title('Efficiency 24 Volt Winding');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('Efficiency');

figure
surf(M(1:10),U(1:44),squeeze(eta(4,1:end,1:10)))
title('Efficiency RE65');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('Efficiency');

%Plot Current
figure
surf(M(1:10),U(1:44),squeeze(I(1,1:end,1:10)))
title('Current 48 Volt Winding');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('Current');

figure
surf(M(1:10),U(1:44),squeeze(I(2,1:end,1:10)))
title('Current 36 Volt Winding');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('Current');

figure
surf(M(1:10),U(1:44),squeeze(I(3,1:end,1:10)))
title('Current 24 Volt Winding');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('Current');

figure
surf(M(1:10),U(1:44),squeeze(I(4,1:end,1:10)))
title('Current RE65');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('Current');

%Plot RPM
figure
surf(M(1:10),U(1:44),squeeze(n(1,1:end,1:10)))
title('RPM 48 Volt Winding');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('RPM');

figure
surf(M(1:10),U(1:44),squeeze(n(2,1:end,1:10)))
title('RPM 36 Volt Winding');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('RPM');

figure
surf(M(1:10),U(1:44),squeeze(n(3,1:end,1:10)))
title('RPM 24 Volt Winding');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('RPM');

figure
surf(M(1:10),U(1:44),squeeze(n(4,1:end,1:10)))
title('RPM RE65');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('RPM');

%Plot Joules per liter
figure
surf(M(1:10),U(24:32),squeeze(joulesPerLiter(1,24:32,1:10)))
title('Liters per Joule 48 Volt Winding');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('Liters per Joule');

figure
surf(M,U,squeeze(joulesPerLiter(2,:,:)))
title('Joules per Liter 36 Volt Winding');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('Joules per Liter');

figure
surf(M,U,squeeze(joulesPerLiter(3,:,:)))
title('Joules per Liter 24 Volt Winding');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('Joules per Liter');

figure
surf(M,U,squeeze(joulesPerLiter(4,:,:)))
title('Joules per Liter RE65');
xlabel('Torque (mNm');
ylabel('Motor Terminal Voltage');
zlabel('Joules per Liter');






