%MZM demo -> sinus als eingang in MZM intensity TF: 2nd and 3rd roder
%nonlinearities in PSD visible

clear; close all; clc;

set(groot,'defaultLegendInterpreter','tex');
set(groot,'defaultAxesTickLabelInterpreter','tex');
set(groot,'defaultTextInterpreter','tex');

%% Fixed parameters
Vpi  = 5.2;     % [V]
f0   = 10e9;     % [Hz]
fs   = 400e9;   % [Hz]
Nper = 500;     % periods for PSD quality

t = (0:1/fs:(Nper/f0 - 1/fs)).';
w = 2*pi*f0;

% PSD settings
nfft = 2^(nextpow2(min(length(t), 2^18))-1);
win  = hann(2^12);
ovl  = round(0.5*numel(win));

N_bessel = 10;

%% UI defaults (normalized)
vb0  = 1.0;   % Vbias/Vpi
vpp0 = 0.5;   % Vpp/Vpi

%% Figure + layout
fig = figure('Color','w','Name','MZM Nonlinearity: Bias & Drive','NumberTitle','off');
tl  = tiledlayout(fig,1,2,'TileSpacing','compact','Padding','compact');

axTF  = nexttile(tl,1); hold(axTF,'on'); grid(axTF,'on');
axPSD = nexttile(tl,2); hold(axPSD,'on'); grid(axPSD,'on');

% Scatter placeholders
hEx   = scatter(axTF, nan, nan, 6, '.', 'DisplayName','Exact');
hTa   = scatter(axTF, nan, nan, 6, '.', 'DisplayName','Taylor (3rd order)');
hJa   = scatter(axTF, nan, nan, 6, '.', 'DisplayName',sprintf('Jacobi--Anger (N=%d)',N_bessel));
hBias = plot(axTF, nan, nan, 'ko', 'MarkerFaceColor','k', 'DisplayName','Bias');

xlabel(axTF,'v/V_\pi'); ylabel(axTF,'P_{out}/P_0');   % <-- TeX (no $...$)
title(axTF,'Transfer characteristic (scatter)');
ylim(axTF,[-0.1 1.1]);
xlim(axTF,[0 2]);
legend(axTF,'Location','best');

% PSD placeholders
hPex = plot(axPSD, nan, nan, 'LineWidth',2.0, 'DisplayName','Exact');
hPta = plot(axPSD, nan, nan, '-', 'LineWidth',1.5, 'DisplayName','Taylor (3rd order)');
hPja = plot(axPSD, nan, nan, '--',  'LineWidth',0.1, 'DisplayName',sprintf('Jacobi--Anger (N=%d)',N_bessel));

xlabel(axPSD,'Frequency [GHz]'); ylabel(axPSD,'PSD [dB/Hz]');
title(axPSD,'Output spectrum (PSD)');
xlim(axPSD,[0 10*f0/1e9]);
legend(axPSD,'Location','best');

%% Sliders + labels
sH = 0.05;  mL = 0.08;  wS = 0.38;  y1 = 0.04;  dy = 0.06;

uicontrol(fig,'Style','text','Units','normalized', ...
    'Position',[mL, y1+dy, wS, 0.03], ...
    'String','v_{bias}/V_{\pi}','HorizontalAlignment','left');

sBias = uicontrol(fig,'Style','slider','Units','normalized', ...
    'Position',[mL, y1+dy-0.02, wS, sH], ...
    'Min',0,'Max',2,'Value',vb0);

tBiasVal = uicontrol(fig,'Style','text','Units','normalized', ...
    'Position',[mL+wS+0.01, y1+dy, 0.08, 0.03], ...
    'String',sprintf('%.3f',vb0),'HorizontalAlignment','left');

uicontrol(fig,'Style','text','Units','normalized', ...
    'Position',[mL, y1, wS, 0.03], ...
    'String','v_{pp}/V_{\pi}','HorizontalAlignment','left');

sVpp = uicontrol(fig,'Style','slider','Units','normalized', ...
    'Position',[mL, y1-0.02, wS, sH], ...
    'Min',0,'Max',2,'Value',vpp0);

tVppVal = uicontrol(fig,'Style','text','Units','normalized', ...
    'Position',[mL+wS+0.01, y1, 0.08, 0.03], ...
    'String',sprintf('%.3f',vpp0),'HorizontalAlignment','left');

