How do I fix "Index in position 2 exceeds array bounds (must not exceed 1)."

clc;
clear;
close all;

%% Data
% Define the parameters
qe = 1;
n = 5e16;
p = 0.525;
E_0 = 10e3;
v_s = 3.5e3;
epsilon_d = 50;
ue = 1;
m0 = 9.1094e-31;
m = 0.44 * m0;
epsilon_0 = 1;
h = 4.1356e-15;
hBar = h / (2 * pi);
b = 0.142e-9;
a = 3 * b / (2 * hBar);
tt = 2.6;
delta = 5.7 / 2;
G = sqrt(delta^2 + 5 * tt^2);

% Define the range for ω_p and B^(-1)
omega_p = linspace(0, 10, 150);
B_inv = linspace(0.04, 1, 150);

k_q = omega_p / v_s;
v = (-2 * tt^2 * a * sin(p * a) / G);
tau = ue * p / (qe * v);
tau_1 = tau / 2;
omega_c = qe * B_inv' * v / p;
gamma = 1 + 1i * omega_c * tau_1;
sigma_g = qe^2 * n * v * tau / p;
Six = p^2 / (2 * pi * hBar^2);
Rix = 2 * p / m / v;
b_g = 2 * pi * hBar^2 * epsilon_0 * v / (qe^2 * p);
beta_g = (k_q' .* v) ./ omega_c;

% Initialize the result matrix
j_x = zeros(numel(omega_p), numel(B_inv));

% Calculate sigma_xx
sigma_xx = zeros(numel(omega_p), numel(B_inv));
for r = -3:3
    J_r = besselj(r, beta_g);
    sigma_xx = sigma_xx + (2 * sigma_g' * J_r) ./ (1 - 1i * (omega_p' - r * omega_c) * tau);
end

% Calculate R_x
R_x = zeros(numel(omega_p), numel(B_inv));
for r = -3:3
    J_r = besselj(r, beta_g);
    R_x = R_x + (omega_c' ./ k_q) .* (J_r .^ 2) ./ (1 - 1i * (omega_p' - r * omega_c) * tau);
end

% Calculate g(omega_q, k_q)
g_omega_q_k_q = 1 + 1i / (epsilon_0 * (epsilon_d + 1)) * (sigma_xx ./ (v_s - R_x));

% Calculate the equation for each combination of omega_p and B^-1
for i = 1:numel(omega_p)
    for j = 1:numel(B_inv)
        sum_term = 0;
        for r = -1:1
            J_r = besselj(r, beta_g(i));
            sum_term = sum_term + J_r / (1 - 1i * (omega_p(i) - r * omega_c(j)) * tau);
        end

        j_x(i, j) = (1 / (2 * qe * n * v(i))) * ((tau * omega_c(j) * p^2) / (8 * pi * (hBar^2))) * (abs((sigma_g(i, j) * E_0) / g_omega_q_k_q(i, j)))^2 ...
            * ((1 / beta_g(i)) / (1 + (omega_c(j)^2 * tau_1^2)))^2 * (sum_term) ...
            * (-1i * gamma(i, j)^2 * (sum_term + 1) * besselj(sum_term + 1, beta_g(i)) + 1i * conj(gamma(i, j))^2 * (sum_term - 1) * besselj(sum_term - 1, beta_g(i)));
    end
end

Index in position 2 exceeds array bounds. Index must not exceed 1.

% Plot the results
[Omega_P, B_inv] = meshgrid(omega_p, B_inv);
surf(Omega_P, B_inv, j_x);
xlabel('\omega_p');
ylabel('B^{-1}');
zlabel('j_x');
title('Plot of j_x against \omega_p and B^{-1}');

title('Plot of j_x against \omega_p and B^{-1}');

Answer 1

[Torsten]

Before entering the loop after which the error message pops up, type

size(gamma)
size(sigma_g) 
size(g_omega_q_k_q)

and check whether the second dimension of these arrays is really >= numel(B_inv) as required because of the following loop.

[Samuel Suakye]

thak you but same problem persist on line 73

[Torsten]

Did you fix the error with the arrays mentionned ? The second dimension of all these arrays is 1 in your code from above. And since I don't know how you modified the code, I cannot give you advice concerning the aforementionned new (or persistent) error.

[Samuel Suakye]

I used the repmat function to replicate the values of gamma and sigma_g to match the size of g_omega_q_k_q

[Samuel Suakye]

% Adjust dimensions of gamma and sigma_g
gamma = repmat(1 + 1i * omega_c * tau_1, numel(omega_p), 1);
sigma_g = repmat(qe^2 * n * v * tau / p, 1, numel(B_inv));

