What to change in the given code for Plotting SPR curve between R and wavelength?

In this code Drude model is used to calculate the complex permitivity of metal. After running this 100% reflectivity is coming.Kindly tell what other changes is required? Here is code:

clear all
close all
c = 299792458; % Speed of light (m/s)
lambda = linspace(200e-9, 1000e-9, 200);
%twopic = 1.883651567308853e+09; % twopic=2*pi*c where c is speed of light
%omegalight = twopic*(lambda.^(-1)); % angular frequency of light (rad/s)
invsqrt2 = 0.707106781186547;  % 1/sqrt(2)
ehbar = 1.51926751447914e+015; % e/hbar where hbar=h/(2*pi) and e=1.6e-19
omegap = 9.03*ehbar;
gamma=1e10;
n_prism = 1.7786;
theta=37.2;
THETA = theta / 180 * pi;
d_gold = 50e-9;
n_subphase=1.33;
en = [n_prism, n_subphase];
ek = [0, 0];
er=en(1)^2-ek(1)^2;
ei=2*en(1)*ek(1);
e(1)=complex(er,ei);
er=en(2)^2-ek(2)^2;
ei=2*en(2)*ek(2);
e(2)=complex(er,ei);
q1=sqrt(e(1)-en(1)^2*sin(THETA)^2)/e(1);
q2=sqrt(e(2)-en(1)^2*sin(THETA)^2)/e(2);
% Pre-allocate epsilon array
epsilon_values = zeros(1, length(lambda)); % Array to store 200 epsilon values
beta_values = zeros(1, length(lambda)); % Array to store 200 beta values
q_values=zeros(1,length(lambda));
em = zeros(length(lambda), 2, 2);
% Loop over each wavelength
for j = 1:length(lambda)
    lambda1 = lambda(j); % Current wavelength
    omegal = 2 * pi * c / lambda1; % Angular frequency for the current wavelength
    
    % Calculate epsilon using the Drude model
    epsilon_D = 1 - (omegap^2) / (omegal^2 + 1i * gamma * omegal);
    epsilon_values(j) = epsilon_D; % Store the epsilon value
    % Calculate beta for the current epsilon
    beta = (d_gold * 2 * pi / lambda1) * sqrt(epsilon_D - en(1)^2 * sin(THETA)^2);
    beta_values(j) = beta; % Store the beta value
    q=sqrt(epsilon_D-en(1)^2*sin(THETA)^2)/epsilon_D;
    q_values(j)=q;
    % Fill the em matrix components
    em(j, 1, 1) = cos(beta);
    em(j, 1, 2) = -1i * sin(beta) / q;
    em(j, 2, 1) = -1i * sin(beta) * q;
    em(j, 2, 2) = cos(beta);
end
% Initialize the total transfer matrix
emtot = eye(2); % 2x2 Identity matrix
% Loop to compute the total transfer matrix emtot
for j = 2:(length(epsilon_values) - 1)
    emtot1 = squeeze(em(j, :, :)); % Extract em(j, :, :)
    emtot = emtot * emtot1;
end
for j = 1:length(lambda)
    rp = ((emtot(1, 1) + emtot(1, 2) * q2) * q1 - (emtot(2, 1) + emtot(2, 2) * q2)) / ...
        ((emtot(1, 1) + emtot(1, 2) * q2) * q1 + (emtot(2, 1) + emtot(2, 2) * q2));
    ref_values(j) = abs(rp)^2; % Reflectance
end
% Plot reflectance vs wavelength
figure;
plot(lambda * 1e9, ref_values); % Convert lambda to nm for plotting
xlabel('Wavelength (nm)', 'FontSize', 12);
ylabel('Reflectance (R)', 'FontSize', 12);
title('Reflectance vs Wavelength', 'FontSize', 14);
grid on;

Answer 1

[Cris LaPierre]

I believe the issue you have is that the values used to compute rp never change, which means it computes the same ref_value in each loop. To plot the Surface Plasmon Resonance (SPR) curve between reflectivity ( R ) and wavelength using the Drude model, you need to modify your code to calculate the reflectivity as a function of wavelength and then plot it.

