%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  hyst.m - Plots voltage-flux-current relationships for a coil
%           wound on a ferromagnetic structure that exhibits
%           hysteresis. The hysteresis loop is formed by piece-
%           wise linear approximation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; 
y1=[-1 -0.9 0 0.9 1]'; x1=[-1 -0.3 0 0.3 1]'; % Base curve for
Icmax=0.1;    %Maximum coercive current       % hysteresis buildup
pctsat=input('Percentage maximum flux = ')/100; Ic=pctsat*Icmax;
if pctsat > 1; pctsat=1; else; end
if pctsat > y1(4); x1max=interp1(y1,x1,pctsat);
   if x1max <= x1(4)+Ic/2; x1max=1.01*(x1(4)+Ic/2); else; end
   x1(1)=-x1max; x1(5)=x1max; y1(1)=-pctsat; y1(5)=pctsat;
   xu=x1+[0 Ic/2 Ic Ic/2 0]'; xd=x1-[0 Ic/2 Ic Ic/2 0]'; yu=y1; yd=y1;
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);
   ydwn=pctsat^2/y1(4);
   if abs(ydwn-pctsat)<1e-4; ydwn=0.85*pctsat; end
   xdwn=interp1(yd,xd,ydwn);
   xu=[-xmaxu -xdwn Ic xmaxu]'; xd=-xu;
   yu=[-pctsat -ydwn 0 pctsat]'; yd=-yu;
end
figure(1);          % Plot hysteresis loop
plot(xu,yu,xd,yd); grid; title('Hysteresis loop');
xlabel('Normalized current'); ylabel('Normalized flux');  pause; clf; 
% Build response for the case of an injected sinusoidal current.
ang=linspace(0,2*pi,1024)'; x=max(xu)*sin(ang); npts=length(x);
y=zeros(npts,1); y(1)=interp1(xu,yu,x(1));
for i=2:npts
   if x(i) > x(i-1)
      y(i)=interp1(xu,yu,x(i));
   else y(i)=interp1(xd,yd,x(i));
   end;
end
% Frequency domain differentiation of flux to yield voltage
dy=fft(y); deg=180/pi;
for i=1:1024; dy(i)=-j*i*dy(i); end; 
dy=ifft(dy); dy=real(dy*max(y)/max(abs(dy)));
plot(ang*deg,x,ang*deg,y,'--',ang*deg,dy,'-.'); grid;
title('Response due to injected sinusoidal current');
ylabel('Current, voltage, flux (normalized)'); xlabel('Degrees');
legend('Coil current','Coil flux','Terminal voltage'); pause
legend off; f=abs(fft(dy)); f(1)=f(1)/2; f=f/max(f)*100;
plot([0:25],f(1:26)); grid
title('Voltage for sinusoidal current');
ylabel('Percent'); xlabel('Harmonic number');pause;
% Build response for the case of an impressed sinusoidal voltage. 
ang=linspace(0,2*pi,1024)'; y=max(yu)*sin(ang); npts=length(y);
x=zeros(npts,1); x(1)=interp1(yu,xu,y(1));
for i=2:npts
   if y(i) > y(i-1)
      x(i)=interp1(yu,xu,y(i));
   else x(i)=interp1(yd,xd,y(i));
   end
end
% Frequency domain differentiation of flux to yield voltage
dy=fft(y); deg=180/pi;
for i=1:1024; dy(i)=-j*i*dy(i); end; 
dy=ifft(dy); dy=real(dy*max(y)/max(abs(dy)));
% Determine fundamental component of current
x1=fft(x); x1=abs(x1(2))*cos(ang+angle(x1(2)))/512;
plot(ang*deg,x,ang*deg,y,'--',ang*deg,dy,'-.',ang*deg,x1,':');
grid; title('Response due to impressed sinusoidal voltage');
ylabel('Current, voltage, flux (normalized)'); xlabel('Angle');
legend('Coil current','Coil flux','Terminal voltage', ...
   ' Fund. current'); pause;
legend off; f=abs(fft(x)); f(1)=f(1)/2; f=f/max(f)*100;
plot([0:25],f(1:26)); grid
title('Current for sinusoidal voltage');
ylabel('Percent'); xlabel('Harmonic number');
legend off;
Im=abs(fft(x)); save Imag Im  % Im saved for use with imthd.m