%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  sinex.m - Plots voltage-flux-current relationships for a
%            magnetic circuit that exhibits saturation. The
%            flux-current plot is formed by piecewise linear
%            approximation. Normalized flux and current 
%            responses formed for sinusoidal voltage input.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
% Set saturation break point for magnetization curve
satphi=input('Percentage maximum flux = ')/100;
sati=input('Percentage maximum current = ')/100;
% Build piecewise linear flux-current arrays
npts=200; n1=fix(npts/2*satphi); n2=npts-n1;
y1=[[linspace(0,satphi,n1)] [linspace(satphi,1,n2+1)]]; y1(n1)=[];
x1=[[linspace(0,sati,n1)] [linspace(sati,1,n2+1)]]; x1(n1)=[];
nz1=fix(0.2*n1); nz2=ceil(0.2*n2);       % Champher corner
ya=y1(n1-nz1); yb=y1(n1+nz2); xa=x1(n1-nz1); xb=x1(n1+nz2);
m=(ya-yb)/(xa-xb); b=(xa*yb-xb*ya)/(xa-xb);
for i=n1-nz1:n1+nz2; y1(i)=m*x1(i)+b; end
curr=[-fliplr(x1) x1]; curr(npts)=[];    % Completed arrays
phi=[-fliplr(y1) y1]; phi(npts)=[];
plot(curr,phi); grid
title('Magnetization curve');
xlabel('Current');ylabel('Flux');pause
% Build response for the case of an impressed sinsuiodal voltage. 
np=1024; ang=linspace(0,2*pi,1024)'; deg=180/pi;
V=max(phi)*cos(ang);          % Voltage time array
Flx=max(phi)*sin(ang);        % Flux time array
I=interp1(phi, curr, Flx);   % Current time array
figure(2);
subplot(2,1,1); plot(ang*deg,I,ang*deg,Flx,'--',ang*deg,V,'-.');
grid; title('Response due to impressed sinusoidal voltage');
ylabel('Current, voltage, flux (normalized)'); xlabel('Angle');
legend('Coil current','Coil flux','Terminal voltage', 0);
fftI=fft(I); f=abs(fftI); f(1)=f(1)/2; f=f/(np/2);  % FFT of current
pause
subplot(2,1,2); plot([0:25],f(1:26)); grid
title('Current for sinusoidal voltage');
ylabel('Current'); xlabel('Harmonic number');