clear all;
close all;
clc;

%% Parameters
Adc = 1e4;
wc = 2*pi*20e3;

Rin = 10e6;
Cin = 1.2e-12;
Rfb = 400e3;
Cfb = 100e-15;
rout = 2.5e6;
Cl = 1.2e-12;
Cparas = 5e-12;
Zbuff = 50;


%% Impedances
s = tf('s');

Zin1 = 1/(1/Rin+s*Cin);
Zin = Zin1 + Zbuff;
Zfb = 1/(1/Rfb+s*Cfb);
Zout = rout;
Zl = 1/(s*Cl);
Zparas = 1/(s*Cparas);

Zin = zpk(Zin);
Zfb = zpk(Zfb);
Zout = zpk(Zout);
Zl = zpk(Zl);
Zparas = zpk(Zparas);
Zbuff = zpk(Zbuff);

Za = 1/(1/Zout+1/Zl);
Za = zpk(Za);
Zb = 1/(1/(Za+Zfb)+1/Zparas);
Zb = zpk(Zb);
Zc = 1/(1/Zin+1/Zparas);
Zc = zpk(Zc);
Zd = 1/(1/(Zc+Zfb)+1/Zl);
Zd = zpk(Zd);

%% Open loop transfer function
A = Adc/(1+s/wc);
A = zpk(A);
A = minreal(A);

%% Loop parameters

ain = (Za/(Za+Zfb))*(Zb/(Zb+Zin));
ain = zpk(ain);
ain = minreal(ain);

a1 = Zd/(Zd+Zout);
a1 = zpk(a1);
a1 = minreal(a1);

bin = Zb/(Zb+Zin);
bin = zpk(bin);
bin = minreal(bin);

b1 = (Zd/(Zd+Zout))*(Zc/(Zfb+Zc));
b1 = zpk(b1);
b1 = minreal(b1);

Ab = -A;


%% Voltage transfer function

% H_Amp = ain+(bin*Ab*a1)/(1-Ab*b1);
% H_Amp = zpk(H_Amp);
% H_Amp = minreal(H_Amp);
% 
% figure(1);
% title("H_amp");
% bode(H_Amp);


%% Tia transfer function

H_tia = (ain*(1-b1*Ab)+bin*Ab*a1)*Zin/(1-b1*Ab-bin);
H_tia = zpk(H_tia);
H_tia = minreal(H_tia);

figure(2);
bode(H_tia);
title('TIA Transfer function', 'FontSize', 14);

%% Loop gain
% 
% LoopGain = -Ab * b1;
% LoopGain = minreal(zpk(LoopGain));
% 
% figure;
% bode(LoopGain);
% title('Loop Gain (Ab * b1)');
% grid on;
% 
% [Gm, Pm, Wcg, Wcp] = margin(LoopGain);
% fprintf('Gain Margin: %.2f dB at %.2f Hz\n', 20*log10(Gm), Wcg/(2*pi));
% fprintf('Phase Margin: %.2f° at %.2f Hz\n', Pm, Wcp/(2*pi));



