fmin=1e3;fmax=1e7;
f=10.^(log10(fmin):0.002:log10(fmax));
f1=10e3;
f2=100e3;
f3=1e6;

beta=
bLog=20*log10(beta);
AvLFdB=20*log10(1/beta);

Av0=10000.*1./(1+j.*f/f1).*1./(1+j.*f/f2).*1./(1+j.*f/f3);
Av0dB=20.*log10(abs(Av0));
Av=Av0./(1+beta.*Av0);
AvdB=20.*log10(abs(Av));
arg=360/2/pi.*angle(Av0);
k=find(arg>0);
arg(k)=arg(k)-360;
k=find(abs(bLog+Av0dB)==min(abs(bLog+Av0dB)));
fMax=f(k)
phiM=arg(k)+180

%Dominating -1BP
R1=
C=
Hd1bp=1./(1+j*2*pi.*f*R1*C);

Av0d1bp=Av0.*Hd1bp;
Av0d1bpdB=20.*log10(abs(Av0d1bp));
Avd1bp=Av0d1bp./(1+beta.*Av0d1bp);
Avd1bpdB=20.*log10(abs(Avd1bp));
argd1bp=360/2/pi.*angle(Av0d1bp);
k=find(argd1bp>0);
argd1bp(k)=argd1bp(k)-360;
k=find(abs(bLog+Av0d1bpdB)==min(abs(bLog+Av0d1bpdB)));
phiMd1bp=argd1bp(k)+180;
k=find(abs(Avd1bpdB-AvLFdB-3)==min(abs(Avd1bpdB-AvLFdB-3)));
fMaxd1bp=f(k);
disp(['Dom. -1 BP: R1=',num2str(R1*1e-3,2),'kOhm, C=',num2str(C*1e9,2),'nF'])
disp(['Phase margin=',num2str(phiMd1bp,2),'deg, fMax=',num2str(fMaxd1bp*1e-3,3),'kHz']);

%Lag Compensation
R1=
R2=
C=
Hlag=(R2+1./j/2/pi./f/C)./(R1+R2+1./j/2/pi./f/C);

Av0lag=Av0.*Hlag;
Av0lagdB=20.*log10(abs(Av0lag));
Avlag=Av0lag./(1+beta.*Av0lag);
AvlagdB=20.*log10(abs(Avlag));
arglag=360/2/pi.*angle(Av0lag);
k=find(arglag>0);
arglag(k)=arglag(k)-360;
k=find(abs(bLog+Av0lagdB)==min(abs(bLog+Av0lagdB)));
phiMlag=arglag(k)+180;
k=find(abs(AvlagdB-AvLFdB-3)==min(abs(AvlagdB-AvLFdB-3)));
fMaxlag=f(k);
disp(['Lag comp: R1=',num2str(R1*1e-3,2),'kOhm, R2=',num2str(R2*1e-3,2),'kOhm, C=',num2str(C*1e9,2),'nF'])
disp(['Phase margin=',num2str(phiMlag,2),'deg, fMax=',num2str(fMaxlag*1e-3,3),'kHz']);

%Lead Compensation
R1=
R2=
C=
Hlead=R2./(R2+((R1/j/2/pi./f/C)./(R1+1/j/2/pi./f/C)));

Av0lead=Av0.*Hlead;
Av0leaddB=20.*log10(abs(Av0lead));
Avlead=Av0lead./(1+beta.*Av0lead);
AvleaddB=20.*log10(abs(Avlead));
arglead=360/2/pi.*angle(Av0lead);
k=find(arglead>0);
arglead(k)=arglead(k)-360;
k=find(abs(bLog+Av0leaddB)==min(abs(bLog+Av0leaddB)));
phiMlead=arglead(k)+180;
k=find(abs(AvleaddB-AvLFdB-3)==min(abs(AvleaddB-AvLFdB-3)));
fMaxlead=f(k);
disp(['Lead comp: R1=',num2str(R1*1e-3,2),'kOhm, R2=',num2str(R2*1e-3,2),'kOhm, C=',num2str(C*1e9,2),'nF'])
disp(['Phase margin=',num2str(phiMlead,2),'deg, fMax=',num2str(fMaxlead*1e-3,3),'kHz']);


figure(1)
title('Open loop gain')
semilogx(f,Av0dB,f,Av0d1bpdB,f,Av0lagdB,f,Av0leaddB)
ax=axis;
axis([fmin fmax 0 ax(4)])
legend('Uncomp','Dominating -1BP','Lag','Lead')
xlabel('f [Hz]')
ylabel('\beta\cdotA_V_0 [dB]')

figure(2)
title('\beta\cdot\A_V_0')
clf
subplot(211)
semilogx(f,bLog+Av0dB,f,bLog+Av0d1bpdB,f,bLog+Av0lagdB,f,bLog+Av0leaddB)
ax=axis;
axis([fmin fmax 0 ceil(max(bLog+Av0dB+10)/10)*10])
legend('Uncomp','Dominating -1BP','Lag','Lead')
xlabel('f [Hz]')
ylabel('\beta\cdotA_V_0 [dB]')

subplot(212)
semilogx(f,arg,f,argd1bp,f,arglag,f,arglead)
set(gca,'Ytick',[-270:45:0])
hold on
line([fmin fmax],[-135 -135])
axis([fmin fmax -270 0])
axis([fmin fmax -180 0])
legend('Uncomp','Dominating -1BP','Lag','Lead')
xlabel('f [Hz]')
ylabel('\phi [\circ]')

figure(3)
semilogx(f,AvdB,f,Avd1bpdB,f,AvlagdB,f,AvleaddB)
legend('Uncomp','Dominating -1BP','Lag','Lead',3)
title('Frequency charateristics')
axis([fmin fmax AvLFdB-30 AvLFdB+10])
xlabel('f [Hz]')
ylabel('A_V [dB]')


