%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   tranmag.m - Calculates magnetization curve for core or
%               shell single-phase transformer.
%               Calls sptdata.m type, winding, turns and
%               dimension data.
%               Joint mmf drop calculated by joint.m.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
sptdata;
joint;
if style == 'shell_design';
   A=d*l/2*SF; kphi=2;
   MLP1=2*h+w+3*l/2; MLP2=w+l/2;
elseif style == 'core__design';
   A=d*l*SF; kphi=1;
   MLP1=2*h+3*l+w; MLP2=w+l;
else; end
npts=200; Im=zeros(1,npts+1); V=Im; B=linspace(0,Bmax,npts);
BR=VR/(4.44e-08*N*A*fR*kphi); B=[B';BR];
for i=1:npts+1
   FJ=interp1(BT, FT, B(i)); % Look up joint mmf
   Im(i)=(Hm7(B(i))*MLP1+Hm7cg(B(i))*MLP2+2*FJ)/N;
   V(i)=4.44e-08*N*fR*B(i)*A*kphi;
end
Xm=0.99*VR/Im(npts+1)*sqrt(2);
if style=='shell_design'
   Pcore=2*Pcm7(BR)*2*0.283*(l+h+w)*l*d*SF; % 2 x Epstein losses
elseif style=='core__design'
   Pcore=2.5*Pcm7(BR)*2*0.283*(2*l+w+h)*l*d*SF;
else; end
Rc=(0.99*VR)^2/Pcore;
disp(' '); disp('   MAGNETIZING REACTANCE & CORE LOSS RESISTANCE');
disp(' '); disp([blanks(7) 'Xm (ohms)' blanks(10) 'Rc (ohms)']);
disp([blanks(6) num2str(Xm) blanks(9) num2str(Rc)]);
pause; clf; figure(1);


plot(Im(1,1:npts),V(1,1:npts),Im(1,npts+1),V(1,npts+1),'o');
grid; title('Transformer magnetization curve');
xlabel('Peak value magnetizing current( I_m ), A');
ylabel('Primary voltage( V_1 ), V');
text(0.50*max(Im), 0.75*max(V),design);
legend('Saturation curve','Rated voltage',4);