%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  hyloop.m - Builds embedded hysteresis loops & calculates areas. 
%             Techniques of this program used by hyst.m.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clf; 
Icmax=0.2;  % Maximum coercive current. Change to alter loop widths.
puFsat=0.7; % Per unit flux at saturation curve knee.
puIsat=0.3; % Per unit current at saturation curve knee.
nloops=5; pctsat=linspace(0,1,nloops+1); % No. loops & peak fluxes
for k=1:nloops+1
   y1=[-1 -puFsat 0 puFsat 1]'; x1=[-1 -puIsat 0 puIsat 1]';
   Ic=pctsat(k)*Icmax;
   if pctsat(k) > y1(4); x1max=interp1(y1,x1,pctsat(k));
      if x1max <= x1(4)+Ic/2; x1max=1.01*(x1(4)+Ic/2); else; end
   x1(1)=-x1max; x1(5)=x1max; y1(1)=-pctsat(k); y1(5)=pctsat(k);
   xu=x1+[0 Ic/2 Ic Ic/2 0]'; xd=x1-[0 Ic/2 Ic Ic/2 0]'; yu=y1; yd=y1;
   wu=0; wd=0;
      for i=1:4  % Calculate hysteresis loop area
      wu=wu+(xu(i)+xu(i+1))/2*(yu(i+1)-yu(i));
      wd=wd+(xd(i)+xd(i+1))/2*(yd(i+1)-yd(i));
      end
   W(k)=wu-wd;
   else
   xu=x1+[0 Ic/2 Ic Ic/2 0]'; xd=x1-[0 Ic/2 Ic Ic/2 0]'; yu=y1; yd=y1;
   xmaxu=interp1(yu, xu, pctsat(k));
   ydwn=pctsat(k)^2/y1(4);
   if abs(ydwn-pctsat(k))<1e-4; ydwn=0.85*pctsat(k); end
   xdwn=interp1(yd, xd, ydwn);
   xu=[-xmaxu -xdwn Ic xmaxu]'; xd=-flipud(xu);
   yu=[-pctsat(k) -ydwn 0 pctsat(k)]'; yd=-flipud(yu);
   wu=0; wd=0;
   for i=1:3     % Calculate hysteresis loop area
      wu=wu+(xu(i)+xu(i+1))/2*(yu(i+1)-yu(i));
      wd=wd+(xd(i)+xd(i+1))/2*(yd(i+1)-yd(i));
      end
   W(k)=wu-wd;
   end
% Plot hysteresis loop
plot(xu,yu,xd,yd); hold on;
end
grid; title('Hysteresis loops');
xlabel('Normalized current'); ylabel('Normalized flux');
hold off
disp(' '); disp(' ');
disp('   Per Unit Energy Loss of Hysteresis Loops (min -> max)');
Energy=[W(2:nloops+1)]; disp(' '); disp(Energy);