Nsamples = 1E2; % Plotting parameter

N = 0.3E26; % Density [m^{-3}] at room temp and pressure
e = 1.6E-19; % Fundamental charge [C]
e0 = 8.9E-12; % Vacuum permittivity [F/m]
me = 9.1E-31; % Electron mass [kg]

wp = sqrt(N*e^2/(e0*me));
gamma = 2*pi*10E9; % damping coefficient
w0 = 2*pi*7.1E15; % undamped frequency
w = 2*pi*linspace(5E15,9E15,Nsamples);

% electric susceptibility from classical dipole model
xe = wp^2./((w0^2-w.^2).^2+2*gamma*w*1j);
n = sqrt(1+xe);
eta = real(n);
kap = imag(n);

plot(1/(2*pi)*w*1E-12,eta); hold on; grid on;
plot(1/(2*pi)*w*1E-12,kap,'r--'); 
xlabel('Drive frequency, f=\omega/2\pi [THz]');
ylabel('\eta = Re(n) [blue]; and \kappa = Im(n) [red]');
title('Problem 5c: Index of refraction n(\omega) = \eta(\omega) + j*\kappa(\omega)');
xbounds = [min(w) max(w)];
xlim(1/(2*pi)*xbounds*1E-12);