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()

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)

[Les Beckham]

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]

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