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); % 2×2 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;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); % 2×2 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; 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

Cart

Cart