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
댓글 수: 1
george veropoulos
2024년 12월 30일
채택된 답변
% 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
% 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))
% Compute Escat for all points (vectorized)
Es_surf =sumEs_surf
m1 = max(abs(E_in_z(xs,ys))); % Vectorized Escat call
% Plot results
phi_d = phi*(180./pi) % Convert phi to degrees
E1 = Escat(xs,ys);
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?
댓글 수: 3
george veropoulos
2024년 12월 31일
이동: Torsten
2024년 12월 31일
Mr. Pavl M.
2025년 1월 2일
편집: Mr. Pavl M.
2025년 1월 2일
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
2025년 1월 4일
추가 답변 (2개)
george veropoulos
2025년 1월 4일
편집: george veropoulos
2025년 1월 4일
0 개 추천
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
댓글 수: 3
Mr. Pavl M.
2025년 1월 4일
편집: Mr. Pavl M.
2025년 1월 7일
Would mind to keep and keep on going objectives, objections, - real, legitimate, good job.
- For you to study and understand better the theory behind the fields, current and EM distributions around and inside natural figures... +
- 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?
- 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
2025년 1월 4일
편집: Walter Roberson
2025년 1월 5일
george veropoulos
2025년 1월 4일
카테고리
도움말 센터 및 File Exchange에서 Variables에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
