- 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.

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

8 views (last 30 days)

Show older comments

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');

##### 0 Comments

### Answers (1)

Ronit
on 21 Jun 2024

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:

Hope this helps!

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!