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

2 views (last 30 days)
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()
Function 'currentMoM' has already been declared within this scope.
% 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)
  2 Comments
Les Beckham
Les Beckham on 7 Jan 2025
What is your question?
Also, please format your code as code. To do this, edit your question and select all of the code, and then click the left button in the CODE section of the editor toolstrip.
george veropoulos
george veropoulos on 7 Jan 2025
Moved: Torsten on 7 Jan 2025
i have two function with name currentMOIM use the Method of Moment. One use the pulse basis and the second the trinagle basis function. In any case i find the scattering field Escat(r, phi) the resluts are agree when the icnoming fileld is z=E0.*exp(-j.*k0.*ra.*cos(phi)) for the MOM pulse and z=E0.*exp(-j.*k0.*ra.*cos(phi)) for the MoM triangle... This is starnge... ?

Sign in to comment.

Answers (0)

Categories

Find more on Matrices and Arrays in Help Center and File Exchange

Products


Release

R2024b

Community Treasure Hunt

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

Start Hunting!