Difference results from MOM method with ttρiangle basis function and pulse

Hi
i find the scattering field in cylinder and i use the MOM method with the triangle basis and pulse. The |Es
| are the same when then incoming field are Einc= E0*exp(-j*k0*r*cos(φ)) in MOM pulse and Einc= E0*exp(j*k0*r*cos(φ)) in MOM triangle.
This the pulse current
function [Is] = currentMoM()
% Load parameters
[f, N, ra, k0, Z0, lambda] = parameter();
% Constants
gamma_const = 1.781;
dfpm = 2*pi/N;
Phi0 = (0:N-1) * dfpm;
% Precompute pairwise phase differences
deltaPhi = zeros(N, N); % Initialize deltaPhi matrix
for i = 1:N
for j = 1:N
deltaPhi(i, j) = Phi0(i) – Phi0(j); % Calculate pairwise phase differences
end
end
% Compute distance matrix R (N x N matrix)
R = zeros(N, N); % Initialize R matrix
for i = 1:N
for j = 1:N
R(i, j) = ra * sqrt(2 – 2 * cos(deltaPhi(i, j))); % Compute distance matrix R
end
end
% Coefficients
coeif1 = (Z0 * k0 / 4) * ra * dfpm;
coeif2 = (Z0 / 2) * sin(k0 * ra * dfpm/ 2);
coeif3 = (gamma_const * k0 * ra * dfpm) / 4;
% Compute lmn matrix
lmn = zeros(N, N); % Initialize lmn matrix
for i = 1:N
for j = 1:N
if i==j
lmn(i, j) = coeif1 * (1 – 1i * (2 / pi) * (log(coeif3) – 1)); % Diagonal elements
else
lmn(i, j) = coeif2 * besselh(0, 2, k0 * R(i, j)); % Compute each element of lmn
end
end
zmn = lmn;
end
% Handle diagonal elements separately

% Set zmn equal to lmn

% Compute gm vector using a for loop
gm = zeros(N, 1); % Initialize gm vector
for i = 1:N
gm(i) = E_in_z(Phi0(i)); % Calculate each element of gm using E_in_z function
end
% Solve for currents using linsolve
%Is = linsolve(zmn, gm);
Is = zmn gm;
%Is=linsolve(zmn, gm)
end
and the triangke
function [Is] = currentMoM()
% Parameter initialization
[~, N, ra, k0, Z0, ~] = parameter();
gamma_const = 1.781;
Eps = 1e-9;
dftm =2.*pi./N;

% Precompute constants
B = (4./(Z0 * k0));
A = (gamma_const*k0*ra)./2;

% Initialize matrices and vectors
lmn = zeros(N, N); % (N x N) matrix for lmn
gm = zeros(1, N); % (1 x N) vector for gm
zmn = zeros(N, N); % (N x N) matrix for zmn

% Compute lmn and zmn matrices using nested for-loops
for index_i = 1:N
for index_j = 1:N
% Precompute integration bounds
xi1_i = (index_i – 1) * dftm;
xi2_i = xi1_i + dftm;
xj1_j = (index_j – 1) * dftm;
xj2_j = xj1_j + dftm;
%log(A) +log(abs(x – y) +Eps))
if index_i == index_j
% Diagonal elements
funa = @(x, y) triangle_basisn(x, index_i) .* triangle_basisn(y, index_j) .*ra .* (1 – 1i * (2./pi).*(log(A) +log(abs(x – y) +Eps))) ;
lmn(index_i, index_j) = integral2(funa, xi1_i, xi2_i, xj1_j, xj2_j,’RelTol’, 1e-6, ‘AbsTol’, 1e-6,’Method’,’iterated’);
else
% Off-diagonal elements
funb = @(x, y)triangle_basisn(x, index_i) .* triangle_basisn(y, index_j) .*ra .* besselh(0, 2, k0.*ra.*sqrt(abs(2 – 2 * cos(x – y))) );
lmn(index_i, index_j) = integral2(funb, xi1_i, xi2_i, xj1_j, xj2_j, ‘RelTol’, 1e-6, ‘AbsTol’, 1e-6);
end

% Assign to zmn
zmn(index_i, index_j) = lmn(index_i, index_j);
end

end

% Compute gm vector using a single for-loop
for index_i = 1:N
xi1 = (index_i – 1) * dftm;
xi2 = xi1 + dftm;
func = @(x)B.* triangle_basisn(x, index_i) .* (E_in_z(x));
gm(index_i) = integral(func, xi1, xi2, ‘RelTol’, 1e-6, ‘AbsTol’, 1e-6);
end

% Solve the linear system
epsilon_reg = 1e-10; % Regularization parameter (can be adjusted)

zmn = zmn + epsilon_reg * eye(N); % Add small value to the diagonal
% Solve the linear system
Is = linsolve(zmn, gm’);

end

the scattering field
function z = Escat(r, phi)
% Load parameters
[f, N, ra, k0, Z0, lambda] = parameter();
[Is] = currentMoM(); % Retrieve the current distribution

% Preallocate for final result
FINAL = zeros(size(r)); % Initialize FINAL to the same size as r

dfpulse = 2 * pi / N;
Phi0 = (0:N-1) * dfpulse; % Angle steps (in radians)
coif = (k0 * Z0 / 4); % Coefficient to scale the results
% Create a column vector of Phi0 values
x1 = Phi0(:);
x2 = x1 + dfpulse; % Integration limits (x1, x2)

% Loop to compute the scattered field for each combination of rho and phi
for jj = 1:N
% Ensure phi and r are compatible in size for broadcasting

