Vm = 100; thetav = 0;         % Voltage amplitude and phase angle
Z = 1.25;  gama = 60;       % Impedance magnitude and phase angle
thetai = thetav - gama;           % Current phase angle in degree
theta = (thetav - thetai)*pi/180;              % Degree to radian
Im = Vm/Z;                                    % Current amplitude
wt=0:.05:2*pi;                                % wt from 0 to 2*pi
v=Vm*cos(wt);                             % Instantaneous voltage
i=Im*cos(wt + thetai*pi/180);             % Instantaneous current
p=v.*i;                                     % Instantaneous power
V=Vm/sqrt(2); I=Im/sqrt(2);             % rms voltage and current
P = V*I*cos(theta);                               % Average power
Q = V*I*sin(theta);                              % Reactive power
S = P + j*Q                                       % Complex power
pr = P*(1 + cos(2*(wt + thetav)));                    % Eq. (2.6)
px = Q*sin(2*(wt + thetav));                          % Eq. (2.8)
PP=P*ones(1, length(wt));  % Average power with length w for plot
xline = zeros(1, length(wt));           % generates a zero vector
wt=180/pi*wt;                       % converting radian to degree
subplot(2,2,1), plot(wt, v, wt, i,wt, xline), grid
title(['v(t)=Vm coswt, i(t)=Im cos(wt +', num2str(thetai), ')'])
xlabel('wt, degree')
subplot(2,2,2), plot(wt, p, wt, xline), grid
title('p(t)=v(t) i(t)'), xlabel('wt, degree')
subplot(2,2,3), plot(wt, pr, wt, PP, wt,xline), grid
title('pr(t)   Eq. 2.6'), xlabel('wt, degree')
subplot(2,2,4), plot(wt, px, wt, xline), grid
title('px(t)  Eq. 2.8'), xlabel('wt, degree')
subplot(111)