%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  imeval.m  -  Calculates induction motor performance
%               given winding & dimensional data
%               Calls imdesign.m for input data
%               Calls Pc.m during Rc calculation
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
imdesign,% imdesign.m contains dimension, winding, & rating data
harris,
V1=VL/sqrt(3)+0*j; ws=2/p*2*pi*f;
Ns1=S1/p/m; Cphi=Ns1*Cs*p/a; del=(D-Dr)/2; 
lam1=pi*D/S1; lam2=pi*Dr/S2;
% Fringing coefficients
qty=lam1*(5*del+bs1); ks=qty/(qty-bs1^2);
qty=lam2*(4.4*del+0.75*br2); kr=qty/(qty-br2^2);

% Stator winding factor
gam=p*pi/S1; kd=sin(Ns1*gam/2)/Ns1/sin(gam/2);
rho=pitch*p/S1*pi; kp=cos(pi/2-rho/2); kw1=kd*kp;

% Flux density checks
phim=0.97e08*a*V1/2.22/kd/kp/p/Ns1/f/Cs;
ts1=lam1-bs1; Bstm=pi*p*phim/2/S1/ts1/la/SF;
Bscm=phim/(Do-D-2*ds1)/la/SF;
tr2=pi*(Dr-2*dr1)/S2-br1; Brtm=pi*p*phim/2/S2/tr2/la/SF;
Brcm=phim/(Dr-Dsh-2*dc-2*dr1)/la/SF;
clc; disp(' '); disp(' '); 
disp('FLUX DENSITY CHECKS - kilolines/sq. in'); disp(' ');
disp(['Bstm = ', num2str(Bstm/1000)]);
disp(['Bscm = ', num2str(Bscm/1000)]);
disp(['Brtm = ', num2str(Brtm/1000)]);
disp(['Brcm = ', num2str(Brcm/1000)]);
% R1 calculation
alf=asin((bs1+0.050)/lam1); tau=pi*D/p;
lalf=rho*tau/(2*pi*cos(alf)); lc=la+1+2*(lalf+ds1);
R1=1.056e-06*lc*Cphi/a/sa1; R1=1.05*R1;  % ac value

% Xm calculation
dele=ks*kr*del;
Xm=1.6713e-07*f*D*la/p/dele*(kw1*S1*Cs/2/a/m)^2;

% Rc calculation

Bt13=pi*p*phim/2/(pi*(D+2*ds1/3)-S1*bs1)/la/SF;
Bscc=0.866*phim/(Do-D-2*ds1)/la/SF;
Wt=(pi/4*((D+2*ds1)^2-D^2)-S1*bs1*ds1)*la*SF*0.283;
Pt=Wt*Pc(Bt13)*4;
Wy=pi/4*(Do^2-(D+2*ds1)^2)*la*SF*0.283;
Py=Wy*Pc(Bscc)*4;
Rc=3*(0.97*VL/sqrt(3))^2/(Pt+Py);

%  X1 components calculation
thetask=p*pi/S1;
Xsk=Xm*thetask^2/12;
pupitch=p*pitch/S1;       % Phase factor
if pupitch <= 0.5; Ks=pupitch;
elseif pupitch <= 2/3; Ks=-0.25+1.5*pupitch;
else; Ks=0.25+0.75*pupitch; end
Fs=(ds1-ds0)/3/bs1+ds0/bs1;
Xss=2.0055e-07*m*f*Cphi^2*la*Ks*Fs/S1;
if S2<S1; ks1=1; else; ks1=S1/S2; end
a1=abs(pi*D/S1-bs1-br2)/2;
Fz=a1/4/del+dr2*(1-ks1/2)/br2;
Xz1=2.0055e-07*m*f*Cphi^2*la*Ks*Fz/S1;
Fes=1+4*log(3.94*D/S1/sqrt(bs1*ds1));
Xes=0.3192e-07*f*(kw1*Cphi)^2*(0.5+lalf)*Ks*Fes/p;

%  X2 components calculation
if S1<S2; kr1=1; else; kr1=S2/S1; end
a2=abs(pi*Dr/S2-bs1-br2)/2;
q2=S2/m/p; gamr=p*pi/S2; kdr=sin(q2*gamr/2)/q2/sin(gamr/2);
rhor=fix(S2/p)*p*pi/S2; kpr=cos((pi-rhor)/2); kwr=kdr*kpr;
Fr=(dr1-dr2)/3/br1+dr2/br2;
Xrs=2.0055e-07*m*f*(kw1*Cphi)^2*la*Fr/kwr^2/S2;
Fz2=a2/4/del+ds0*(1-kr1/2)/bs1;
Xz2=2.0055e-07*m*f*(kw1*Cphi)^2*la*Fz2/kwr^2/S2;
fer1=S2*Der/m/p*(1+4*log(3.94*Dr/S2/sqrt(wer*der)));
fer2=4*br*p*(1+4*log(3.94*Dr/S2/sqrt(br1*dr1)));
Xer=7.98e-09*f*(kw1*Cphi)^2/p^2*(fer1+fer2);

