a0 = 1e5;
w1 = 1e4;
w2 = 1e6;
s = tf('s');
a = a0/(1+s/w1)/(1+s/w2)
h = bode(a,'r');
a_norm = a / dcgain(a);
step(a_norm,'r')
title('Normalized Open-Loop Step Response');
ylabel('Normalized Amplitude');
A  =  1 / b
A0 = 10;
b = 1 / A0;    % approximation for ab>>1
R1 = 10000;
R2 = R1 * (1/b - 1)
A = feedback(a,b);
bodemag(a,'r',A,'b');
legend('Open-Loop Gain (a)','Closed-Loop Gain (A)')
ylim([0,110]);
L = a * b;
S = 1 / (1 + L);
S = feedback(1,L);
bodemag(A,'b',S,'g')
legend('Closed-Loop Gain(A)', 'System Sensitivity(S)',4)
step(A)
margin(L)
[Gm,Pm,Wcg,Wcp] = margin(L);
C = 1/(R2*Wcp)

###########################################
## this causes an error
K = R1/(R1+R2);
C = [1:.2:3]*1e-12;
for n = 1:length(C)
    b_array(:,:,n) = tf([K*R2*C(n) K],[K*R2*C(n) 1]);
end
###########################################

A_array = feedback(a,b_array);
L_array = a*b_array;

step(A,'b:',A_array,'b',[0:.005:1]*1.5e-6);
title('Closed-Loop Step Response (Compensated)');

[Gm,Pm,Wcg,Wcp] = margin(L_array);
plot(C*1e12,Pm,'g');
ax = gca;
xlim([0.8 3.6]);
ylim([45 60]);
ax.Box = 'on';
xlabel('Compensation Capacitor, C (pF)');
ylabel('Phase Margin (deg)')

A_comp = A_array(:,:,6);
step(A,'b:',A_comp,'b')
legend('Uncompensated (0 pF)','Compensated (2 pF)')
bode(a,'r',A,'b:',A_comp,'b')
legend('a(s)','A(s)','A\_comp(s)');
