%------------------------------------------------------------
% Tony Hyun Kim
% Problem 5e, Hydrogen absorption
%------------------------------------------------------------
clear all; close all;
Nsamples = 1E4; % 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)); % Plasma frequency
gamma = 2*pi*10E9; % damping coefficient
w0 = 2*pi*7.1E15; % undamped frequency
c0 = 3E8; % vacuum speed of light
w = 2*pi*linspace(7.09995E15,7.10005E15,Nsamples);

a = (gamma*wp^2/c0)*w.^2./((w0^2-w.^2).^2+4*gamma^2*w.^2);
plot(1/(2*pi)*w*1E-12,a);
xlabel('Drive frequency, f=\omega/2\pi [THz]');
ylabel('Absorption \alpha [m^{-1}]');
title('Problem 5b: \gamma = 2\pi(10GHz)');
xbounds = [min(w) max(w)];
xlim(1/(2*pi)*xbounds*1E-12);

% Analytic results
a_max = 1/4*wp^2/(c0*gamma);
w_L = sqrt((w0^2+2*gamma^2)-2*gamma*sqrt(gamma^2+w0^2));
w_R = sqrt((w0^2+2*gamma^2)+2*gamma*sqrt(gamma^2+w0^2));
FWHM = w_R-w_L;
fprintf(1,'Analytic a_max: %f\n',a_max);
fprintf(1,'Analytic w_L: %f\n',w_L);
fprintf(1,'Analytic w_R: %f\n',w_R);
fprintf(1,'Analytic FWHM (Hz): %f\n',1/(2*pi)*FWHM);
hold on;
plot(1/(2*pi)*xbounds*1E-12,1/2*a_max*[1 1],'r-'); grid on;