%% Store handles in fig.UserData (so callback can always access them)
S = struct();
S.Vpi = Vpi; S.f0 = f0; S.fs = fs; S.w = w; S.t = t;
S.win = win; S.ovl = ovl; S.nfft = nfft;
S.N_bessel = N_bessel;

S.axTF = axTF; S.axPSD = axPSD;
S.hEx = hEx; S.hTa = hTa; S.hJa = hJa; S.hBias = hBias;
S.hPex = hPex; S.hPta = hPta; S.hPja = hPja;

S.sBias = sBias; S.sVpp = sVpp;
S.tBiasVal = tBiasVal; S.tVppVal = tVppVal;

fig.UserData = S;

%% Continuous update while dragging
addlistener(sBias,'Value','PostSet',@(~,~)updatePlots(fig));
addlistener(sVpp ,'Value','PostSet',@(~,~)updatePlots(fig));

% Initial draw
updatePlots(fig);

%% ===== Callback (separate function at end of script) =====
function updatePlots(fig)
    S = fig.UserData;

    % Read slider values (normalized)
    vb_n  = S.sBias.Value;     % Vbias/Vpi
    vpp_n = S.sVpp.Value;      % Vpp/Vpi

    % Update value labels
    S.tBiasVal.String = sprintf('%.3f', vb_n);
    S.tVppVal.String  = sprintf('%.3f', vpp_n);

    % Convert to volts / amplitude
    Vpi   = S.Vpi;
    Vbias = vb_n * Vpi;
    Vpp   = vpp_n * Vpi;
    Vm    = Vpp/2;

    t = S.t; w = S.w;

    % Drive
    v = Vbias + Vm*cos(w*t);
    x = v./Vpi;

    % Exact intensity
    P_exact = cos((pi/2)*x).^2;

    % Taylor 3rd order around Vbias
    k  = (pi/2)/Vpi;
    vb = Vbias;
    g0 = cos(k*vb)^2;
    g1 = -k*sin(2*k*vb);
    g2 = -2*k^2*cos(2*k*vb);
    g3 =  4*k^3*sin(2*k*vb);
    dv = v - vb;
    P_taylor = g0 + g1*dv + 0.5*g2*dv.^2 + (1/6)*g3*dv.^3;

    % Jacobi–Anger / Bessel series (truncated)
    a = pi*(Vbias/Vpi);
    b = pi*(Vm/Vpi);
    N = S.N_bessel;

    P_ja = 0.5*ones(size(t));
    P_ja = P_ja + 0.5*cos(a)*besselj(0,b);

    for m = 0:floor((N-1)/2)
        n = 2*m + 1;
        P_ja = P_ja - (0.5*2)*sin(a)*besselj(n,b).*cos(n*w*t);
    end
    for m = 1:floor(N/2)
        n = 2*m;
        P_ja = P_ja - (0.5*2)*cos(a)*besselj(n,b).*cos(n*w*t);
    end

    % Update TF scatter
    S.hEx.XData = x;  S.hEx.YData = P_exact;
    S.hTa.XData = x;  S.hTa.YData = P_taylor;
    S.hJa.XData = x;  S.hJa.YData = P_ja;

    xb = Vbias/Vpi;
    pb = cos((pi/2)*xb)^2;
    S.hBias.XData = xb; S.hBias.YData = pb;

    xpad = 0.05*(max(x)-min(x) + eps);
    % xlim(S.axTF,[min(x)-xpad, max(x)+xpad]);
    ylim(S.axTF,[-0.1 1.1]);

    % PSDs
    fs = S.fs;
    [Se,f] = pwelch(P_exact-mean(P_exact), S.win, S.ovl, S.nfft, fs, 'onesided');
    [St,~] = pwelch(P_taylor-mean(P_taylor), S.win, S.ovl, S.nfft, fs, 'onesided');
    [Sj,~] = pwelch(P_ja-mean(P_ja),         S.win, S.ovl, S.nfft, fs, 'onesided');

    S.hPex.XData = f/1e9; S.hPex.YData = 10*log10(Se + realmin);
    S.hPta.XData = f/1e9; S.hPta.YData = 10*log10(St + realmin);
    S.hPja.XData = f/1e9; S.hPja.YData = 10*log10(Sj + realmin);

    xlim(S.axPSD,[0 10*(S.f0)/1e9]);
    ylim(S.axPSD,[-180 -80]);

    drawnow limitrate;
end