close all
clear all

TargetCurve = load('PowerVsRPM.mat');

Ntuned = 0*1500;

N = 10:6000;

Rfets = 0.025;
%Rfets = 0.6;

name{1} = 'AKM52-H';                    Kv(1) = 92.7; R(1) = 2.35; L(1) = 11.9;
name{2} = 'AKM52-K';                    Kv(2) = 60.1; R(2) = 0.93; L(2) = 5;
name{3} = 'AKM52-L';                    Kv(3) =	48.3; R(3) = 0.61; L(3) = 3.24;
name{4} = 'AKM52-M';                    Kv(4) =	42.4; R(4) = 0.49; L(4) = 2.5;
name{5} = 'AKM52-S-6';                  Kv(5) =	40.4; R(5) = 0.51; L(5) = 2.27;	
name{6} = 'AKM52-S';                    Kv(6) =	23.7; R(6) = 0.16; L(6) = 0.78;
		
name{7} = 'AKM53-L';                    Kv(7) =	63.6; R(7) = 0.69; L(7) = 3.64;
name{8} = 'AKM53-M';                    Kv(8) =	54.7; R(8) = 0.51; L(8) = 2.7;
name{9} = 'AKM53-P';                    Kv(9) =	38.4; R(9) = 0.28; L(9) = 1.33;	
name{10} = 'AKM53 (16Tof 19.5bi-MW)';    Kv(10) = 23.7; R(10) = 0.13; L(10) = 0.5;
		
name{11} = 'AKM54-N';                    Kv(11) =	51.3; R(11) = 0.33; L(11) = 1.8;
name{12} = 'AKM54 (23Tof 20bi-MW)';     Kv(12) = 45.3; R(12) = 0.29; L(12) = 1.37;
	
name{13} = 'AKM54 (16Tof 19.5bi-MW)';   Kv(13) = 31.5; R(13) = 0.19; L(13) = 0.66;

factor = 1;
name{14} = ['AKM62Q div ' num2str(factor)];   Kv(14) = 35.5/factor; R(14) = 0.24/factor^2; L(14) = 1.6/factor^2;
factor = 1;
name{15} = ['AKM65P div ' num2str(factor)];   Kv(15) = 80.5/factor; R(15) = 0.37/factor^2; L(15) = 2.8/factor^2;


factor = 4;
name{18} = ['AKM74Q div ' num2str(factor)];   Kv(18) = 129/factor; R(18) = 0.25/factor^2; L(18) = 3.8/factor^2;
factor = 2;
name{19} = ['AKM74Q div ' num2str(factor)];   Kv(19) = 80.5/factor; R(19) = 0.37/factor^2; L(19) = 2.8/factor^2;


name{20} = ['AKM64 (18tpc of 18bi)'];   Kv(20) = 40; R(20) = 0.15; L(20) = 0.87;
name{21} = ['AKM63 (24tpc of 19bi)'];   Kv(21) = 40; R(21) = 0.19; L(21) = 1.22;
name{22} = ['AKM62 (36tpc of 21-bifilar)'];   Kv(22) = 40; R(22) = 0.29; L(22) = 2.01;




for mot_num = [20] %7:13


for VoltageTarget = [250];



Rwinding = R(mot_num);
Lwinding = L(mot_num)/1000;


Cwinding = (1/Lwinding)*(12/(2*pi*Ntuned))^2;

Vg = N*1.35*Kv(mot_num)/1000;


figure('Position',[1000 66 917 1432],'PaperSize',[11 17],'PaperPosition',[.5 .5 10.5 16.5]);

i = 0;
for PowerTargetFactor = [.25 .5 .75 1 2 4];
i = i+1;
legend_str{i} = [num2str(PowerTargetFactor) ' W/RPM'];

PowerTarget = PowerTargetFactor * N;

if(Ntuned > 0)
Rl = ((Rwinding+2*Rfets).^2 +(2*pi*Lwinding*N/12-12./(2*pi*N*Cwinding)).^2).^.5;
else
Rl = ((Rwinding+2*Rfets).^2 +(2*pi*Lwinding*N*5/60).^2).^.5;
end
Resistance = Rwinding+2*Rfets
Reactance = 2*pi*Lwinding*N*5/60

Reff = (VoltageTarget.^2)./PowerTarget;

ReqBoost = VoltageTarget./Vg;

ReqDutyCycle = 0.5*Vg.*(1+(1-4*Rl*(VoltageTarget^2)./(Reff.*Vg.^2)).^.5)/VoltageTarget;

Rl2 = Rwinding+2*Rfets;

Eff = 1./(1+Rl./(Reff.*(ReqDutyCycle.^2)));

ah1 = subplot(5,1,1)
plot(N,PowerTarget./N)
hold all
xlabel('RPM');
ylabel('Specific Power (W/RPM)')

