N = 0.3E26;     % Hydrogen density [m^{-3}]

e0 = 1.6E-19;   % Electron charge [C]

eps0 = 8.9E-12; % Vacuum permittivity [F/m]

me = 9.1E-31;   % Electron mass [kg]



wp = sqrt(N*e0^2/(eps0*me)); % Plasma frequency

gamma = 2*pi*(100E9); % Damping freq [s^-1]

w0 = 2*pi*(7.1E15); 



L = 3;

Npoints = 1e3;

w = linspace(w0-L*gamma,w0+L*gamma,Npoints);

f = 1/(2*pi)*w;

fBounds = [min(f) max(f)];



eta = 1 + 1/2*wp^2*(w0^2-w.^2)./((w0^2-w.^2).^2+4*gamma^2*w.^2);

kap = -wp^2*gamma*w./((w0^2-w.^2).^2+4*gamma^2*w.^2);

plot(f,eta); hold on; grid on;

plot(f,kap,'r--');

xlim(fBounds);

xlabel('Drive frequency f=1/(2\pi)\omega [THz]');

ylabel('Real and imaginary parts of the index of refraction: n = \eta + j\kappa');

title('6.621 PS#1, Problem 5(f)');

legend('Real part, \eta','Imag part, \kappa');