[Torsten]

Like this ? "gamma" still has a wrong 2nd dimension and sigma_g has a wrong 1st dimension.

clc;
clear;
close all;
%% Data
% Define the parameters
qe = 1;
n = 5e16;
p = 0.525;
E_0 = 10e3;
v_s = 3.5e3;
epsilon_d = 50;
ue = 1;
m0 = 9.1094e-31;
m = 0.44 * m0;
epsilon_0 = 1;
h = 4.1356e-15;
hBar = h / (2 * pi);
b = 0.142e-9;
a = 3 * b / (2 * hBar);
tt = 2.6;
delta = 5.7 / 2;
G = sqrt(delta^2 + 5 * tt^2);
% Define the range for ω_p and B^(-1)
omega_p = linspace(0, 10, 150);
B_inv = linspace(0.04, 1, 150);
k_q = omega_p / v_s;
v = (-2 * tt^2 * a * sin(p * a) / G);
tau = ue * p / (qe * v);
tau_1 = tau / 2;
omega_c = qe * B_inv' * v / p;
gamma = 1 + 1i * omega_c * tau_1;
sigma_g = qe^2 * n * v * tau / p;
Six = p^2 / (2 * pi * hBar^2);
Rix = 2 * p / m / v;
b_g = 2 * pi * hBar^2 * epsilon_0 * v / (qe^2 * p);
beta_g = (k_q' .* v) ./ omega_c;
% Initialize the result matrix
j_x = zeros(numel(omega_p), numel(B_inv));
% Calculate sigma_xx
sigma_xx = zeros(numel(omega_p), numel(B_inv));
for r = -3:3
    J_r = besselj(r, beta_g);
    sigma_xx = sigma_xx + (2 * sigma_g' * J_r) ./ (1 - 1i * (omega_p' - r * omega_c) * tau);
end
% Calculate R_x
R_x = zeros(numel(omega_p), numel(B_inv));
for r = -3:3
    J_r = besselj(r, beta_g);
    R_x = R_x + (omega_c' ./ k_q) .* (J_r .^ 2) ./ (1 - 1i * (omega_p' - r * omega_c) * tau);
end
% Calculate g(omega_q, k_q)
g_omega_q_k_q = 1 + 1i / (epsilon_0 * (epsilon_d + 1)) * (sigma_xx ./ (v_s - R_x));
% Adjust dimensions of gamma and sigma_g
gamma = repmat(1 + 1i * omega_c * tau_1, numel(omega_p), 1);
sigma_g = repmat(qe^2 * n * v * tau / p, 1, numel(B_inv));
% Calculate the equation for each combination of omega_p and B^-1
size(gamma)

ans = 1×2

22500 1

size(sigma_g)

ans = 1×2

1 150

size(g_omega_q_k_q)

ans = 1×2

150 150

numel(B_inv)

ans = 150

numel(omega_p)

ans = 150

for i = 1:numel(omega_p)
    for j = 1:numel(B_inv)
        sum_term = 0;
        for r = -1:1
            J_r = besselj(r, beta_g(i));
            sum_term = sum_term + J_r / (1 - 1i * (omega_p(i) - r * omega_c(j)) * tau);
        end
        j_x(i, j) = (1 / (2 * qe * n * v(i))) * ((tau * omega_c(j) * p^2) / (8 * pi * (hBar^2))) * (abs((sigma_g(i, j) * E_0) / g_omega_q_k_q(i, j)))^2 ...
            * ((1 / beta_g(i)) / (1 + (omega_c(j)^2 * tau_1^2)))^2 * (sum_term) ...
            * (-1i * gamma(i, j)^2 * (sum_term + 1) * besselj(sum_term + 1, beta_g(i)) + 1i * conj(gamma(i, j))^2 * (sum_term - 1) * besselj(sum_term - 1, beta_g(i)));
    end
end

Index in position 2 exceeds array bounds. Index must not exceed 1.

Index in position 2 exceeds array bounds. Index must not exceed 1.
% Plot the results
[Omega_P, B_inv] = meshgrid(omega_p, B_inv);
surf(Omega_P, B_inv, j_x);
xlabel('\omega_p');
ylabel('B^{-1}');
zlabel('j_x');
title('Plot of j_x against \omega_p and B^{-1}');
title('Plot of j_x against \omega_p and B^{-1}');

[Samuel Suakye]

yeah..

[Samuel Suakye]

please any trick on that?

[Torsten]

There is no trick. You have to know what you want to calculate.

[Samuel Suakye]

I want to calculate for j_x and plot a 3D graph