Problem in boundary condition

Hi

im trying to solve a proble in EM with method of MAS ( METHOD OF AUXILARY SOURCE)

I have a infiite cylinder ( PERFECT CONDUCTOR) and in incominf fileld strike the cylinder FROM THE Ein+Escat=0 in the surface of cylinder i solve th system and i fins the current Is ( the function is the currentMAS()) . After i find the Escat field form the epresion of MAS theort. After the i check the

sumEs_surf =abs(Etotal(xs,ys))

Etotal= Escat(x,y)+E_in_z(x,y)

xs=ra cos(phi)

ys=cos(phi)

rho = ra;

phi =2*pi*(0:Nc-1)/Nc% Create phi values

% Compute Cartesian coordinates (xa, ya)

xs = rho .* cos(phi); % x-coordinates

ys = rho .* sin(phi); % y-coordinates

the sum is not equal to zero !! the real part is zero but the imaginary not

function [Is] = currentMAS()
%UNTITLED Summary of this function goes here
%   Detailed explanation goes her
[~,N,Nc,a,ra,k0,Z0,~] = parameter();
% Preallocate arrays for efficiency
xc = ra * cos(2 * pi * (0:N-1)/N);   % Observation points (x-coordinates)
yc = ra* sin(2 * pi*(0:N-1)/N);   % Observation points (y-coordinates)
xa = a * cos(2 * pi * (0:Nc-1) / Nc); % Source points (x-coordinates)
ya = a * sin(2 * pi * (0:Nc-1) / Nc); % Source points (y-coordinates)

% Compute distance matrix R between source and observation points
[Xa, XC] = meshgrid(xa, xc); % Observation (columns), Source (rows)
[Ya, YC] = meshgrid(ya, yc);
Ra= sqrt((Xa - XC).^2 + (Ya - YC).^2);

% Compute matrix WI based on the Green's function
WI = (1 / 4) * k0 * Z0 .* besselh(0, 2, k0 .* Ra);

% Compute the source field at each source point
Source = arrayfun(@(xi, yi)E_in_z(xi,yi),xc,yc);
%Source = E_in_z(xc, 0); 

