%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  dcOCsat.m - plots the magnetization curve for a dc
%              machine. Calls dcdata.m for input data.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
dcdata,
kE=1.666667e-10*p*N*2*nc/a; kT=1.173e-9*p*N*2*nc/a;
phipR=TdR/kT/IaR;
% Build apparent tooth flux density array
Bt=linspace(0,180000,100);
for i=1:length(Bt); Ht(i)=hm22(Bt(i)); end
lam3=pi*(d-4/3*ds)/N;
kt=3.2*(SF*lam3/(lam3-bs)-1);
for i=1:length(Ht)
   Ba(i)=Bt(i)+kt*interp1(Bt,Ht,Bt(i));
end
taut=fix(psi*N/p)+0.5;        % Teeth per pole span
At3=SF*la*(lam3-bs)*taut;     % Total tooth area @ 1/3 depth
rc=(d-2*ds-dshft)/2;
Ac=2*SF*la*rc;                % Rotor core area
Af=wf*tf;                     % Frame area
taup=taut*p/N*pi/4*(d+2*del); % Pole arc
Ask=SF*(la+0.125)*wsk;        % Pole shank area
Ash3=SF*(la+0.125)*(2*wsk+taup)/3; % Shoe area @ 1/3 depth
lam=pi*d/N;                   % Slot pitch
qty=lam*(5*del+bs); ks=qty/(qty-bs^2); % Carter coefficient
phip=linspace(0,1.2*phipR,200); 
phip=[phip phipR]; m=length(phip);
for i=1:m
   ATt=interp1(Ba,Ht,phip(i)/At3)*ds;
   ATc=hm22(phip(i)/Ac)*pi/2/p*(d-ds-rc/2);
   ATf=h1010(LF*phip(i)/Af/2)*pi/2/4*(ODf-tf);
   ATsh=hm22(LF*phip(i)/Ash3)*lsh;
   ATsk=hm22(LF*phip(i)/Ask)*(ODf-d-2*del-2*lsh-2*tf)/2;
   ATg=p*ks*del*phip(i)/pi/0.665/d/la/3.2;
   ATp(i)=ATg+ATsk+ATsh+ATf+ATc+ATt;
end
plot(ATp(1:m-1),phip(1:m-1)/1000,ATp(m),phip(m)/1000,'o');
title('Magnetization curve for dc machine'); grid
xlabel('MMF per pole, A-t');
ylabel('Flux per pole, kilolines');
% Generate & save OCC for use in performance study
EIfdes=[kE*nmR*phip(1:m-1)' ATp(1:m-1)'/Nf];
save EIfdes EIfdes -ascii -double