%% ============================================================
%  SNR vs Optical Input Power (dBm) — Shot vs Thermal vs Combined
%  - Generate optical field as sine with target RMS power (verified)
%  - Magnitude-square detection -> optical power
%  - Photocurrent = Rd * P
%  - Add shot noise + thermal noise (white, PSD-based)
%  - Compute SNR for: shot-only, thermal-only, combined
% ============================================================


% Constants
k = Constant.Boltzmann;
q = Constant.ElementaryCharge;

% Receiver / PD parameters
T  = 20 + 273.15;      % K
R  = 50;               % Ohm (front-end/load)
Rd = 0.7;              % A/W (responsivity)
Be = 100e9;            % Hz  (electrical noise bandwidth)
Id = 0;                % A   (dark current, optional)

% Sampling for time-domain demo (needs to be >> Be)
fs   = 1e12;           % Hz
N    = 2^14;           % samples
t    = (0:N-1).'/fs;

% Choose an electrical tone within bandwidth (arbitrary for demo)
f0 = 10e9;             % Hz

% Thermal current PSD (two-sided) and variance in Be
Si_th    = 4*k*T/R;            % A^2/Hz (two-sided)
sigma2_th = Si_th * Be;        % A^2
sigma_th  = sqrt(sigma2_th);   % A_rms

% Optical input power sweep (in dBm)
P_dBm = linspace(-40, 10, 300);
P_W   = 10.^((P_dBm - 30)/10); % W

% Pre-allocate results
SNR_shot_dB = zeros(size(P_dBm));
SNR_th_dB   = zeros(size(P_dBm));
SNR_tot_dB  = zeros(size(P_dBm));

% --- Main loop over optical input power
for ii = 1:numel(P_W)
    Pavg = P_W(ii);

    % Optical field with RMS power = Pavg:
    % Let x(t) be the optical field amplitude such that |x|^2 has mean Pavg.
    % Use a sinusoid: x(t) = A*sin(2*pi*f0*t), then mean(|x|^2)=A^2/2.
    A = sqrt(2*Pavg);                      % -> mean(|x|^2) = Pavg

    x = A * sin(2*pi*f0*t);                % "optical field" (real for simplicity)

    % Verify RMS/mean power numerically (optional)
    P_meas = mean(abs(x).^2);              % should be ~ Pavg

    a = P_meas - Pavg;
    assert(a<1);

    % Magnitude-square detection -> optical power waveform
    Popt = abs(x).^2;                      % W (instantaneous)

    % Photocurrent waveform (includes DC + 2f0 component for this demo)
    I_sig = Rd * Popt;                     % A

    % Define "signal power" as the mean-squared photocurrent due to signal
    % (for this sine-squared waveform, it's OK for a demo SNR definition)
    P_sig = mean(I_sig.^2);

    % -------- Shot noise (white) --------
    % Two-sided PSD: Si_shot = 2*q*(I_photo + I_dark)
    % For this demo, use average current to set the white-noise level:
    Ibar = mean(I_sig) + Id;
    Si_shot   = 2*q*Ibar;                  % A^2/Hz
    sigma2_sh = Si_shot * Be;              % A^2
    sigma_sh  = sqrt(sigma2_sh);           % A_rms

    n_sh = sigma_sh * randn(N,1);          % time-domain shot noise

    % -------- Thermal noise (white Gaussian) --------
    n_th = sigma_th * randn(N,1);          % time-domain thermal noise

    % Noise powers (mean-square) in time domain
    Pn_sh = mean(n_sh.^2);
    Pn_th = mean(n_th.^2);
    Pn_tot = mean((n_sh + n_th).^2);       % ~ Pn_sh + Pn_th (independent)

    % SNRs
    SNR_shot = P_sig / Pn_sh;
    SNR_th   = P_sig / Pn_th;
    SNR_tot  = P_sig / Pn_tot;

    SNR_shot_dB(ii) = 10*log10(SNR_shot);
    SNR_th_dB(ii)   = 10*log10(SNR_th);
    SNR_tot_dB(ii)  = 10*log10(SNR_tot);

    % Optional sanity check (can comment out)
    % if ii == 1
    %     fprintf('Check Pavg target/meas: %.3g W / %.3g W\n', Pavg, P_meas);
    % end
end

% Plot comparison
figure(1); clf; hold on;
plot(P_dBm, SNR_shot_dB, 'LineWidth', 1.5);
plot(P_dBm, SNR_th_dB,   'LineWidth', 1.5);
plot(P_dBm, SNR_tot_dB,  'LineWidth', 1.5);
grid on;
xlabel('Optical input power [dBm]');
ylabel('SNR [dB]');
title('SNR vs Optical Input Power: Shot vs Thermal vs Combined');
legend('Shot-noise only','Thermal-noise only','Shot + Thermal','Location','best');

% Print example at 0 dBm
[~,idx0] = min(abs(P_dBm - 0));
fprintf('At 0 dBm:\n');
fprintf('  SNR (shot only)    = %.2f dB\n', SNR_shot_dB(idx0));
fprintf('  SNR (thermal only) = %.2f dB\n', SNR_th_dB(idx0));
fprintf('  SNR (combined)     = %.2f dB\n', SNR_tot_dB(idx0));

% Also print NEP based on thermal PSD (constant NEP_th)
NEP_th = sqrt(Si_th)/Rd; % W/sqrt(Hz)
fprintf('Thermal NEP = %.3g W/sqrt(Hz)\n', NEP_th);