r_squared = r.^2; % Compute r squared
ra_squared = ra.^2;
% Compute the integrand term

% Define the integrand with the properly dimensioned term
integrand = @(y)besselh(0, 2, k0 * sqrt(r_squared + ra_squared – 2*r.*ra .*cos(phi-y)));

% Perform the numerical integration for each jj
FINAL = FINAL -coif * Is(jj).* ra .* integral(integrand, x1(jj), x2(jj),’ArrayValued’,true);
end

% Return the result
z = FINAL;
end
function z = Escat(r,phi)
[~,N,ra,k0,Z0,~] =parameter();
Is=currentMoM();
%Phi0=zeros(N);
FINAL=0;
A=(k0.*Z0./4);
dftm=2.*pi./N;
%index=1:N
%Phi0=(index-2)*dftm;
Phi0=(0:N-1).*dftm;
x1 = Phi0(:); % Create a column vector of x1 values
x2 = x1 + dftm; % x2 values depend on x1
for jj=1:N
F=@(y)besselh(0,2,k0*sqrt(r.^2 + ra.^2-2.*r.*ra.*cos(phi-y))).* triangle_basisn(y,jj);

% Perform the integral
%pts = optimset(‘TolFun’,1e-6); % Adjust tolerance
%integral_result = integral(integrand, lower_limit, upper_limit, ‘ArrayValued’, true);

integral_result=integral(F,x1(jj),x2(jj),’RelTol’, 1e-6, ‘AbsTol’, 1e-6,’ArrayValued’,false);
% Accumulate the result
FINAL = FINAL-A.*Is(jj).*ra.*integral_result;
end
z=FINAL;
end
function [f,N,ra,k0,Z0,lambda] = parameter()
%UNTITLED Summary of this function goes here
c0=3e8;
Z0=120.*pi;
ra=1;
N=60;
f=300e6;
lambda=c0./f;
k0=2*pi./lambda;
end
and teh triangle basis function
function z = triangle_basisn(phi, k)
% Get number of segments (N) and compute segment width (dftm)
[~, N, ~, ~, ~, ~] = parameter();
df = 2*pi/ N; % Segment width

% Compute the left and right boundaries for the triangle basis
left_bound = (k – 2) * df;
right_bound = (k – 1) * df;
center_bound = k * df;

% Check if phi lies within the support of the triangular basis function
if phi >= left_bound & phi <= right_bound
% Left part of the triangle
z = (phi – left_bound) / df;
elseif phi > right_bound & phi <= center_bound
% Right part of the triangle
z = (center_bound – phi) / df;
else
% Outside the support of the triangle
z = 0;
end
end

function z=E_in_z(phi)
%amplitude
[f,N,ra,k0,Z0] = parameter();
E0=1;
%z=E0.*exp(j.*k0.*ra.*cos(phi)+j.*k0.*ra.*sin(phi));
z=E0.*exp(+j.*k0.*ra.*cos(phi));
end
for the triangle

function z=E_in_z(phi)
%amplitude
[f,N,ra,k0,Z0] = parameter();
E0=1;
%z=E0.*exp(j.*k0.*ra.*cos(phi)+j.*k0.*ra.*sin(phi));
z=E0.*exp(-j.*k0.*ra.*cos(phi));
end
for the pulse

thank you
George

the main script is

clear all
clc
[f,N,ra,k0,Z0,lambda] = parameter();
%ph_i=pi/2;
rho=6*lambda ;
phi=-pi:pi/180:pi;
Es=zeros(size(phi));
phi_d=zeros(size(phi));
%phi=0:pi/200:pi/3
tic
parfor jj=1:length(phi)
Es(jj)=abs(Escat(rho,phi(jj)));
phi_d(jj)=phi(jj).*(180/pi);
end

plot(phi_d,Es,’b–‘,’LineWidth’,2)Hi
i find the scattering field in cylinder and i use the MOM method with the triangle basis and pulse. The |Es
| are the same when then incoming field are Einc= E0*exp(-j*k0*r*cos(φ)) in MOM pulse and Einc= E0*exp(j*k0*r*cos(φ)) in MOM triangle.
This the pulse current
function [Is] = currentMoM()
% Load parameters
[f, N, ra, k0, Z0, lambda] = parameter();
% Constants
gamma_const = 1.781;
dfpm = 2*pi/N;
Phi0 = (0:N-1) * dfpm;
% Precompute pairwise phase differences
deltaPhi = zeros(N, N); % Initialize deltaPhi matrix
for i = 1:N
for j = 1:N
deltaPhi(i, j) = Phi0(i) – Phi0(j); % Calculate pairwise phase differences
end
end
% Compute distance matrix R (N x N matrix)
R = zeros(N, N); % Initialize R matrix
for i = 1:N
for j = 1:N
R(i, j) = ra * sqrt(2 – 2 * cos(deltaPhi(i, j))); % Compute distance matrix R
end
end
% Coefficients
coeif1 = (Z0 * k0 / 4) * ra * dfpm;
coeif2 = (Z0 / 2) * sin(k0 * ra * dfpm/ 2);
coeif3 = (gamma_const * k0 * ra * dfpm) / 4;
% Compute lmn matrix
lmn = zeros(N, N); % Initialize lmn matrix
for i = 1:N
for j = 1:N
if i==j
lmn(i, j) = coeif1 * (1 – 1i * (2 / pi) * (log(coeif3) – 1)); % Diagonal elements
else
lmn(i, j) = coeif2 * besselh(0, 2, k0 * R(i, j)); % Compute each element of lmn
end
end
zmn = lmn;
end
% Handle diagonal elements separately