energy=0:0.001:2;                        %in eV
energy0=1;                              %resonant energy in eV
delta_e=0.01;                            %resonance width in eV
c=3e14;                                 %c in micron/s
lambda=1.23984./energy;                 %lambda in microns
lambda0=1.23984./energy0;               %resonant wavelength in microns
delta_lambda=lambda0/energy0*delta_e;   %delta_lambda in microns
omega=2*pi*c./lambda;                   %(angular) frequency in Hz
omega0=2*pi*c/lambda0;                  %(angular) resonant frequency in Hz
gama=omega0/lambda0*delta_lambda;       %gamma in Hz

% 1 cavity Reflection & Transmission 
Tc=zeros(2,2);

for n=1:length(omega)
    
Tc(1,1)=1-i*gama/(omega(n)-omega0);Tc(1,2)=-i*gama/(omega(n)-omega0);
Tc(2,1)=i*gama/(omega(n)-omega0);Tc(2,2)=1+i*gama/(omega(n)-omega0);

b1(n)=-Tc(2,1)/Tc(2,2);                 % Reflected field: 1 cavity
b2(n)=1/Tc(2,2);                        % Transmitted field: 1 cavity

end

% Plot Refectivity 1 cavity. As a function of energy (eV)

figure(1);

subplot(2,3,1);plot(energy,abs(b1).^2); xlabel('Energy (eV)');ylabel('R');
subplot(2,3,4);plot(energy,abs(b2).^2); xlabel('Energy (eV)');ylabel('T');

% Plot Refectivity 1 cavity. As a function of lambda (microns)

subplot(2,3,2);plot(lambda,abs(b1).^2); xlabel('Wavelength (\mu m)');ylabel('R');
subplot(2,3,5);plot(lambda,abs(b2).^2); xlabel('Wavelength (\mu m)');ylabel('T');

% Plot Refectivity 1 cavity. As a function of omega (Hz)

subplot(2,3,3);plot(omega,abs(b1).^2); xlabel('Frequency (Hz)');ylabel('R');
subplot(2,3,6);plot(omega,abs(b2).^2); xlabel('Frequency (Hz)');ylabel('T');