%% Dual Phase Margin Analysis: H(s) and Loop Gain

% --- Parameters ---
Rin  = 10e9;        % Input resistor
Cin  = 1.2e-12;     % Input capacitor
Rfb  = 10e6;        % Feedback resistor
rout = 5e6;         % Output resistance
A0   = 1e5;         % Op-amp DC gain
wp   = 2*pi*10e6;   % Op-amp pole (rad/s)
s    = tf('s');

% --- Frequency sweep of Cfb ---
Cfb_vals = logspace(-14, -11, 100);  % Sweep 10 fF to 10 pF
PM_H = nan(size(Cfb_vals));          % Phase margin via H(s)
PM_LG = nan(size(Cfb_vals));         % Phase margin via loop gain

for k = 1:length(Cfb_vals)
    Cfb = Cfb_vals(k);

    try
        % Impedances
        Zin = 1 / (1/Rin + s*Cin);   % Rin || 1/(s*Cin)
        Zfb = 1 / (1/Rfb + s*Cfb);   % Rfb || 1/(s*Cfb)
        Zl  = 1 / (1/rout + s*Cl);   % rout || 1/(s*Cl)
        Ztotal = Zl + Zin + Zfb;

        % Amplifier model
        A = A0 / (1 + s/wp);
        Ab = -A;

        % Coefficients
        ain = rout / Ztot;
        a1  = (Zin + Zfb) / Ztot;
        bin = (Zfb + rout) / Ztot;
        b1  = Zin / Ztot;

        % --- Closed-loop TF ---
        H = ain + (bin * Ab * a1) / (1 - b1 * Ab);

        % --- Open-loop gain ---
        LoopGain = b1 * Ab;

        % --- Phase margins ---
        [~, PM_H(k)] = margin(H);
        [~, PM_LG(k)] = margin(LoopGain);

    catch
        warning('Unstable at Cfb = %.2e F', Cfb);
    end
end

% --- Plotting ---
figure;
semilogx(Cfb_vals*1e12, PM_H, 'b', 'LineWidth', 2); hold on;
semilogx(Cfb_vals*1e12, PM_LG, 'r--', 'LineWidth', 2);
yline(45, 'k:', '45° Target');
xlabel('C_{fb} (pF)');
ylabel('Phase Margin (°)');
title('Phase Margin vs C_{fb}');
legend('Closed-loop H(s)', 'Open-loop -b_1A_b', 'Target 45°');
grid on;
