c = 2.998e8; % [m/s]
d = 150e-6;  % [m]
R1 = 0.04;

N = 1000;
lambda = linspace(798,803,N)*1e-9; % [m]

% Derived quantities
t0 = 2*d/c;
w = 2*pi*c./lambda;
sR1 = sqrt(R1);

% Group delay (use analytic formula)
gd_num = t0*(1-R1^2)-2*sR1*(1-R1)*t0*cos(w*t0);
gd_den = (-2*sR1+(1+R1)*cos(w*t0)).^2+((1-R1)*sin(w*t0)).^2;
gd = gd_num./gd_den;
clear gd_num gd_den;
subplot(121);
plot(lambda*1e9,gd*1e12);
xlim([798 803]);
ylabel('Group delay [ps]');
xlabel('Wavelength [nm]');
grid on;
title('Problem 3d: d=150um, R1=0.04');

% Group delay dispersion (numerical)
subplot(122);
dw = abs(w(2)-w(1));
dT = diff(gd);
gdd = dT/dw;
plot(lambda(1:(N-1))*1e9,gdd*1e24);
xlabel('Wavelength [nm]');
ylabel('Group delay dispersion [ps^2]');
xlim([798 803]);

grid on;