% Minimal MZM transfer-function demo (sinusoidal drive) — aligned with your notation
%
% Implements exactly:
%   E_out(t) = E0 * exp(j*w0*t) * exp(-j*w0*L*n_eff/c0) * 1/2 * [ exp(-j*phi1(t)) + rho*exp(-j*phi2(t)) ]
% with  phi_{1,2}(t) = pi * v_{1,2}(t)/Vpi
%
% Push-pull:
%   v1(t) = +v_drive(t)/2 ,  v2(t) = -v_drive(t)/2   =>  phi1 = +pi/2 * v_drive/Vpi, phi2 = -pi/2 * v_drive/Vpi
%
% And the ideal TF (rho=1):
%   E_out/E_in = exp(-j*w0*L*n_eff/c0) * cos( (pi/2) * v_drive/Vpi )
%
% Note: E_in(t) = E0 * exp(j*w0*t) in this script.

% Parameters
c0    = physconst('lightspeed');      % [m/s]
lambda0 = 1310e-9;      % [m]
omega0  = 2*pi*c0/lambda0;

L     = 5e-3;           % [m] effective phase section length (set as needed)
n_eff = 2.2;            % [-] effective index (set as needed)

E0    = 1;              % field amplitude (arbitrary)
Vpi   = 3.2;            % [V] half-wave voltage (your V_pi)

% Drive
f0    = 1e9;            % [Hz]
fs    = 100e9;          % [Hz]
Nper  = 1;              % number of periods
Vpp   = 0.5*Vpi;        % [V] peak-to-peak of v_drive(t)

biasV = 2;            % [V] differential bias added to v_drive

% Analytic
v_ = linspace(-1,2, 2001);
% Field transfer function (amplitude)
Field_mzm_analytic = cos((pi/2)*v_);
% Power transfer function (intensity)
P_mzm_analytic = Field_mzm_analytic.^2;

% Imbalance factor in YOUR notation:
rho   = 1;            % rho=1 -> ideal balanced MZM (collapses to ideal TF)

% Time axis + differential drive voltage v_drive(t)
T  = Nper/f0;
t  = (0:1/fs:T-1/fs).';

v_drive = biasV + (Vpp/2)*sin(2*pi*f0*t);  % v_drive(t) (peak = Vpp/2)

% Push-pull branch voltages (consistent with v_drive = v1 - v2)
v1 = +0.5*v_drive;   % arm 1
v2 = -0.5*v_drive;   % arm 2

% Phases phi1, phi2
phi1 = pi * v1 / Vpi;
phi2 = pi * v2 / Vpi;

% Fields: E_in and E_out (exactly your Eq. (mzm_e_field))
E_in  = E0 .* exp(1i*omega0*t);

common_phase = exp(-1i * (omega0*L*n_eff/c0));  % exp(-j*omega0*L*n_eff/c0)

E_out = E0 .* exp(1i*omega0*t) .* common_phase .* 0.5 .* ...
        ( exp(-1i*phi1) + rho .* exp(-1i*phi2) );

% Transfer function (numerical): E_out/E_in
H_num = E_out ./ E_in;

% Power (normalized)
Pnorm_num   = abs(H_num).^2;     % since |E_out/E_in|^2



% Ideal TF (analytic) for comparison (rho=1, push-pull)
H_ideal = common_phase .* cos( (pi/2) * (v_drive./Vpi) );

Pnorm_ideal = abs(H_ideal).^2;
Pnorm_math = cos( (pi/2) * (v_drive./Vpi) ).^2;





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

% Normalized voltage axis (multiples of Vpi)
v_norm = v_drive./Vpi;

colfield = [0,0,0]; %is black
colpow = linspecer(2);
colpow = colpow(1,:);
colvdrive = linspecer(2);
colvdrive = colvdrive(2,:);

%% SIGNAL IN
figure(1); clf
plot(v_norm,t*1e9, 'LineWidth', 1.0,'Color',colvdrive); grid on;
ylabel('t [ns]'); xlabel('v_{drive}(t)/V_\pi'); 
title('Drive voltage (normalized)');
xlim([min(v_) max(v_)]);

%% IN/OUT (static transfer) — normalized x-axis + analytic curve
figure(2); clf
plot(v_, Field_mzm_analytic, 'LineWidth', 1.2,'LineStyle','--','Color',colfield); hold on;% analytic power TF
plot(v_, P_mzm_analytic, 'LineWidth', 1.2, 'Color',colpow); hold on;% analytic power TF
% show input time signal
plot(v_norm,-1+t*1e9, 'LineWidth', 1.0,'Color',colvdrive); grid on;
% show output time signal
plot(2+t*1e9, Pnorm_num,   'LineWidth', 1.0,'DisplayName','Intensity', 'Color',colvdrive); hold on;
plot(2+t*1e9, real(H_ideal), '--', 'LineWidth', 1.0,'DisplayName','Field','Color',colfield); hold on;
scatter(v_norm, Pnorm_num, 12, '.', 'LineWidth', 1,'MarkerEdgeColor',colvdrive); 
scatter(biasV./Vpi,(cos((pi/2)*biasV./Vpi)^2),10,'Marker','o');
line([min(v_drive), min(v_drive)]./Vpi,[(cos((pi/2)*min(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--');
line([max(v_drive) max(v_drive)]./Vpi,[(cos((pi/2)*max(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--');
xline([min(v_norm) max(v_norm)])

grid on;
xlabel('v_{drive}(t)/V_\pi'); ylabel('|E_{out}/E_{in}|^2');
% legend
xlim([min(v_) max(v_)+1]);
ylim([-1 1]);

% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mzm.tex');


%%
% % FIELD TF (only field here; do not mix power into this figure)
figure(3); clf
% plot(t*1e9, real(H_num), 'LineWidth', 1.0); hold on;
% plot(t*1e9, real(H_ideal), '--', 'LineWidth', 1.0,'DisplayName','Field','Color',colfield); hold on;
plot(t*1e9, Pnorm_num,   'LineWidth', 1.0,'DisplayName','Intensity', 'Color',colpow); hold on;
grid on;
xlabel('t [ns]'); ylabel('Re\{E_{out}/E_{in}\}');
legend
mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mzm_out.tex');