clear all; close all;

%% Werte der Bauteile und Konstanten
L   = 4.7e-6;
Rl  = 5.0e-3;
Co  = 47.0e-6;
Rc  = 5.0e-3;
Vi0 = 10.0;
Ro  = 8.6;
D0  = 0.1674;

%% Berechnung
f_vect = logspace(1, 8, 1000);
vout_del_d_vect = [];
da = 1.0e-6; % Abweichung für numerische Ableitung
for f = f_vect;
  s = 1j*2.0*pi*f;
  
  % positive Abweichung
  dd = +da;
  K = (1-(D0+dd))/(s*L+Rl);
  M = (Rc*Ro)/(Rc+Ro);
  Q = ((1-(D0+dd))^2)/(s*L+Rl);
  T = (Rc*Ro)/(Rc+Ro);
  H = 1/(1+s*Rc*Co);
  W = Ro/(Rc+Ro);
  vout_pos = K * M * ((1 + Q*T - H*W)^(-1)) * Vi0;
  
  % negative Abweichung
  dd = -da;
  K = (1-(D0+dd))/(s*L+Rl);
  M = (Rc*Ro)/(Rc+Ro);
  Q = ((1-(D0+dd))^2)/(s*L+Rl);
  T = (Rc*Ro)/(Rc+Ro);
  H = 1/(1+s*Rc*Co);
  W = Ro/(Rc+Ro);
  vout_neg = K * M * ((1 + Q*T - H*W)^(-1)) * Vi0;
  
  vout_del_d = (vout_pos - vout_neg) / (2.0 * da); % numerische Ableitung berechnen
  vout_del_d_vect = vertcat(vout_del_d_vect, vout_del_d);
end;

%% Plot
figure;
hold on;
[ax, hline1, hline2] = plotyy(f_vect, 20*log10(abs(vout_del_d_vect)), f_vect, (180/pi)*unwrap(angle(vout_del_d_vect)), @semilogx, @semilogx);
set(hline1, 'Color', 'red');
set(hline2, 'Color', 'blue');
set(ax(1), 'ycolor', 'red');
set(ax(2), 'ycolor', 'blue');
ylabel(ax(1), "Gain [dB]");
ylabel(ax(2), "Phase [deg]");
xlabel("Frequency [Hz]");
grid on;