clear;
clc;

% Sweep the degenerate pump frequency around its nominal wavelength.
pump_detuning_hz = (-800:0.01:800) .* 1e9;

% Degenerate FWM setup: two pump photons at f_p and one signal at f_s
% generate an idler at f_i = 2*f_p - f_s.
pump_wavelength_nm = 1310;
signal_wavelength_nm = 1308;
zdw_wavelength_nm = 1310;

f_pump_nominal = wavelength2frequency(pump_wavelength_nm, 'nm');
f_signal_scalar = wavelength2frequency(signal_wavelength_nm, 'nm');
f_zdw_scalar = wavelength2frequency(zdw_wavelength_nm, 'nm');

f_pump = f_pump_nominal + pump_detuning_hz;
f_signal = f_signal_scalar .* ones(size(f_pump));
f_zdw = f_zdw_scalar .* ones(size(f_pump));
f_idler = 2 .* f_pump - f_signal;

% Fiber parameters
dispersion_slope_ps_nm2_km = 0.07;
attenuation_db_per_km = 0.21;
fiber_length_km = 10;

% Launch powers and nonlinear coefficient
pump_power_dbm = 10;
signal_power_dbm = 10;
pump_power_w = dbm2watt(pump_power_dbm);
signal_power_w = dbm2watt(signal_power_dbm);
gamma_w_inv_km_inv = 1.3;
degeneracy_factor = 3;

[P_fwm, eta, delta_beta, L_eff_m] = calcFwmPower( ...
    f_pump, f_pump, f_signal, f_zdw, ...
    dispersion_slope_ps_nm2_km, attenuation_db_per_km, fiber_length_km, ...
    pump_power_w, pump_power_w, signal_power_w, ...
    gamma_w_inv_km_inv, degeneracy_factor);

f_pump_thz = f_pump .* 1e-12;
f_zdw_thz = f_zdw_scalar .* 1e-12;
f_signal_thz = f_signal_scalar .* 1e-12;
idler_power_dbm = 10 .* log10(max(P_fwm, realmin) ./ 1e-3);

figure;
tiledlayout(2,1);

ax1 = nexttile;
plot(ax1, f_pump_thz, eta, 'LineWidth', 2);
hold(ax1, 'on');
xline(ax1, f_zdw_thz, '--r', 'ZDW', 'LineWidth', 1.2, ...
    'LabelOrientation', 'horizontal', 'LabelVerticalAlignment', 'bottom');
xline(ax1, f_signal_thz, '--k', 'Signal', 'LineWidth', 1.2, ...
    'LabelOrientation', 'horizontal', 'LabelVerticalAlignment', 'middle');
ylabel(ax1, 'FWM efficiency');
grid(ax1, 'on');
title(ax1, 'FWM Efficiency and Idler Power versus Pump Frequency');

ax2 = nexttile;
plot(ax2, f_pump_thz, idler_power_dbm, 'LineWidth', 2);
hold(ax2, 'on');
xline(ax2, f_zdw_thz, '--r', 'ZDW', 'LineWidth', 1.2, ...
    'LabelOrientation', 'horizontal', 'LabelVerticalAlignment', 'bottom');
xline(ax2, f_signal_thz, '--k', 'Signal', 'LineWidth', 1.2, ...
    'LabelOrientation', 'horizontal', 'LabelVerticalAlignment', 'middle');
xlabel(ax2, 'Pump frequency (THz)');
ylabel(ax2, 'FWM idler power (dBm)');
grid(ax2, 'on');

fprintf('Pump wavelength        : %.3f nm -> %.6f THz\n', ...
    pump_wavelength_nm, f_pump_nominal .* 1e-12);
fprintf('Signal wavelength      : %.3f nm -> %.6f THz\n', ...
    signal_wavelength_nm, f_signal_scalar .* 1e-12);
fprintf('ZDW wavelength         : %.3f nm -> %.6f THz\n', ...
    zdw_wavelength_nm, f_zdw_scalar .* 1e-12);
fprintf('Pump launch power      : %.2f dBm -> %.4g W\n', ...
    pump_power_dbm, pump_power_w);
fprintf('Signal launch power    : %.2f dBm -> %.4g W\n', ...
    signal_power_dbm, signal_power_w);
fprintf('Peak FWM efficiency    : %.4g\n', max(eta));
fprintf('Peak FWM idler power   : %.4g W (%.2f dBm)\n', ...
    max(P_fwm), 10 .* log10(max(P_fwm) ./ 1e-3));
fprintf('Effective fiber length : %.4f km\n', L_eff_m ./ 1e3);
fprintf('Idler wavelength range : %.3f nm to %.3f nm\n', ...
    min(frequency2wavelength(f_idler, 'nm')), max(frequency2wavelength(f_idler, 'nm')));
fprintf('Max |delta beta|       : %.4g 1/m\n', max(abs(delta_beta)));

function wavelength = frequency2wavelength(frequency, outputUnit)
    c = physconst('LightSpeed');
    wavelength = c ./ frequency;

    switch lower(outputUnit)
        case 'm'
        case 'nm'
            wavelength = wavelength .* 1e9;
        otherwise
            error('Unsupported output unit "%s". Use "m" or "nm".', outputUnit);
    end
end

function power_w = dbm2watt(power_dbm)
    power_w = 1e-3 .* 10.^(power_dbm ./ 10);
end