I need help on how to plot Surface Plasmon Polaritons (SPPs) for Multilayer Systems

12 次查看(过去 30 天)
I have derived the surface plasmon polariton (SPP) relation for the system under consideration. The system consists of three layers: semi-infinite layers of InSb as the upper layer and SiC (lower layer), separated by a vacuum gap of 10 nm. I have written a code to calculate the parallel wavevector (G) and plot it against frequency (ω). However, the results seem incorrect, as I only obtain numerical values for some frequency points. Could you help me identify and correct the error? Here is my code:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
clc
% Physical parameters
hbar = 1.054571596e-34;
h = 2 * pi * hbar;
global c0;
c0 = 2.99792458e+8;
kb = 1.3806503e-23;
u0 = 4 * pi * 10^-7;
ep0 = 8.854 * 10^-12;
qe = 1.602176462e-19;
T_h0 = [350];
T_l = [300];
T_SiC = T_h0;
Tl_SiC = T_l;
T_InSb = T_h0;
Tl_InSb = T_l;
d = 10e-9;
N = 200;
% Initialization of different parameters
w1 = 0.01e14;
w2 = 5e14;
w = linspace(w1, w2, N);
lambda = (2 * pi * c0) ./ w;
w_ev = w / 1.5193E15;
ep1 = 1+0*w; %gap
ep2 = epsilon_SiC1w_T(w, Tl_SiC);
ep3 = epsilon_InSbw_T(w, T_InSb);
% Initialize arrays to store solutions
G_solutions = zeros(N, 1);
abs_G_solutions = zeros(N, 1);
for i = 1:N
syms G
k0 = w(i) / c0; % vacuum
k1=k0; %gap
k2 = sqrt(G^2 - ep2(i) * k0.^2); % SiC material
k3 = sqrt(G^2 - ep3(i) * k0.^2); % InSb material
% Dispersion equation for three-layer system: InSb/gap/SiC
equation = exp(-2 * k1 * d) - (((k2 * ep1 + k1 * ep2(i)) .* (k3 * ep1 + k1 * ep3(i))) ./ ...
((k2 * ep1 - k1 * ep2(i)) .* (k3 * ep1 - k1 * ep3(i)))) == 0;
% Solve the equation
solutions = solve(equation, G);
numerical_solutions = vpa(solutions, 10);
% Check if numerical_solutions is empty
if ~isempty(numerical_solutions)
G_solutions(i) = double(numerical_solutions(1)); % Assuming first solution is desired
abs_G_solutions(i) = abs(G_solutions(i));
else
G_solutions(i) = NaN; % Assign NaN if no solution found
abs_G_solutions(i) = NaN;
end
end
% Display the numerical solutions
disp('Numerical solutions for G at different frequencies:');
disp(G_solutions);
disp('Absolute values of G at different frequencies:');
disp(abs_G_solutions);
% Plot the scatter plot
figure(1)
scatter(w, abs_G_solutions, 'filled')
xlabel('Frequency (w)')
ylabel('Absolute value of G (|G|)')
title('Scatter plot of |G| vs Frequency (w)')
grid on
figure(2)
set(gcf, 'Units', 'pixels', 'Position', [100, 100, 500, 300]);
scatter(w,abs_G_solutions,'blue','filled','LineWidth',3);
xlabel('Energy (eV)','FontSize',14);
ylabel('Absolute value of G (|G|)')
set(gca,'linewidth',3);
set(gca,'FontSize',14);
set(0, 'DefaultAxesFontName', 'Times New Roman');

回答(1 个)

Ronit
Ronit 2024-6-21
Hello Ambali,
To address the issues, you're encountering with plotting Surface Plasmon Polaritons (SPPs) for a multilayer system consisting of InSb and SiC layers separated by a vacuum gap, it seems that your MATLAB code has several key areas for improvement.
% Updated syms declaration to specify G as positive for physical solutions
syms G positive
% Updated vpasolve call with a specified range for G
solutions = vpasolve(equation, G, [0, inf]); % Limiting G to positive values
% Store the first solution if it exists
if ~isempty(solutions)
G_solutions(i) = double(solutions(1));
abs_G_solutions(i) = abs(G_solutions(i));
else
G_solutions(i) = NaN; % Assign NaN if no solution is found
abs_G_solutions(i) = NaN;
end
Explanation for the above changes:
  • Positive ‘G’ Declaration: By declaring ‘G’ as positive in the symbolic variable declaration, the solver is directed to look for positive solutions, aligning with the physical expectation that the wavevector component should be real and positive for propagating modes.
  • Specifying Range in ‘vpasolve’: The updated ‘vpasolve’ call includes a range [0, inf] for ‘G’. This helps the solver by providing a physical range of interest, potentially improving the accuracy and efficiency of the solution process.
  • Handling No Solution Cases: The conditional check for empty solutions ‘(if ~isempty(solutions))’ ensures that the code can gracefully handle scenarios where no solution is found by the solver, assigning ‘NaN’ to the ‘G_solutions’ and ‘abs_G_solutions’ arrays for those frequencies. This prevents the code from crashing or throwing errors when solutions are not found.
Hope this helps!

类别

Help CenterFile Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

产品


版本

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by