Here is an example I generated using your question in the MATLAB AI Chat Playground.

% Parameters for the Drude model
epsilon_inf = 1; % High-frequency permittivity
omega_p = 1e15; % Plasma frequency (rad/s)
gamma = 1e13; % Damping constant (rad/s)

% Wavelength range (in meters)
lambda = linspace(200e-9, 1000e-9, 200);

% Calculate the complex permittivity
omega = 2 * pi * 3e8 ./ lambda; % Angular frequency
epsilon = epsilon_inf - (omega_p^2 ./ (omega.^2 + 1i * gamma * omega));

% Calculate reflectivity R using Fresnel equations
n0 = 1; % Refractive index of air
n1 = sqrt(epsilon); % Refractive index of metal
R = abs((n0 - n1) ./ (n0 + n1)).^2; % Reflectivity

% Plotting the SPR curve
figure;
plot(lambda * 1e9, R, 'LineWidth', 2); % Convert wavelength to nm for plotting
xlabel('Wavelength (nm)');
ylabel('Reflectivity (R)');
title('SPR Curve');
grid on;

●

[Cris LaPierre]

Said another way, you compute a total transfer matrix, then use that to compute reflectance for each wavelength. Since emtot is not wavelength specific, you always get the same reflectance value. I'm not an expert in this topic, but it seems if you instead use em, that has a different transfer matrix for each wavelength, so you would obtain wavelength-specific reflectance values.

Here is your code with that change made.

c = 299792458; % Speed of light (m/s)
lambda = linspace(200e-9, 1000e-9, 200);
ehbar = 1.51926751447914e+015; % e/hbar where hbar=h/(2*pi) and e=1.6e-19
omegap = 9.03*ehbar;
gamma=1e10;
n_prism = 1.7786;
theta=37.2;
THETA = theta / 180 * pi;
d_gold = 50e-9;
n_subphase=1.33;
en = [n_prism, n_subphase];
ek = [0, 0];
er=en(1)^2-ek(1)^2;
ei=2*en(1)*ek(1);
e(1)=complex(er,ei);
er=en(2)^2-ek(2)^2;
ei=2*en(2)*ek(2);
e(2)=complex(er,ei);
q1=sqrt(e(1)-en(1)^2*sin(THETA)^2)/e(1);
q2=sqrt(e(2)-en(1)^2*sin(THETA)^2)/e(2);

% Pre-allocate epsilon array
epsilon_values = zeros(1, length(lambda)); % Array to store 200 epsilon values
beta_values = zeros(1, length(lambda)); % Array to store 200 beta values
q_values=zeros(1,length(lambda));
em = zeros(length(lambda), 2, 2);

omegal = 2 * pi * c ./ lambda; % Angular frequency
% Calculate epsilon using the Drude model
epsilon_values = 1 - (omegap^2) ./ (omegal.^2 + 1i * gamma * omegal);

% Calculate beta for the current epsilon
beta_values = (d_gold * 2 * pi ./ lambda) .* sqrt(epsilon_values - en(1)^2 * sin(THETA)^2);
q_values=sqrt(epsilon_values-en(1)^2*sin(THETA)^2)./epsilon_values;

% Fill the em matrix components
em = [cos(beta_values);
    -1i * sin(beta_values) .* q_values;
    -1i * sin(beta_values) ./ q_values;
    cos(beta_values)];
em = reshape(em,2,2,[]);

% Compute reflectance for each wavelength
ref_values = ((em(1, 1,:) + em(1, 2,:) * q2) * q1 - (em(2, 1,:) + em(2, 2,:) * q2)) ./ ...
        ((em(1, 1,:) + em(1, 2,:) * q2) * q1 + (em(2, 1,:) + em(2, 2,:) * q2));
ref_values = abs(squeeze(ref_values)).^2;

% Plot reflectance vs wavelength
figure;
plot(lambda * 1e9, ref_values); % Convert lambda to nm for plotting
xlabel('Wavelength (nm)', 'FontSize', 12);
ylabel('Reflectance (R)', 'FontSize', 12);
title('Reflectance vs Wavelength', 'FontSize', 14);
grid on;