Hi, I saw you code and executed MATLAB code in R2024a and I find it correct.
function [Is] = current()
[f,N,Nc,a1,a2,ra,k0,Z0,lambda ,gap_angle] = parameter();
% Angular positions φ_n
n=1:N;
phi_n = -gap_angle+ ((2*n - 1) * pi / N); % Compute φ_n
% Initialize matrices
A = zeros(4*N, 2*N); % Coefficient matrix (4N equations, 2N unknowns)
b = zeros(4*N, 1); % Right-hand side vector
% Fill the coefficient matrix A and vector b
for m = 1:N
for n = 1:N
% Compute angular difference and distances R1 and R2
phi_diff = phi_n(m) - phi_n(n);
R1 = sqrt(ra^2 + a1^2 - 2*ra*a1*cos(phi_diff));
R2 = sqrt(ra^2 + a2^2 - 2*ra*a2*cos(phi_diff));
% Bessel functions
H0_R1 = besselh(0, 2, k0*R1); % H_0^(2)(k0*R1)
H0_R2 = besselh(0, 2, k0*R2); % H_0^(2)(k0*R2)
H1_R1 = besselh(1, 2, k0*R1); % H_1^(2)(k0*R1)
H1_R2 = besselh(1, 2, k0*R2); % H_1^(2)(k0*R2)
% Equation (5): Rows 1:N
A(m, n)=(k0*Z0/4) * H0_R1;
% Equation (6): Rows N+1:2N
A(m+N, n+N) = H0_R2;
% Equation (7): Rows 2N+1:3N
A(m+2*N, n) = (k0*Z0/4) * H0_R1;
A(m+2*N, n+N) = -(k0*Z0/4) * H0_R2;
% Equation (8): Rows 3N+1:4N
A(m+3*N, n) = (1j*k0/4) * ((ra - a1*cos(phi_diff)) / R1) * H1_R1;
A(m+3*N, n+N) = -(1j*k0/4) * ((ra - a2*cos(phi_diff)) / R2) * H1_R2;
end
% Right-hand side vector b
b(m) =E_inc_z(ra, phi_n(m)); % Incident E_z for Eq (5)
b(m+N) = 0; % Zero for Eq (6)
b(m+2*N) =E_inc_z(ra, phi_n(m)); % Incident E_z for Eq (7)
b(m+3*N) = H_inc_phi(ra, phi_n(m)); % Incident H_phi for Eq (8)
end
% Solve the system using linsolve
x = linsolve(A, b);
% Extract unknowns w_n^(1) and w_n^(2)
w1 = x(1:N);
w2 = x(N+1:2*N);
Is = [w1; w2]
% Display results
%disp('w_n^(1):');
%disp(w1);
%disp('w_n^(2):');
%disp(w2);
end % added by Sam (Editor)
function [f, N, Nc, a1, a2, ra, k0, Z0, lambda, gap_angle] = parameter()
% Define or load the parameters here
f = 10e9; % Frequency (Hz)
N = 10; % Number of elements (example)
Nc = 4; % Some parameter Nc (example)
a1 = 0.5; % Radius 1 (example)
a2 = 0.6; % Radius 2 (example)
ra = 0.8; % Radius a (example)
k0 = 2 * pi * f / 3e8; % Wave number (example with speed of light)
Z0 = 377; % Impedance of free space (Ohms)
lambda = 3e8 / f; % Wavelength (meters)
gap_angle = pi/4; % Gap angle (example)
end
function E = E_inc_z(ra, phi)
% Define the incident electric field E_z (for example, a plane wave)
E = exp(-1j * ra * cos(phi)); % Example expression, modify as per your requirements
end
function H = H_inc_phi(ra, phi)
% Define the incident magnetic field H_phi (for example, a plane wave)
H = exp(-1j * ra * sin(phi)); % Example expression, modify as per your requirements
end
Is = current();
disp(Is);
Below are the key points for further explanation:
- the result is a 20×1 matrix, which matches the output size expecting for w1 and w2 combined (since w1 and w2 are each of size N×1N and N= 10 in setup).
- The first 10 values (from 1 to 10) are likely associated with the first unknowns (w1), and the second 10 values (from 11 to 20) are likely associated with the second unknowns (w2).
- For equation 5, the term (k0*Z0/4) * H0_R1 is correctly placed in the matrix A for the first set of equations. The right-hand side b(m) is set to E_inc_z(ra, phi_n(m)), which matches the equation.
- For equation 6, The term H0_R2 is correctly placed in the matrix A for the second set of equations. The right-hand side b(m+N) is set to 0, which matches the equation.
- For equation 7, The terms (k0*Z0/4) * H0_R1 and -(k0*Z0/4) * H0_R2 are correctly placed in the matrix A. The right-hand side b(m+2*N) is set to E_inc_z(ra, phi_n(m)), which matches the equation.
- For equation 8, The terms (1j*k0/4) * ((ra - a1*cos(phi_diff)) / R1) * H1_R1 and -(1j*k0/4) * ((ra - a2*cos(phi_diff)) / R2) * H1_R2 are correctly placed in the matrix A. The right-hand side b(m+3*N) is set to H_inc_phi(ra, phi_n(m)), which matches the equation.