I1m=zeros(1,100); PF=I1m; TTd=I1m; eta=I1m; Ps=I1m;
wm=linspace(0,ws-0.01,100); npts=length(wm);
for i=1:npts;      %  Speed loop - performance calculations
s=(ws-wm(i))/ws; fr=s*f;
if krc==1; kalf=0.3505; kbet=1; else; kalf=0.2700; kbet=1.64; end
alfd=kalf*sqrt(fr*b2c/br1)*d2c; 
K1=alfd*(sinh(2*alfd)+sin(2*alfd))/(cosh(2*alfd)-cos(2*alfd));
J1=3*(sinh(2*alfd)-sin(2*alfd))/(2*alfd*(cosh(2*alfd)-cos(2*alfd)));
%  R2 calculation with skin effect adjustment
R2=kbet*1.169e-06*m*(kw1*Cphi/kwr)^2*(lb*K1/S2/sb+2*Der*Kring/pi/p^2/sr);
X1=Xsk+Xes+Xss+Xz1;
X2=Xer+Xz2+Xrs*J1;      % Adjustment for skin effects
Z1=R1+X1*j; Z2=R2/s+X2*j; ZM=(0+Rc*Xm*j)/(Rc+Xm*j);
Z=Z1*Z2+Z1*ZM+Z2*ZM; I1=V1*(Z2+ZM)/Z; I1m(1,i)=abs(I1);  % Stator current
I2=ZM/(Z2+ZM)*I1; E=I2*Z2;
PF(1,i)=cos(atan2(imag(Z/(Z2+ZM)),real(Z/(Z2+ZM)))); % Power factor
TTd(1,i)=3*abs(I2)^2*R2/s/ws;                        % Developed torque
Pfw=0.03*HP*746*(wm(i)/ws)^(2.8);                    %  F&W losses
Ps(1,i)=TTd(1,i)*wm(i)-Pfw; if Ps(1,i)<0; Ps(1,i)=0; end      
Pin=3*real(V1*conj(I1)); eta(1,i)=Ps(1,i)/Pin*100;   % Efficiency 
end;
nm=wm*30/pi; Ps=Ps/746;
% Plot performance
figure(1);
subplot(2,1,1), plot(nm,Ps); grid; title('Output power');
xlabel('Speed, rpm'); ylabel('Power, hp');
subplot(2,1,2), plot(nm,TTd);grid; title('Developed torque');
xlabel('Speed, rpm'); ylabel('Developed torque, N-m');
figure(2);
subplot(2,1,1), plot(nm,I1m);grid; title('Input current');
xlabel('Speed, rpm'); ylabel('Stator current, A'); 
subplot(2,1,2), plot(nm,PF);grid; title('Input power factor');
xlabel('Speed, rpm'); ylabel('Input PF');
figure(3);
subplot(2,1,1), plot(nm,eta);grid; title('Efficiency');
xlabel('Speed, rpm'); ylabel('Efficiency, %');
% Linear interprolation for rated conditions
for i=1:npts-1; if Ps(npts-i) >= HP; k=npts-i; break; end; end;
if i==npts-1; disp('DESIGN CANNOT PRODUCE RATED OUTPUT');
else;
nmR=nm(k+1)-(HP-Ps(k+1))*(nm(k+1)-nm(k))/(Ps(k)-Ps(k+1));
delnm=nm(k+1)-nm(k); delnm1=nm(k+1)-nmR;
I1R=I1m(k+1)+(I1m(k)-I1m(k+1))*delnm1/delnm;
TTdR=TTd(k+1)+(TTd(k)-TTd(k+1))*delnm1/delnm-0.03*HP*746* ...
(nmR*pi/30)^(1.8)/ws^(2.8);
etaR=eta(k+1)+(eta(k)-eta(k+1))*delnm1/delnm;
PFR=PF(k+1)+(PF(k)-PF(k+1))*delnm1/delnm;
disp(' '); disp(' ');disp('RATED CONDITIONS');disp(' ');
disp(['Speed(rpm) = ',num2str(nmR)]);
disp(['Current(A) = ',num2str(I1R)]);
disp(['Torque(N_m) = ',num2str(TTdR)]);
disp(['Efficiency(%) = ',num2str(etaR)]);
disp(['Power factor = ',num2str(PFR)]);
end