% Example 10.7

clc
clear

% system parameters
R = 12.5e-3;
L = 1.2;
V0 = 15;
omega = 120*pi;

% part (a)

% dc current
Idc = 35;

% dc voltage
Vdc_a = R*Idc;

% Power
P = Vdc_a*Idc;

%Calculate the delay angle
alpha_da = acos(pi*R*Idc/(2*V0));

%Now calculate the load voltage
for n = 1:1300
   theta(n) = 2*pi*(n-1)/1000;
   t(n) = theta(n)/omega; 
   vs(n) = V0*sin(theta(n));

   if theta(n) < alpha_da
	vL(n) = -vs(n);
   elseif (theta(n) < pi + alpha_da)
	vL(n) = vs(n);
   elseif theta(n) < 2*pi + alpha_da
	vL(n) = -vs(n);
   elseif theta(n) < 3*pi + alpha_da
	vL(n) = vs(n);
   elseif theta(n) < 4*pi + alpha_da
	vL(n) = -vs(n);
   else 
	vL(n) = vs(n);
   end

end

plot(1000*t,vL)
xlabel('time [msec]')
ylabel('Load voltage [V]')

pause

% part (b)

% delay angle
alpha_db = 0.9*pi;

% Find the new dc voltage
Vdc_b = (2*V0/pi)*cos(alpha_db);

% Time constants
tau = L/R;

% Initial current
im0 = Idc;

% Calculate the time at which the current reaches zero
tzero = -tau*log((-Vdc_b/R)/(im0-Vdc_b/R));

% Now plot the load voltage
for n = 1:1300
   theta(n) = 2*pi*(n-1)/1000;
   t(n) = theta(n)/omega; 
   vs(n) = V0*sin(theta(n));

   if theta(n) < alpha_db
	vL(n) = -vs(n);
   elseif (theta(n) < pi + alpha_db)
	vL(n) = vs(n);
   elseif theta(n) < 2*pi + alpha_db
	vL(n) = -vs(n);
   elseif theta(n) < 3*pi + alpha_db
	vL(n) = vs(n);
   elseif theta(n) < 4*pi + alpha_db
	vL(n) = -vs(n);
   else 
	vL(n) = vs(n);
   end

end

plot(1000*t,vL)
xlabel('time [msec]')
ylabel('Load voltage [V]')

fprintf('part (a):')
fprintf('\n  Vdc_a = %g [mV]',1000*Vdc_a)
fprintf('\n  Power = %g [W]',P);
fprintf('\n  alpha_d = %g [rad] = %g [degrees]',alpha_da,180*alpha_da/pi)
fprintf('\npart (b):')
fprintf('\n  alpha_d = %g [rad] = %g [degrees]',alpha_db,180*alpha_db/pi)
fprintf('\n  Vdc_b = %g [V]',Vdc_b)
fprintf('\n  Current will reach zero at %g [sec]',tzero)
fprintf('\n')