jdx = find(ReqDutyCycle>1);
ah2 = subplot(5,1,2)
plot(N,ReqBoost,N(jdx),ReqBoost(jdx),'rx')
hold all
xlabel('RPM');
ylabel('Required Boost Factor')

idx = find(imag(ReqDutyCycle) > 0);
ah3 = subplot(5,1,3)
plot(N,ReqDutyCycle,N([idx jdx]),ReqDutyCycle([idx jdx]),'r.')
hold all
xlabel('RPM');
ylabel('Required Boost Duty Cycle')

ah4 = subplot(5,1,4)
plot(N,Eff,N([idx jdx]),Eff([idx jdx]),'r.')
hold all
xlabel('RPM');
ylabel('Efficiency')
end


%Compute same for three power targets (P, Pmin, Pmax)

legend_str{i+1} = ['Nominal Power Line'];
legend_str{i+2} = ['Max Power Line'];
legend_str{i+3} = ['Min Power Line'];

arg = TargetCurve.P./TargetCurve.RPM; arg(1) = 0;
Ptarg(1,:) = interp1(TargetCurve.RPM,arg,N).*N;
arg = TargetCurve.Pmax./TargetCurve.RPM ; arg(1) = 0;
Ptarg(2,:) = interp1(TargetCurve.RPM,arg,N).*N;
arg = TargetCurve.Pmin./TargetCurve.RPM ; arg(1) = 0;
Ptarg(3,:) = interp1(TargetCurve.RPM,arg,N).*N;

%figure
for i = 1:3

    
if(Ntuned > 0)
Rl = ((Rwinding+2*Rfets).^2 +(2*pi*Lwinding*N/12-12./(2*pi*N*Cwinding)).^2).^.5;
else
Rl = ((Rwinding+2*Rfets).^2 +(2*pi*Lwinding*N*5/60).^2).^.5;
end

Reff = (VoltageTarget.^2)./Ptarg(i,:)
ReqBoost = VoltageTarget./Vg;
ReqDutyCycle = 0.5*Vg.*(1+(1-4*Rl*(VoltageTarget^2)./(Reff.*Vg.^2)).^.5)/VoltageTarget;

Rl2 = Rwinding+2*Rfets;
Eff = 1./(1+Rl./(Reff.*(ReqDutyCycle.^2)));
   

ah1 = subplot(5,1,1)
plot(N,Ptarg(i,:)./N,'k--')
hold all
xlabel('RPM');
ylabel('Specific Power (W/RPM)')


jdx = find(ReqDutyCycle>1);
ah2 = subplot(5,1,2)
plot(N,ReqBoost,'k-');
hold all
plot(N(jdx),ReqBoost(jdx),'rx')
hold all
xlabel('RPM');
ylabel('Required Boost Factor')

idx = find(imag(ReqDutyCycle) > 0);
ah3 = subplot(5,1,3)
plot(N,ReqDutyCycle,'k--');
hold all;
plot(N([idx jdx]),ReqDutyCycle([idx jdx]),'r.')
hold all
xlabel('RPM');
ylabel('Required Boost Duty Cycle')

ah4 = subplot(5,1,4)
plot(N,Eff,'k--');
hold all;
plot(N([idx jdx]),Eff([idx jdx]),'r.');
hold all



end


ah5 = subplot(5,1,5)
plot(N,(Vg-VoltageTarget)/Rwinding);
hold all
xlabel('RPM');
ylabel('I')

linkaxes([ah1 ah2 ah3 ah4 ah5],'x');

title_str = [name{mot_num} '  ' ...
             'VoltageTarget = ' num2str(VoltageTarget) ' V' '   ' ...
             'Kv = '  num2str(Kv(mot_num)) ' VperkRPM' '   ' ...
             'Rw = '  num2str(Rwinding) ' Ohms' '   ' ...
             'Rfets = '  num2str(Rfets) ' Ohms' '   ' ...
             ];
         
ah1 = subplot(5,1,1)
xlabel('RPM');
ylabel('Specific Power (W/RPM)')
title(title_str);
legend(legend_str,'Location','NorthEast');


ah2 = subplot(5,1,2)
xlabel('RPM');
ylabel('Required Boost Factor')
ylim([0 50]);



ah3 = subplot(5,1,3)
xlabel('RPM');
ylabel('Required Boost Duty Cycle')

ah4 = subplot(5,1,4)
xlabel('RPM');
ylabel('Efficiency')

ah5 = subplot(5,1,5)
grid on
xlabel('RPM');
xlabel('RPM');
ylabel(['(N*Kv - ' num2str(VoltageTarget) ')/Rwinding'])
ylim([0 600]);
print('-dpdf',[title_str '.pdf']);

end
end
