% phase_margin_2
% does multiple plots
kvals=[.125 .287 .55 .92 1.5 2.6 5.7 1000];
% following wvals for wp1 always make ft = 1MHz
wvals=2*pi*1E+6*[2.917 2.004 1.544 1.302 1.155 1.063 1.015 1];

for kt=8:-1:1
    kf=kvals(kt);

% set up symbolic variable s for transfer function
s = tf('s');
% frequencies 1Hz to 100MHz
f=logspace(0,8,161);

% Amplifier
wp_1  = wvals(kt);  %       Unity gain freq of op-amp (dominant pole)
% second pole factor of kf higher than w1
wp_2  =  kf*wp_1;

% denominator polynomial for finding roots numerically
den_poly=[(1/(wp_1*wp_2)) (1/wp_1)  1];
% pole locations (in Hz)
pole_loc=(1/(2*pi))*roots(den_poly);

% Op amplifier open loop transfer function
A=1/((s/wp_1)*(1+s/wp_2));

% Magnitude and phase of A
AOL=squeeze(freqresp(A,2*pi*f));  % need to squeeze out extra dimensions
magAOL=abs(AOL);
phaseAOL=(180/pi)*unwrap(angle(AOL));


fT=interp1(magAOL,f,1)

% phase at unity gain frequency
pfT=interp1(f,phaseAOL,fT)

% phase margin
pmargin=pfT+180;




% Closed loop response
G=A/(A+1);

GCL=squeeze(freqresp(G,2*pi*f));  % need to squeeze out extra dimensions
magGCL=abs(GCL);
phaseGCL=(180/pi)*unwrap(angle(GCL));

% open separate figure window for each plot
figure(kt)
% plot
loopgain_plot
% requires keystroke to continue
pause
end  

