% Demonstration of sinusoidal steady-state power.
%
% ELEC 226, Rich Kozick

% Magnitude of voltage and current

Vm = 4; % V
Im = 3; % A

% Phase of voltage and current, in degrees

PhiVd = 30;
PhiId = 0;

PhiV = PhiVd/180*pi;  % Convert to radians
PhiI = PhiId/180*pi;

% Generate voltage and current waveforms, with frequency 1 Hz

T = 5;       % Time interval, in seconds
dt = 0.01;   % Spacing between time samples
t = 0:dt:T;  % Sampling times

v = Vm*cos(2*pi*t + PhiV);
i = Im*cos(2*pi*t + PhiI);
vec_ones = ones(size(v));   % Vector of 1's, same size as v and i
p = v .* i;     % Instantaneous power
Pav = Vm * Im / 2 * cos(PhiV - PhiI) * vec_ones;    % Average power
Q = Vm * Im / 2 * sin(PhiV - PhiI) * vec_ones;      % Reactive power
aS = sqrt(Pav.^2 + Q.^2);                             % Apparent power
fprintf('P_{av} = %.2f (W), Q = %.2f (VAR), |S| = %.2f (VA)\n', Pav(1), Q(1), aS(1));

figure(1)
plot(t, v, 'r-', t, i, 'g-', t, p, 'b:', t, Pav, 'b-', t, Q, 'c-', ...
    t, aS, 'k-', t, Pav+aS, 'm:', t, Pav-aS, 'm:')
legend('v(t)', 'i(t)', 'p(t) = v(t) * i(t)', 'P_{av}', 'Q', '|S| = ( P_{av}^2 + Q^2 )^{1/2}', 'P_{av} +/- |S|')
xlabel('TIME (sec)')
ylabel('V / A / W / VARs / VA')
title('Voltage, Current, and Powers')
axis([0, T, min([-Vm,-Im,Pav(1)-aS(1)])-1, Pav(1)+aS(1)+1]);