%Is = linsolve(WI, Source');
%lambda = 1e-1;  % Regularization factor
%Is = (WI' * WI + lambda * eye(size(WI, 2))) \ (WI' * Source');
Is =(pinv(WI)* Source') 
end

function z = Escat(x, y)
% Get the current values from the input parameters
[Is] = currentMAS()
[f,N,Nc,a,ra,k0,Z0,lambda] = parameter();
% Precompute coordinates of the points
xx = a * cos(2 * pi * (0:Nc-1) / Nc); % Length Nc
yy = a * sin(2 * pi * (0:Nc-1) / Nc); % Length Nc

% Initialize the S_UM for each (x, y) pair
E_scat = zeros(size(x)); % Result will be the same size as x and y

% Compute the distances and the contributions for each (x, y) pair
for jj = 1:Nc
    RR = sqrt((x - xx(jj)).^2 + (y - yy(jj)).^2); % Calculate distance for all (x, y)
    
    % Sum the contributions (note the broadcasting)
    E_scat =  E_scat - (1/4).*Is(jj) .* k0 * Z0 .* besselh(0, 2, k0 .* RR);
end

% Output the computed value
z =  E_scat;
end
function [f,N,Nc,a,ra,k0,Z0,lambda] = parameter()
%UNTITLED Summary of this function goes here
c0=3e8;
Z0=120.*pi;
ra=1;
a=0.85;
N=120;
Nc=120;
f=300e6;
lambda=c0./f;
k0=2*pi./lambda;
end

function y = Etotal(x,y)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
y=Escat(x,y)+E_in_z(x,y);
end
function   z=E_in_z(x,y)
%amplitude 
[f,N,Nc,a,ra,k0,Z0] = parameter();
E0=1;
z=E0.*exp(1i.*k0*x);
end

I Atach the mathematical description

tnank you

i wish a happy new year

[george veropoulos]

% Define phi and rho
clear all
[f,N,Nc,a,ra,k0,Z0,lambda] = parameter();
rho = ra;
phi =2*pi*(0:Nc-1)/Nc% Create phi values
% Compute Cartesian coordinates (xa, ya)
xs = rho .* cos(phi); % x-coordinates
ys = rho .* sin(phi); % y-coordinates
% Compute Escat for all points (vectorized)
Es_surf =abs(Etotal(xs,ys))
%max(abs(E_in_z(xs,ys))); % Vectorized Escat call
% Plot results
phi_d = phi*(180./pi)  % Convert phi to degrees
figure;
plot(phi_d, Es_surf, 'b', 'LineWidth', 2);
set(gca,'FontSize',19,'FontName','Times')
xlabel('$\phi$','Interpreter','latex')
ylabel('$BERROR$','Interpreter','latex' )
%legend('surface point','collocation point')
print('MAS_closed.jpg','-djpeg')

Answer 1

[Pavl M.]

% See your question initial body corrected:

%Hi,
%I'm trying to solve a problem in EM, (particullary Electric field?) distribution with
% method of MAS (METHOD OF AUXILARY SOURCE)
%I have  an infiite circular cylinder(PERFECT CONDUCTOR) and in incoming fileld strike
%the cylinder FROM THE Ein+Escat=0 in the surface of cylinder 
%I solve the system and I find the current Is (the function is the  currentMAS()). 
%After I find the Escat field form the expresion of MAS theory.
%After then I check the  sumEs_surf =abs(Etotal(xs,ys))  
%Etotal= Escat(x,y)+E_in_z(x,y)

%What if by integration, like in:
%https://physics.stackexchange.com/questions/113858/electric-field-at-surface-side-of-cylinder
%E = rho*r/(4*epsilon0)?
%E = E0*exp(j*k*b*cos(phi_i-phi_m))

%are you sure about the theoretic MAS formulas, besides implementation in
%script-code?
%What if N!=Nc?

%Where is transverse magnetic TM plane wave depicted?
%What with condition numbers?
%What are the placements of auxilliary sources(ASs) and collocation points(CPs)?
%Have you found the resonances points?
%Compute the distances and the contributions for each (x, y) pair
%Sum the contributions (note the broadcasting)
%What do you mean by the broadcasting?        
%How much field it absorbs, how much field it scatters(=radiates)?

%Who are they? What others?
%ver. 04 Jan 2025, To be continued...
clear; close all; clc;

function [f,N,Nc,a,ra,k0,Z0,lambda,E0] = parameter()
%

c0=3e8; %not to change, constant of light

Z0=120.*pi;
ra=1;
a=0.6; %ASs circle radius
N=40;
Nc=40;
E0 = ra/a;
f=500e6; %initially was 300e6

lambda=c0./f; %not to alter
k0=2*pi./lambda; %not to modify

hz = 2*pi*ra/(lambda*N); %normalized discretization step
phi1 = pi/N;
B = 10^7; %fixed MAS matrix condition number threshold
%arithmetic precision = 

end

function [Is] = currentMAS()
%

[~,N,Nc,a,ra,k0,Z0,~] = parameter();

    % Preallocate arrays for efficiency
    xc = ra * cos(2 * pi * (0:N-1)/N);   % Observation points (x-coordinates)
    yc = ra* sin(2 * pi*(0:N-1)/N);   % Observation points (y-coordinates)
    xa = a * cos(2 * pi * (0:Nc-1) / Nc); % Source points (x-coordinates)
    ya = a * sin(2 * pi * (0:Nc-1) / Nc); % Source points (y-coordinates)

    % Compute distance matrix R between source and observation points
    [Xa, XC] = meshgrid(xa, xc); % Observation (columns), Source (rows)
    [Ya, YC] = meshgrid(ya, yc);
    Ra= sqrt((Xa - XC).^2 + (Ya - YC).^2);

    % Compute matrix WI based on the Green's function
    WI1 = (1/4)*k0*Z0.*besselh(0,2,k0.*Ra);

    % Compute the source field at each source point
    Source1 = arrayfun(@(xi, yi)E_in_z(xi,yi),xc,yc);
    %Source = E_in_z(xc, 0);

    %Is = linsolve(WI, Source');
    %lambda = 1e-1;  % Regularization factor
    %Is = (WI' * WI + lambda * eye(size(WI, 2))) \ (WI' * Source');

    Is1 =(pinv(WI1)* Source1') %current

     % Initialize the Green's function matrix WI
    WI2 = zeros(N, Nc);
    
    %Interaction matrix:
    %Zmn = besselh(nu,K,Z,scale)

    % Compute WI matrix without meshgrid, Boundary condition formula,
    % are you sure about it's source and validity, (it looks okay in accord with the .pdf in my sight)?
    %Is1 it better than Is1, it looks more as per the .pdf than Is1?
    for i = 1:N
        %for j = 1:Nc
            % Compute the distance between observation point i and source point j
            R = sqrt((xa(i)-xc).^2+(ya(i)-yc).^2);
            % Fill the Green's function matrix
            WI2(i,:) = (1/4)*k0*Z0.*besselh(0,2,k0*R);
        %end
    end

    % Compute the source field at each observation point
    Source2 = zeros(N, 1);
    for i = 1:N
        Source2(i) = E_in_z(xc(i), yc(i));
    end

    % Solve the system of equations to find the current
    Is2 = linsolve(WI2, Source2)

    figure
    plot(abs(Is1))
    title('abs(Is1)')
    figure
    plot(abs(Is2))
    title('abs(Is2)')
    figure
    plot(phase(Is1))
    title('phase(Is1)')
    figure
    plot(phase(Is2))
    title('phase(Is2)')

    Is = (Is1+Is2)/2
    
    %ToDo:
    %Ra
    %Source1
    
    %figure
    %quiver(Xa,Ya,Ra,4)
    %figure
    %quiver(Ya,YC,Source1,3,"x")
    %figure
    %quiver(Xa,XC,Source2,0,"x")
    %figure
    %quiver(Ya,YC,Source2,3,"x")
    %figure
    %contour(Xa,XC,abs(Source1))   
    %figure
    %contour(Ya,YC,abs(Source1))  
    %figure
    %contour(Xa,XC,abs(Source2)) 
    %figure
    %contour(Ya,YC,abs(Source2))

end

function z = Escat(x, y)

    % Get the current values from the input parameters
    [Is] = currentMAS()

    figure
    plot(abs(Is))
    title('abs(Is)')

    figure
    plot(phase(Is))
    title('phase(Is)')

    [f,N,Nc,a,ra,k0,Z0,lambda,E0] = parameter();

    % Precompute coordinates of the points
    xx = a * cos(2 * pi * (0:Nc-1) / Nc); % Length Nc
    yy = a * sin(2 * pi * (0:Nc-1) / Nc); % Length Nc

    %[Xa, XC] = meshgrid(xa, xx); % Observation (columns), Source (rows)
    %[Ya, YC] = meshgrid(ya, yy);

    % Initialize the S_UM for each (x, y) pair
    E_scat = zeros(size(x)); % Result will be the same size as x and y

    %Compute the distances and the contributions for each (x, y) pair
    %As per formula of Boundary Conditions in surface of cyllinder, are you
    %sure the formula in .pdf is correct? What about potentials?
    for jj = 1:Nc
        RR = sqrt((x - xx(jj)).^2 + (y - yy(jj)).^2); % Calculate distance for all (x, y)
        % Sum the contributions (note the broadcasting)
        E_scat =  E_scat-(1/4).*Is(jj).*k0*Z0.*besselh(0,2,k0.*RR);
    end

    % Output the computed value
    z =  E_scat;
end

function y = Etotal(x,y)
%
y=Escat(x,y)+E_in_z(x,y);
end

function   z=E_in_z(x,y)
%amplitude
[f,N,Nc,a,ra,k0,Z0,E0] = parameter();
z=E0.*exp(1i.*k0*x);
end

[f,N,Nc,a,ra,k0,Z0,lambda] = parameter();

rho = ra;
phi =2*pi*(0:Nc-1)/Nc% Create phi values

phi = 1×40

0 0.1571 0.3142 0.4712 0.6283 0.7854 0.9425 1.0996 1.2566 1.4137 1.5708 1.7279 1.8850 2.0420 2.1991 2.3562 2.5133 2.6704 2.8274 2.9845 3.1416 3.2987 3.4558 3.6128 3.7699 3.9270 4.0841 4.2412 4.3982 4.5553

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Compute Cartesian coordinates (xa, ya)
xs = rho .* cos(phi); % x-coordinates
ys = rho .* sin(phi); % y-coordinates

%[Xa, XC] = meshgrid(xa, xs); % Observation (columns), Source (rows)
%[Ya, YC] = meshgrid(ya, ys);

sumEs_surf = abs(Etotal(xs,ys))

Is1 = 
   0.0005 + 0.0007i
   0.0007 + 0.0002i
   0.0009 + 0.0005i
  -0.0003 + 0.0006i
  -0.0006 - 0.0013i
   0.0023 - 0.0001i
  -0.0034 + 0.0034i
  -0.0025 - 0.0081i
   0.0103 + 0.0043i
  -0.0137 + 0.0071i
   0.0024 - 0.0175i
   0.0134 + 0.0122i
  -0.0146 + 0.0046i
   0.0031 - 0.0110i
   0.0065 + 0.0049i
  -0.0029 + 0.0028i
  -0.0011 - 0.0017i
   0.0009 - 0.0010i
   0.0003 + 0.0002i
  -0.0005 - 0.0001i
  -0.0007 - 0.0001i
  -0.0005 - 0.0001i
   0.0003 + 0.0002i
   0.0009 - 0.0010i
  -0.0011 - 0.0017i
  -0.0029 + 0.0028i
   0.0065 + 0.0049i
   0.0031 - 0.0110i
  -0.0146 + 0.0046i
   0.0134 + 0.0122i
   0.0024 - 0.0175i
  -0.0137 + 0.0071i
   0.0103 + 0.0043i
  -0.0025 - 0.0081i
  -0.0034 + 0.0034i
   0.0023 - 0.0001i
  -0.0006 - 0.0013i
  -0.0003 + 0.0006i
   0.0009 + 0.0005i
   0.0007 + 0.0002i

Is2 = 
  -0.0007 - 0.0001i
  -0.0005 - 0.0001i
   0.0003 + 0.0002i
   0.0009 - 0.0010i
  -0.0011 - 0.0017i
  -0.0029 + 0.0028i
   0.0065 + 0.0049i
   0.0031 - 0.0110i
  -0.0146 + 0.0046i
   0.0134 + 0.0122i
   0.0024 - 0.0175i
  -0.0137 + 0.0071i
   0.0103 + 0.0043i
  -0.0025 - 0.0081i
  -0.0034 + 0.0034i
   0.0023 - 0.0001i
  -0.0006 - 0.0013i
  -0.0003 + 0.0006i
   0.0009 + 0.0005i
   0.0007 + 0.0002i
   0.0005 + 0.0007i
   0.0007 + 0.0002i
   0.0009 + 0.0005i
  -0.0003 + 0.0006i
  -0.0006 - 0.0013i
   0.0023 - 0.0001i
  -0.0034 + 0.0034i
  -0.0025 - 0.0081i
   0.0103 + 0.0043i
  -0.0137 + 0.0071i
   0.0024 - 0.0175i
   0.0134 + 0.0122i
  -0.0146 + 0.0046i
   0.0031 - 0.0110i
   0.0065 + 0.0049i
  -0.0029 + 0.0028i
  -0.0011 - 0.0017i
   0.0009 - 0.0010i
   0.0003 + 0.0002i
  -0.0005 - 0.0001i

Is = 
  -0.0001 + 0.0003i
   0.0001 + 0.0001i
   0.0006 + 0.0003i
   0.0003 - 0.0002i
  -0.0008 - 0.0015i
  -0.0003 + 0.0014i
   0.0016 + 0.0042i
   0.0003 - 0.0096i
  -0.0022 + 0.0045i
  -0.0001 + 0.0096i
   0.0024 - 0.0175i
  -0.0001 + 0.0096i
  -0.0022 + 0.0045i
   0.0003 - 0.0096i
   0.0016 + 0.0042i
  -0.0003 + 0.0014i
  -0.0008 - 0.0015i
   0.0003 - 0.0002i
   0.0006 + 0.0003i
   0.0001 + 0.0001i
  -0.0001 + 0.0003i
   0.0001 + 0.0001i
   0.0006 + 0.0003i
   0.0003 - 0.0002i
  -0.0008 - 0.0015i
  -0.0003 + 0.0014i
   0.0016 + 0.0042i
   0.0003 - 0.0096i
  -0.0022 + 0.0045i
  -0.0001 + 0.0096i
   0.0024 - 0.0175i
  -0.0001 + 0.0096i
  -0.0022 + 0.0045i
   0.0003 - 0.0096i
   0.0016 + 0.0042i
  -0.0003 + 0.0014i
  -0.0008 - 0.0015i
   0.0003 - 0.0002i
   0.0006 + 0.0003i
   0.0001 + 0.0001i

Is = 
  -0.0001 + 0.0003i
   0.0001 + 0.0001i
   0.0006 + 0.0003i
   0.0003 - 0.0002i
  -0.0008 - 0.0015i
  -0.0003 + 0.0014i
   0.0016 + 0.0042i
   0.0003 - 0.0096i
  -0.0022 + 0.0045i
  -0.0001 + 0.0096i
   0.0024 - 0.0175i
  -0.0001 + 0.0096i
  -0.0022 + 0.0045i
   0.0003 - 0.0096i
   0.0016 + 0.0042i
  -0.0003 + 0.0014i
  -0.0008 - 0.0015i
   0.0003 - 0.0002i
   0.0006 + 0.0003i
   0.0001 + 0.0001i
  -0.0001 + 0.0003i
   0.0001 + 0.0001i
   0.0006 + 0.0003i
   0.0003 - 0.0002i
  -0.0008 - 0.0015i
  -0.0003 + 0.0014i
   0.0016 + 0.0042i
   0.0003 - 0.0096i
  -0.0022 + 0.0045i
  -0.0001 + 0.0096i
   0.0024 - 0.0175i
  -0.0001 + 0.0096i
  -0.0022 + 0.0045i
   0.0003 - 0.0096i
   0.0016 + 0.0042i
  -0.0003 + 0.0014i
  -0.0008 - 0.0015i
   0.0003 - 0.0002i
   0.0006 + 0.0003i
   0.0001 + 0.0001i

sumEs_surf = 1×40

0.5196 0.4767 0.3057 0.0564 0.4890 0.5405 0.0765 0.5995 0.0566 0.5986 0.0000 0.5986 0.0566 0.5995 0.0765 0.5405 0.4890 0.0564 0.3057 0.4767 0.5196 0.4767 0.3057 0.0564 0.4890 0.5405 0.0765 0.5995 0.0566 0.5986

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Compute Escat for all points (vectorized)
Es_surf =sumEs_surf

Es_surf = 1×40

0.5196 0.4767 0.3057 0.0564 0.4890 0.5405 0.0765 0.5995 0.0566 0.5986 0.0000 0.5986 0.0566 0.5995 0.0765 0.5405 0.4890 0.0564 0.3057 0.4767 0.5196 0.4767 0.3057 0.0564 0.4890 0.5405 0.0765 0.5995 0.0566 0.5986

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

m1 = max(abs(E_in_z(xs,ys))); % Vectorized Escat call

% Plot results

phi_d = phi*(180./pi)  % Convert phi to degrees

phi_d = 1×40

0 9.0000 18.0000 27.0000 36.0000 45.0000 54.0000 63.0000 72.0000 81.0000 90.0000 99.0000 108.0000 117.0000 126.0000 135.0000 144.0000 153.0000 162.0000 171.0000 180.0000 189.0000 198.0000 207.0000 216.0000 225.0000 234.0000 243.0000 252.0000 261.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

E1 = Escat(xs,ys);

Is1 = 
   0.0005 + 0.0007i
   0.0007 + 0.0002i
   0.0009 + 0.0005i
  -0.0003 + 0.0006i
  -0.0006 - 0.0013i
   0.0023 - 0.0001i
  -0.0034 + 0.0034i
  -0.0025 - 0.0081i
   0.0103 + 0.0043i
  -0.0137 + 0.0071i
   0.0024 - 0.0175i
   0.0134 + 0.0122i
  -0.0146 + 0.0046i
   0.0031 - 0.0110i
   0.0065 + 0.0049i
  -0.0029 + 0.0028i
  -0.0011 - 0.0017i
   0.0009 - 0.0010i
   0.0003 + 0.0002i
  -0.0005 - 0.0001i
  -0.0007 - 0.0001i
  -0.0005 - 0.0001i
   0.0003 + 0.0002i
   0.0009 - 0.0010i
  -0.0011 - 0.0017i
  -0.0029 + 0.0028i
   0.0065 + 0.0049i
   0.0031 - 0.0110i
  -0.0146 + 0.0046i
   0.0134 + 0.0122i
   0.0024 - 0.0175i
  -0.0137 + 0.0071i
   0.0103 + 0.0043i
  -0.0025 - 0.0081i
  -0.0034 + 0.0034i
   0.0023 - 0.0001i
  -0.0006 - 0.0013i
  -0.0003 + 0.0006i
   0.0009 + 0.0005i
   0.0007 + 0.0002i

Is2 = 
  -0.0007 - 0.0001i
  -0.0005 - 0.0001i
   0.0003 + 0.0002i
   0.0009 - 0.0010i
  -0.0011 - 0.0017i
  -0.0029 + 0.0028i
   0.0065 + 0.0049i
   0.0031 - 0.0110i
  -0.0146 + 0.0046i
   0.0134 + 0.0122i
   0.0024 - 0.0175i
  -0.0137 + 0.0071i
   0.0103 + 0.0043i
  -0.0025 - 0.0081i
  -0.0034 + 0.0034i
   0.0023 - 0.0001i
  -0.0006 - 0.0013i
  -0.0003 + 0.0006i
   0.0009 + 0.0005i
   0.0007 + 0.0002i
   0.0005 + 0.0007i
   0.0007 + 0.0002i
   0.0009 + 0.0005i
  -0.0003 + 0.0006i
  -0.0006 - 0.0013i
   0.0023 - 0.0001i
  -0.0034 + 0.0034i
  -0.0025 - 0.0081i
   0.0103 + 0.0043i
  -0.0137 + 0.0071i
   0.0024 - 0.0175i
   0.0134 + 0.0122i
  -0.0146 + 0.0046i
   0.0031 - 0.0110i
   0.0065 + 0.0049i
  -0.0029 + 0.0028i
  -0.0011 - 0.0017i
   0.0009 - 0.0010i
   0.0003 + 0.0002i
  -0.0005 - 0.0001i

Is = 
  -0.0001 + 0.0003i
   0.0001 + 0.0001i
   0.0006 + 0.0003i
   0.0003 - 0.0002i
  -0.0008 - 0.0015i
  -0.0003 + 0.0014i
   0.0016 + 0.0042i
   0.0003 - 0.0096i
  -0.0022 + 0.0045i
  -0.0001 + 0.0096i
   0.0024 - 0.0175i
  -0.0001 + 0.0096i
  -0.0022 + 0.0045i
   0.0003 - 0.0096i
   0.0016 + 0.0042i
  -0.0003 + 0.0014i
  -0.0008 - 0.0015i
   0.0003 - 0.0002i
   0.0006 + 0.0003i
   0.0001 + 0.0001i
  -0.0001 + 0.0003i
   0.0001 + 0.0001i
   0.0006 + 0.0003i
   0.0003 - 0.0002i
  -0.0008 - 0.0015i
  -0.0003 + 0.0014i
   0.0016 + 0.0042i
   0.0003 - 0.0096i
  -0.0022 + 0.0045i
  -0.0001 + 0.0096i
   0.0024 - 0.0175i
  -0.0001 + 0.0096i
  -0.0022 + 0.0045i
   0.0003 - 0.0096i
   0.0016 + 0.0042i
  -0.0003 + 0.0014i
  -0.0008 - 0.0015i
   0.0003 - 0.0002i
   0.0006 + 0.0003i
   0.0001 + 0.0001i

Is = 
  -0.0001 + 0.0003i
   0.0001 + 0.0001i
   0.0006 + 0.0003i
   0.0003 - 0.0002i
  -0.0008 - 0.0015i
  -0.0003 + 0.0014i
   0.0016 + 0.0042i
   0.0003 - 0.0096i
  -0.0022 + 0.0045i
  -0.0001 + 0.0096i
   0.0024 - 0.0175i
  -0.0001 + 0.0096i
  -0.0022 + 0.0045i
   0.0003 - 0.0096i
   0.0016 + 0.0042i
  -0.0003 + 0.0014i
  -0.0008 - 0.0015i
   0.0003 - 0.0002i
   0.0006 + 0.0003i
   0.0001 + 0.0001i
  -0.0001 + 0.0003i
   0.0001 + 0.0001i
   0.0006 + 0.0003i
   0.0003 - 0.0002i
  -0.0008 - 0.0015i
  -0.0003 + 0.0014i
   0.0016 + 0.0042i
   0.0003 - 0.0096i
  -0.0022 + 0.0045i
  -0.0001 + 0.0096i
   0.0024 - 0.0175i
  -0.0001 + 0.0096i
  -0.0022 + 0.0045i
   0.0003 - 0.0096i
   0.0016 + 0.0042i
  -0.0003 + 0.0014i
  -0.0008 - 0.0015i
   0.0003 - 0.0002i
   0.0006 + 0.0003i
   0.0001 + 0.0001i

E2 = E_in_z(xs,ys);

figure
plot(phi_d,abs(E1))
title('Absolute Scat. Field')

figure
plot(phi_d,phase(E1))
title('Phase Scat. Field')

figure
quiver(real(E1),imag(E1),3)
title('Directions Scat. Field')

figure
plot(phi_d,abs(E2))
title('Absolute Vertical Field')

figure
plot(phi_d,phase(E2))
title('Phase Vertical. Field')

figure
quiver(real(E2),imag(E2),3)
title('Directions vertical field')

figure;
plot(phi_d, Es_surf, 'b', 'LineWidth', 2);
hold on
plot(m1)
set(gca,'FontSize',19,'FontName','Times')
xlabel('$\phi$','Interpreter','latex')
ylabel('$Berror$','Interpreter','latex' )
%legend('surface point','collocation point')
print('MAS_closed.jpg','-djpeg')

My this doesn't solve it yet, while sparks more interest to the field/domain of physics and more benign attention to the question and shows that I don't really knowledgeable in MAS methode, which is true in order not to mess, not to undermine and show experience in other fields and to show that jack-of-all-trades is always maximum a jack-of just-single 1 trade. Not to gamble for sure. Of course why, for what we do not pretend nor try to produce under jack-of-all-trades or over-scattered umbrella,it is yet slightly gathered(limited), let it focus more, I just wanted to help to enhance, make it better in confidence.

Novel contribution to this question is that to further develop quiver() and contour() plots for better perception of the field.

What with MoM and Gauss, SIBC methods?

You noticed low(very small) values of Berror, which is a difference(error) value between E1-E2, so logically it may be deemed allright to be small?

It is also not that paid assignement and so let you see who others will do it better on so or other terms. You are OK.

As an additional sense, to realize me with as much luggage, things to come to settle around Athens in secret place or other not previous bad, but to move forward, where there are people are needed and economically benign not to charge you much as to understand how difficult poority from outside caused to us lawfully to overcome previous difficulties would be fine so it would be also fairly well to check applications of the MAS methode and fields utilization for transfering machinery not only for the cyllinder likely used in antennas.

I have completed this 2 loops to 1 loop re-scripting:

% Initialize the Green's function matrix WI

WI2 = zeros(N, Nc);

% Compute WI matrix without meshgrid, Boundary condition formula,

% are you sure about it's source and validity, (it looks okay in accord with the .pdf in my sight)?

%Is1 it better than Is1, it looks more as per the .pdf than Is1?

for i = 1:N

% Compute the distance between observation point i and source point j

R = sqrt((xa(i)-xc).^2+(ya(i)-yc).^2);

% Fill the Green's function matrix

WI2(i,:) = (1/4)*k0*Z0.*besselh(0,2,k0*R);

end

What is epsilon0 (permiability) in this case?

Additional pecular interesting point of this question is that E field is incoming (explicit) and not by internal cyllinder

charge distribution?

Complex systems are very sensitive to initial conditions, script implements that 1/4 * bessel(...) formula 1 to 1 as from .pdf, was it cheked for validity? What in case of other N, Nc ( I tested 2 cases of different N,Nc by re-running modified scripts and I have not achieved gain in E),finite, limited cyllinder?

[george veropoulos]

The problem is that the total field I. Thdkthis. Point are not equal to zero

[Pavl M.]

Dear George Veropoulos,

OK. What is sure, there are many healthy points and flows (see phase plot) and a little bit of ill, sick with due,(E1,E2, quiver usage to amend).

You wrote: "The problem is that the total field I. Thdkthis. Point are not equal to zero"

To continue more clear, what do you mean by total field this Es_surf = Etotal(xs,ys) ?

What do you mean by I, which current?

What do you mean by Thd?

Who, why, for what wants to suppress the field around that figure(cyllinder)? What is more correct is to amplify better to amplify good, not parazitic, make stronger and more beutiful shaped utile field and so associated with EM field E current ->.

What do you mean by point are not equal to zero?

[george veropoulos]

Dear Pavl

yes this i mean Es_surf =abs(Etotal(xs, ys)) because the cylinder is perfect conductor this sum must ne equal to zero.

rho = ra;
%phi = 2*pi*(0:Nc-1)/Nc; % Create phi values
phi =0:2*pi/(2*N):2*pi
% Compute Cartesian coordinates (xs, ys)
xs = rho .* cos(phi); % x-coordinates
ys = rho .* sin(phi); % y-coordinates

When i change the msgrid with double for loop this quantity goes to zero!

George

Answer 2

[george veropoulos]

the problem was in [Xa, XC] = meshgrid(xa, xc); % Observation (columns), Source (rows)

[Ya, YC] = meshgrid(ya, yc);

may be give wrong results I use a double For loop and the abs( Etotl) is very smail

[Pavl M.]

Would mind to keep and keep on going objectives, objections, - real, legitimate, good job.

1. For you to study and understand better the theory behind the fields, current and EM distributions around and inside natural figures... +
2. Where was abs(Etotl) very small? In my answer Es_surf =sumEs_surf is not very small, min = 0, max = 2, see plot Berror vs phi.Where to trans-script/transform the double for loop to loopless matrix based operations, maybe can do. Can you click on my the an answer accepted?
3. I don't know many distributions. Let who others else skillfull to solve it let's see..

For 2D or 3D field plot?

What is other, there, is it so simple? Here it is not rudimentary. Gamble not. Bet not of course. Kindly see the recent my update to the answer here - magnitude improved, not small contribution. Is it elementry there? On it quiver built-in command execution will be better than on my so far answer of course. (so far the quiver results of the filed of course yet are very not satisfiable.)

So far in code there were xs, ys, xx, yy. What are the xa, xc, ya, yc? - Found, clarified.

[george veropoulos]

function [Is] = currentMAS()
% Load parameters
[~, N, Nc, a, ra, k0, Z0, ~] = parameter();
% Preallocate arrays for observation and source points
xc = ra * cos(2*pi*(0:N-1)/N);
yc = ra * sin(2*pi*(0:N-1)/N);
xa =a * cos(2*pi*(0:Nc-1)/Nc);
ya =a * sin(2*pi*(0:Nc-1)/Nc);
% Initialize the Green's function matrix WI
WI = zeros(N, Nc);

% Compute WI matrix without meshgrid
for i = 1:N
    for j = 1:Nc
        % Compute the distance between observation point i and source point j
        R = sqrt((xa(i) - xc(j))^2 + (ya(i) - yc(j))^2);
        % Fill the Green's function matrix
        WI(i, j) = (1/4) * k0 * Z0 * besselh(0, 2, k0 * R);
    end
end

% Compute the source field at each observation point
Source = zeros(N, 1);
for i = 1:N
    Source(i) = E_in_z(xc(i), yc(i));
end

% Solve the system of equations to find the current
Is = linsolve(WI, Source);
end

[george veropoulos]

the sum Einc(xs,ys)+Escat(xs,ys)=0 overt the cylinder surface must be equal to zero xs=ra*co(phi) ys=ra*sin(phi) ra is the radius of cylinder the Einc is the incoming EM field and Escat is the scattering EM field

Answer 3

[george veropoulos]

thank you !! very useful answers