% Power Factor Correction (lag only, with power waves)
% Carlton Cross  28 May 10  Rev 1 Dec 11

clear all
close all

Vline = 220;            % line voltage, Vrms
P = 30000;              % Machine power, W
Pphase = P/3;           % Machine power per phase, W
pftxt = ['lag'];        % annotation text (only lag is possible)
pf = 0.83;               % motor power factor (lag)
pfd = 0.95;              % desired power factor (lag)
Zline = 0.2 + 0.1j;     % line impedance (series)

w = 2*pi*60;            % line frequency (60 Hz)
Vp = Vline/sqrt(3);     % line voltage
Sm = Pphase*(1 + j*tan(acos(pf)));   % motor power
Sd = Pphase*(1 + j*tan(acos(pfd)));  % corrected motor power
Sc = Sd - Sm;                   % capacitor reactive power
C = Sc/(-j*w*Vp^2);
Im = conj(Sm/Vp);
Ic = j*w*C*Vp;
Iline = Im + Ic;
Plossold = 3*real(Zline)*abs(Im)^2
Plossnew = 3*real(Zline)*abs(Iline)^2;
Psave =  Plossold - Plossnew;
Psavepercent = 100*Psave/P;

T = 2*pi/w;
npts = 500;
t = linspace(-0.2*T, 1.2*T, npts);
vp = sqrt(2)*Vp*sin(w*t);                           % motor voltage
im = sqrt(2)*abs(Im)*sin(w*t + angle(Im));          % motor current
ic = sqrt(2)*abs(Ic)*sin(w*t + angle(Ic));          % capacitor current
iline = sqrt(2)*abs(Iline)*sin(w*t + angle(Iline)); % line current
Pmt = vp.*im;                                       % instantaneous power


% Set up figure positions.
scrnz = get(0,'ScreenSize');
scale = scrnz(3)/scrnz(4);
figure('Position',[0.02*scrnz(3) 0.1*scrnz(4) 0.8*scrnz(3) 0.8*scrnz(4)])
figure(1)
plot(t,vp,'k', t,im,'b--', t,ic,'r-.', t,iline,':m', 'linewidth',1.5)
v = axis;   % get the graph edges
xmin = v(1); xmax = v(2); ymin = v(3);  ymax = v(4);
hold on
plot(t,Pmt/500,'c','linewidth',1.5)
plot([xmin xmax],[0 0],'k','linewidth',1.5)   % horizontal axis
plot([0 0],[ymin ymax],'k','linewidth',1.5)   % vertical axis
hold off
grid
set(gca,'xtick',[0 T/4 T/2 3*T/4 T])
set(gca,'ytick',[])
legend('Line Voltage, V','Motor Current, A',...
    'Capacitor Current, A','Line Current, A','Instantaneous Power, VA')
title('\fontsize{14}Waveforms and Power-Factor')
ylabel('Magnitude, V, A, VA')
xlabel('Time, s')
vx = axis;   % axis gets the graph end points
xmin = vx(1); xmax = vx(2); ymin = vx(3); ymax = vx(4);
xrel = 0.31;                    % establish a relative position for text
yrel = 0.3;
vspace = (ymax-ymin)/25;        % spacing for 25 lines in the box
px = xrel*(xmax-xmin) + xmin;   % x coordinate for text
py = yrel*(ymax-ymin) + ymin;   % y coordinate for text
text(px,(py-0*vspace), ['Motor Power Factor = ' num2str(pf) ' ' pftxt],...
    'backgroundcolor',[1 1 1])
text(px,(py-1*vspace), ['Desired Power Factor = ' num2str(pfd) ' ' pftxt],...
    'backgroundcolor',[1 1 1])
text(px,(py-2*vspace), ['Motor Power = ' num2str(round(3*Pphase)) ' W'],...
    'backgroundcolor',[1 1 1])
text(px,(py-3*vspace), ['Line Impedance = ' num2str(Zline) ' \Omega'],...
    'backgroundcolor', [1 1 1])
text(px,(py-4*vspace), ['Old Line Loss = ' num2str(round(Plossold)) ' W'],...
    'backgroundcolor',[1 1 1])
text(px,(py-5*vspace), ['New Line Loss = ' num2str(round(Plossnew)) ' W'],...
    'backgroundcolor',[1 1 1])
text(px,(py-6*vspace), ['Power Saved = ' num2str(round(Psave)) ' W', ...
    '  ' num2str(round(100*Psavepercent)/100) ' %'],    'backgroundcolor',[1 1 1])
text(xmin,(ymin-2.5*vspace),...
    ['Carlton Cross  From file p:\\matlab\\e228\\exmpl\\powerfactordemo.m' ...
    '    plotted ' date]);
